Fracture toughness of hierarchical self-similar honeycombs

International Journal of Solids and Structures 152–153 (2018) 151–160

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Fracture toughness of hierarchical self-similar honeycombs Michael Ryvkin∗, Raz Shraga The Iby and Aladar Fleischman Faculty of Engineering, School of Mechanical Engineering, Tel Aviv University, Ramat Aviv 69978, Israel

a r t i c l e

i n f o

Article history: Received 30 January 2018 Revised 11 May 2018 Available online 23 August 2018 Keywords: Self-similar hierarchy Periodic honeycomb Brittle fracture Discrete Fourier transform

a b s t r a c t The influence of self-similar hierarchy on the brittle fracture behavior of a two-dimensional honeycomb is investigated. The honeycomb walls are modeled as Bernoulli-Euler beam elements rigidly connected at the nodal points, and zero, first and second order hierarchies are addressed. The stress state in the vicinity of a semi-infinite crack is determined, and the fracture toughness is calculated in terms of tensile strength of parent material. The growth of hierarchical order leads to increase in the number of nodal degrees of freedom to be taken into account, and the computational cost of the calculations is reduced by employing a novel analysis method based on the discrete Fourier transform. An increase in hierarchical order enhances the fracture toughness after optimization of geometrical parameters at fixed relative density. This effect becomes more pronounced for honeycombs of lower relative density. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction An efficient way to improve the mechanical performance of solid structural element is to replace the homogeneous material with a material with microstructure. The iteration of this procedure leads to the hierarchical materials possessing remarkable strength, stiffness and toughness. There are a numerous examples of implementation of this design strategy both in nature (e.g., Chen and Pugno (2013); Gibson (2012); Lakes (1993); Gao (2006)) and engineering practice (e.g., Sen and Buehler, 2011; Carpinteri and Pugno, 2008; Lakes, 1993; Taylor et al., 2011; Kooistra et al., 2007; Fan et al., 2008; Mousanezhad et al., 2015a; Mousanezhad et al., 2015b). Hierarchical materials can be broken down into the two main groups: dense composite materials and cellular ones (Lakes, 1993). Representatives of the first group are the high-rank sequentially laminated composites that have stiffness approaching the theoretical bounds (see Milton (2002) and references there). The bone microstructure represents an example of a high toughness dense hierarchical bio-composite (Currey, 1984). Cellular hierarchical bio-composites (e.g., Gibson, 2012) serve as a source of inspiration for many studies stimulated by recent technology advances that provide an attractive opportunity to manufacture new microarchitectured materials (e.g., Meza et al., 2015; Arabnejad et al., 2016). There are two main ways to introduce hierarchy in a periodic structure. The first way, which is illustrated in (Fig. 1a), is to con-

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (M. Ryvkin).

https://doi.org/10.1016/j.ijsolstr.2018.06.022 0020-7683/© 2018 Elsevier Ltd. All rights reserved.

sider parent solid material of structural elements as an effective one representing higher order microstructure (e.g., Vigliotti and Pasini, 2013; Zhu and Wang, 2013). Another way (Fig. 1b) is to introduce scaled periodic pattern of the structure at the nodal points where several elements are meet and in this case the original and the scaled patterns may be of the same length scale (e.g., Ajdari et al., 2012; Li et al., 2017). Using this approach, it is possible to obtain fractal-like self-similar honeycombs which are the subject of the present investigation. The mechanical behavior of self-similar hexagonal honeycombs was studied intensively in several recent works. Their stiffness and strength were investigated in Ajdari et al. (2012); Haghpanah et al. (2013) and Oftader et al. (2014), and linear and non-linear elastic behavior of a spider-web honeycomb was considered in Mousanezhad et al. (2015a). The influence of the selfsimilar geometry of two-dimensional phononic crystals on their wave propagation behavior was examined in Mousanezhad et al. (2015). However, contrary to non-hierarchical honeycombs (see, for example, Fleck and Qiu, 2007; Kucherov and Ryvkin, 2014), the brittle fracture behavior of hierarchical ones has not been investigated yet and it is the subject of the present paper. The paper is organized as follows. In the next section, the microstructure of the honeycombs considered is described, and the general scheme of the fracture toughness evaluation is outlined. Section 3 gives the analysis procedure in detail, and it includes also verification examples. Section 4 presents the results of the parametric study. In the final Section 5, concluding remarks are drawn, and the components of the homogeneous K−field employed in the analysis are specified in the Appendix.

152

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

placement step i are defined by referring to the initial honeycomb

γi =

bi b0

and

ηi =

ti , t0

(1)

where ti and bi are the thickness and length of beam elements, respectively. The geometrical constraint of non-overlapping for different hexagons of the nth order in a self-similar hierarchical honeycomb has the form (Haghpanah et al., 2013)

γi ≤ γ j , i < j , n 

γi ≤

i=1 n 

1 , 2

γi ≤ γ j , 0 < j < n .

(2)

i= j+1

Fig. 1. Two ways of introducing structural hierarchy.

Similar to the original H0 honeycomb, all these microstructures are isotropic on a macroscale thanks to the elastic symmetry, and the Young modulus and Poisson ratio of homogenized material are denoted as E∗ and ν ∗ , respectively. Another essential honeycomb material characteristic is its relative density defined as area fraction of the solid parent material. In the present work, we carry out calculations for honeycombs with zero, first and second order hierarchies. Their relative densities are given as

2 t0 3 b0

ρ0 = √

(3) 2 t0 3 b0

ρ1 = (4η1 γ1 − 2γ1 + 1 ) √

(4) 2 t0 . 3 b0

ρ2 = 2[(12η2 − 4η1 − 2 )γ2 + 2γ1 (2η1 − 1 ) + 1] √

(5)

2.2. The fracture toughness evaluation

Fig. 2. Self-similar hexagonal hierarchical honeycombs.

2. Problem formulation 2.1. Self-similar hierarchical honeycombs Consider two-dimensional hexagonal honeycomb composed of beam elements rigidly connected at nodal points (Fig. 2a). The honeycomb has unit out-of-plane thickness, and a case of the plane stress deformation is addressed. The beams satisfy the BernoulliEuler assumptions and can undergo bending and axial deformations. The length and thickness of the beams are denoted as b0 and t0 respectively, and a periodic pattern of the honeycomb, which is referred as zero-order one H0 , is shown in Fig. 2b. Note that it is possible to use the more primitive unit cell as used in the studies of elastic (e.g., Vigliotti and Pasini, 2012) and fracture properties (Lipperman et al., 2007) of this honeycomb. However, the pattern considered preserves the structure of its boundary with six beam elements approaching it for honeycombs of the higher hierarchical orders, and it is convenient for the development of a uniform analysis procedure. Higher order self-similar hierarchical honeycombs are obtained from the basic zero-order one by sequential replacement of each node where three elements join with a scaled hexagon. The periodic cells for the first and second order hierarchies H1 and H2 are shown in Fig. 2c,d. The non-dimensional parameters for each re-

The subject of the present investigation is the fracture toughness of the self-similar hierarchical honeycombs. Assuming that brittle fracture and crack advance occur when the maximum local tensile stress in some beam in the crack-tip vicinity attains the tensile strength of the parent material σ fs , it is possible to derive the fracture toughness of a honeycomb KC theoretically by a numerical experiment (see Schmidt and Fleck, 2001; Fleck and Qiu, 2007). At the outer boundary of the honeycomb, the conditions applied correspond to the K-field for a traction-less semi-infinite crack in the homogeneous material with the effective properties E∗ , ν ∗ . After evaluation of the stress state of the honeycomb, the maximal tensile stress in the crack tip vicinity σmax is determined by comparison for all beams the expression

σmax =

Mt Nax + . 2I t

(6)

Here, M is an absolute value of the bending moment at a beam extremity, t is the beam thickness, I is moment of inertia of beam’s cross-section and Nax is the axial force in the beam element. Then, in view of the problem linearity, one obtains

KC = K

σfs , σmax

(7)

where K is the amplitude of the homogeneous K− field solution, which is presented in the Appendix. The analysis scheme is illustrated in Fig. 3. Note that the rectangular analysis domain f must include a sufficiently large number of repetitive cells to provide true values of stresses in the zone of interest at the crack-tip vicinity. Consequently, the condition

nf =

Lf 1, l

(8)

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

153

defined by two translations and one rotation, their positive directions are shown in Fig. 4b. The honeycomb equilibrium is provided by the equilibrium equations of the cells

Kuk = Fk ,

k = {k1 , k2 },

(9)

ki = −Ni , −Ni + 1, ., Ni − 1; i = 1, 2.

Fig. 3. The fracture toughness evaluation in a honeycomb.

where Lf is the length of the rectangle side, and l is the characteristic honeycomb periodicity parameter must be fulfilled. The actual value of the ratio nf is to be determined by numerical experiments. In the next section we present the evaluation of the stress field in rectangular honeycomb domain with the K-field boundary conditions. 3. Analysis method 3.1. Problem in local coordinates Consider a rectangular honeycomb domain −Li ≤ Xi ≤ Li , i = 1, 2 with a crack (Fig. 4). For definiteness, hereafter in this section, Mode I crack in H1 hierarchical honeycomb is addressed, the analysis steps remain the same for any hierarchical order and Mode II crack as well. Each honeycomb node where several beam elements are rigidly connected has two translational and one rotational degrees of freedom, and the stress state of the honeycomb is completely defined by the nodal displacements vector. The size of this vector, i.e., the total number of degrees of freedom to be included in the analysis increases rapidly with the order of hierarchy. Therefore, the Representative Cell Method, RCM (see Ryvkin and Nuller, 1997; Ryvkin, 2008), based on the discrete Fourier transform is employed. This method allows for a significant reduction in the computational cost of the analysis; it was employed successfully for the solution of similar problems for finite periodic composites (Ryvkin and Hadar, 2015; Ryvkin et al., 2017) and beam lattices (Cherkaev and Ryvkin, 2017). Therefore, the solution procedure is presented in brief. The honeycomb is viewed as an assemblage of bonded identical cells indicated in Fig. 4 by the dashed lines. Since the cell boundaries cross the beams, the additional nodes are introduced at the intersection points. The cells are numbered by the use of vector index k = {k1 , k2 }, where ki = −Ni , −Ni + 1, . . . , Ni − 1; i = 1, 2. This figure illustrates the cell numbering for the case N1 = N2 = 2. Note, that the actual values of these parameters used in the calculations are much higher. The crack faces are tractions-free, and the sought stress state is caused by the K-field conditions applied at the rectangle boundaries. In accordance with the RCM approach, the problem is considered in local coordinate systems (x1 , x2 ), which are introduced in all cells identically at their centers (see Fig. 4b). Consequently, the numbering of 30 cell nodes is also uniform, and the numbering index of boundary nodes is shown in the figure. Thus, the identity of each honeycomb node is determined by index k, which defines the cell location and index s, 1 ≤ s ≤ 30, that specifies the node number within the cell. A node displacement vector uks ≡ {u1 , u2 , uθ }ks is

Here, uk = {u1 , u2 , ., u30 }k is the vector of nodes displacements having 30 × 3 = 90 entries, K is the 90 × 90 stiffness matrix of a cell, which is independent of k, and Fk is a vector of external forces applied to the nodes Fk = {F1 , F2 , ., F30 }k , Fs = {F1(s) , F2(s) , Ms }. The positive directions of the force components are the same as of the generalized displacements shown in Fig. 4b. These forces express the interaction between the neighboring cells, consequently, only the forces applied to boundary nodes 1, 2, 3, 4, 5 and 6 are non-zero. Together with the nodal displacements, they define the boundary values vectors k V1+ = {u3 , u4 , F3 , F4 }k ,

V1k− = {u1 , u2 , −F1 , −F2 }k , k V2+ = {u6 , F6 }k ,

V2k− = {u5 , −F5 }k ,

k = {k1 , k2 }.

(10)

The minus sign is introduced for the force terms at the boundaries xi = −li , i = 1, 2 for convenience in order to unify the formulation of continuity conditions at the interfaces between the cells. These conditions are read as follows k1 ,k2 V1+ − V1k−1 +1,k2 = 0,

(11)

k1 = −N1 , −N1 + 1, ., N1 − 2; k2 = −N2 , −N2 + 1, ., N2 − 1, k1 ,k2 V2+ − V2k−1 ,k2 +1 = Gkcr1 δk2 ,−1 ,

k1 = −N1 , −N1 + 1, ., N1 − 1; k2 = −N2 , −N2 + 1, ., N2 − 2,

(12)

where the first and the second equations express the conditions at the interfaces parallel to the X2 and X1 axis, respectively. The nonzero right-hand side of the second equation, including Kronecker delta symbol δ ij , yields from the discontinuity of the displacements at the crack line between the cells with k2 = −1 and k2 = 0. For the considered Mode I crack cr k1 Gkcr1 = {gcr , u , gF }

{gcru }k1 = {ucr1 , ucr2 , ucrθ }k1 , {gcrF }k1 = {F1cr , F2cr , Mcr }k1 .

(13)

{ucr }k1 2

Here, jumps in the crack opening displacements and rotak1 are unknown a priory and should be defined by tion angles {ucr } θ an iterative procedure, and the jumps in the horizontal displacements and force values vanish {ucr }k1 = 0, {gcr }k1 = {0, 0, 0}k1 . 1 F The conditions at the rectangular domain boundaries are applied by specifying the jumps in displacements and force components for the opposite nodes N1 −1,k2 1 ,k2 V1+ − V1−N = Gk12 , for the boundaries X1 = ±L1 , − k1 ,N2 −1 V2+ − V2k−1 ,−N2 = Gk21 , for the boundaries X2 = ±L2 ,

ki = −Ni , −N1 , ., Ni − 1; i = 1, 2,

(14)

where

Gk3i−i = {gu(3−i ) , gF(3−i ) }ki , i = 1, 2,

{gu(1) }k2 = {uK3 − uK1 , uK4 − uK2 }k2 , {gF(1) }k2 = {FK3 − FK1 , FK4 − FK2 }k2 , {gu(2) }k1 = {uK6 − uK5 }k1 , {gF(2) }k1 = {FK6 − FK5 }k1 . Contrary to specification of the K-field by prescribing the displacements at the remote boundary (e.g., Fleck and Qiu, 2007) here both kinematic and force parameters are used, the uniqueness of such formulation was demonstrated by Ryvkin and Hadar (2016). Note, that the employed connection between the first and the last cell in

154

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

Fig. 4. Periodic domain with cell numbering (a) and the representative cell (b).

each direction may be viewed as a closure of a periodic strip into a cyclic circle. This fact was mentioned in Karpov et al. (2002) where application of the discrete Fourier transform for the analysis of periodic beam structures was considered. The magnitudes of the boundary displacements and forces are derived from the K-field solution, which is presented in the Appendix. For example, for node number 6 located at cell with k1 = j at the upper boundary of the domain

{uK6 } j = {uK1 , uK2 , uKθ }(X1(6) , X2(6) ) ,

{FK6 } j = {F1K(6) , F2K(6) , M(K6) } j ,

{F2K(6) } j = {

} =

M(K6) j

  

2l1 ( j+2 ) 2l1 ( j+1 ) 2l1 ( j+2 ) 2l1 ( j+1 ) 2l1 ( j+2 ) 2l1 ( j+1 )

K σ12 dX1 ,

K 22 dX1 .

(16)

The problem (9), (11)–(14) formulated for the honeycomb consisting of 4N1 N2 cells is reduced to the problem for the single representative cell by application of the double discrete Fourier transform N 2 −1

f (k1 ,k2 ) ei(ϕ1 k1 +ϕ2 k2 ) ,

k1 =−N1 k2 =−N2

where

ϕr =

π qr Nr

,

qr = −Nr , −Nr + 1, ., Nr − 1 , r = 1, 2. (17)

After some manipulation (details of this procedure can be found in Ryvkin et al. (2017) where the general approach to solution of periodic problems of the considered type is presented), one obtains

KuF (ϕ1 , ϕ2 ) = FF (ϕ1 , ϕ2 ), VF1+ (ϕ1 , ϕ2 ) − γ1−1 VF1− (ϕ1 , ϕ2 ) = γ1N1 −1 GF12 (ϕ2 ), VF2+ (ϕ1 , ϕ2 ) − γ2−1 VF2− (ϕ1 , ϕ2 ) = γ2N2 −1 GF21 (ϕ1 ) − γ2−1 GFcr1 (ϕ1 ),

r = 1, 2,

(18)

f ( k j ) e iϕ j k j ,

k j =−N j

where ϕ j =

πqj Nj

, q j = −N j , −N j + 1, ., N j − 1 ; j = 1, 2. (19)

The unknowns in the structural mechanics problem obtained for the representative cell are the complex valued transforms of nodal displacements uF and forces acting at cell boundaries FF . After the solution of the problem for 4N1 N2 combinations of the discrete transform parameters (ϕ 1 , ϕ 2 specified in (17)) it is possible to derive the actual displacements and forces by the inverse transform formula. For example,

where

N 1 −1



{F}k1 ,k2 =

3.2. The representative cell problem

f F ( ϕ1 , ϕ2 ) =

N j −1

f Fj ( ϕ j ) =

K σ22 dX1 ,

(X1(6) − l1 /2 )σ

γr = exp(iϕr ),

(15)

where X1(6 ) = 2l1 ( j + 1 ) + l1 /2 and X2(6 ) = L2 are the global coordinates of the node. The force components are the resultants of the K-field stresses acting at the cell boundary 2l1 ( j + 1 ) ≤ X1 ≤ 2l1 ( j + 2 ) , X2 = L2

{F1K(6) } j =

where

and superscript Fj , j = 1, 2 denotes the one-dimensional Fourier transform

ϕr =

N 1 −1 1 4N1 N2

π qr Nr

N 2 −1

FF (ϕ1 , ϕ2 )e−i(ϕ1 k1 +ϕ2 k2 ) ,

q1 =−N1 q2 =−N2

,

qr = −Nr , −Nr + 1, ., Nr − 1 , r = 1, 2. (20)

The stress state in all honeycomb beams is determined then by a standard procedure using local stiffness matrices. Note, that when at least one of the transform parameters ϕ r vanishes, then γr = 1, and, for the adopted jumps formulation of boundary conditions, the solution of the system (18) becomes nonunique due to the possibility of a rigid body motion in the corresponding direction. Consequently, in these cases, a degree of freedom in this direction is fixed. The obtained solution depends upon the crack opening and jumps in the rotation angles at the crack faces defined by vectors k1 {ucr }k1 and {ucr 2 θ } (see (13)). These vectors are to be adjusted to provide zero values for the corresponding interface forces

{F2cr }k1 = 0 , {Mcr }k1 = 0 , k1 = −N1 , −N1 + 1, ., −1

{F1cr }k1

(21)

(third force component vanishes due to the problem symmetry). The iterations conducted in order to fulfill this condition are organized as following. At the first step, the jump vectors k1 are derived from the K-field solution for the ho{ucr }k1 , {ucr 2 θ} mogenized material presented in the Appendix. Then the representative cell problem is solved, and the interface forces at the

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

crack line are determined by (17). The deviation from the condition (21) is analyzed for each force component of the obtained solution, and a change in the prescribed jumps in the displacement values is introduced in order to diminish the force values. Then the solution for the representative cell with the new jumps is obtained, and the procedure is repeated until (21) is fulfilled with the required accuracy. 3.3. Effective properties The remote displacements in the K-field solution, which are necessary for formulation of the boundary conditions, depend upon the effective elastic properties of the honeycomb material. The considered self-similar honeycombs possess threefold symmetry; consequently, their in-plane isotropic elastic behavior can be defined by two effective moduli E∗ and ν ∗ . Their values for the first and second order self-similar honeycombs being considered were determined by Ajdari et al. (2012). The authors derived the analytic expressions for the case when just the bending compliance of honeycomb walls was considered, and obtained numerical results using the finite elements method when bending, axial and shear deformations were taken into account. In the present analysis, we address bending and axial deformation of beam elements, consequently, the effective properties were established for this type of elastic behavior. The effective moduli are found from the periodic solution for the bi-axial tension of the rectangular honeycomb domain without the crack. This solution is obtained by the use of the advocated analysis method (9)–(17). To this end, the jump of displacements u2 in X2 − direction defined in (14) is fixed, and all other displacement jumps as well as the force jumps and the jumps at the crack line (13) are set to zero. The average normal stress in X1 − (X2 −) direction is calculated by dividing the corresponding boundary force resultant for the nodes 1 and 2 (node 6) by the cell size 2l2 (2l1 ), see Fig. 4. The average deformations are defined from the known displacements jumps, and the sought values of elastic moduli are derived from the Hooke’s law. The obtained results are in complete agreement with the data reported by Ajdari et al. (2012). 3.4. Model parameters The theoretical fracture toughness problem is considered for a semi-infinite crack in an unbounded honeycomb, and the finite domain formulation illustrated in Fig. 3 must include the domain size analysis. The application of the K-field conditions at the rectangle boundaries at the finite distance from the crack tip leads to a deviation of the calculated stress values from the theoretical ones. The deviation may be significant within a boundary layer; properties of the layers of this type that appears at the boundaries of finite size lattices were investigated in Fleck and Qiu (2007), and Phani and Fleck (2008). The analysis goal is to establish such a domain size that the boundary layer does not reach the zone of interest near the crack tip where the maximal stresses defining the fracture toughness are developed. The numerical tests were carried out for the H0 honeycombs with different relative densities. The value of nf defined in (8) was gradually increased and it was found that for n f = 70 the influence of the domain size becomes negligibly small. Consequently, we use the domains with N1 = N2 = 35 including 4900 repetitive cells for the fracture toughness analysis. Since the number of nodes in a hierarchical honeycomb of order n scales with 3n a high number of degrees of freedom for the considered model makes the direct numerical simulation rather computationally expensive. The increasing complexity of the microstructure of hierarchical honeycombs is illustrated in Fig. 5, where the data for the representative cells and total degrees of freedom numbers in the model are presented. This

155

observation explains the choice of the discrete Fourier transform based representative cell method as the analysis tool. The method allows to address the structural mechanics problem for the representative cell only and, consequently, to reduce the computational cost of the analysis. For example, in the present problem for the second order hierarchy one has to resolve the system for the representative cell with 270 degrees of freedom 70 × 35 = 2450 times in order to get the stress state at each iteration step. The advantage of this procedure is that it requires less computer resources than the direct finite element modeling of the entire domain including about 1.2 million degrees of freedom. On the other hand, there is a disadvantage of the suggested method in view of iterative approaching of the solution, which is absent in the direct modeling. The balance between these two tendencies may be different for different problems, but it plausible to assume that for the higher hierarchical orders the suggested method will be preferred. 3.5. Verification of the approach The suggested methodology allows to determine the fracture toughness for the honeycombs of arbitrary hierarchical order. For the basic hexagonal layout with H0 hierarchy, it is possible to verify the employed approach by a comparison with the results known from the literature. The dependence of the normalized Mode I fracture toughness upon the relative density related parameter t0 /b0 (see (3)) is presented in Fig. 6. A very good agreement of our results with the results of Fleck and Qiu (2007), and Lipperman et al. (2007) is observed. 4. Results 4.1. First order hierarchy In this section the influence of the microstructure geometry of the hierarchical honeycombs on their brittle fracture behavior is examined. We begin with the first order hierarchy. It is assumed that for each layout all the beam elements have equal thickness, i.e., η1 = 1, t0 = t1 , and dependence on the parameter γ 1 is considered. In accordance with (2) it varies for this hierarchy in the region 0 ≤ γ 1 ≤ 1/2. It is assumed further that the relative density of honeycomb is fixed, and this constraint is met by varying the beams’ thickness, which is defined from (4)

t1 =

√ 3ρ1 b0 . 2(2γ1 + 1 )

(22)

Thus, the material volume in a unit area is fixed, and the uniform thickness of beams is defined by the total length of the layout net in it. The central question in the analysis is a relation between the calculated fracture toughness and the toughness of the basic H0 honeycomb with γ1 = 0 (see Fig. 1a). For this honeycomb the ratio t0 /b0 = 0.1, which corresponds to the relative density ρ = 0.115, is chosen. The results for the Mode I and Mode II fracture toughness are presented in Fig. 7a, they are normalized using the corresponding values for H0 hierarchy. Note that the absolute values of fracture toughness of hierarchical honeycomb for Mode I are found to be significantly higher than for Mode II. This phenomenon is well known from the investigations of non-hierarchical (zero order) honeycombs (Fleck and Qiu, 2007; Lipperman et al., 2007). It is seen, that in the small region in the vicinity of the optimal value of parameter γ 1 the fracture toughness of the H1 −hierarchical honeycomb is higher than of the basic H0 one. The increase, however, is relatively small: 2.1% for Mode I and 3.6% for Mode II. The optimal γ 1 -values are 0.35 and 0.39 for the Mode I and Mode II, respectively. A possible explanation for this fact is as follows. For these values, the original H0 honeycomb composed of the beams with length b0 is replaced by the hierarchical one

156

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

Fig. 5. Computational model parameters for hierarchical honeycombs.

Fig. 6. Mode I fracture toughness of the basic hexagonal honeycomb vs. beam element thickness.

where all the beams (including both introduced hexagons and the remaining parts of the H0 layout between them, see Fig. 2c) have approximately the same length b0 /3. Consequently, the maximal size of the longest beam of the honeycomb approaches its minimal value, and, therefore, the mechanical behavior of this bending dominated honeycomb is improved. This feature affects the stiffness of the honeycomb also rather than just the fracture toughness. In fact, as it is seen in Fig. 7b where the graph for effective Young modulus E∗ is presented, the above values of γ 1 are very close to the value corresponding to the maximal honeycomb stiffness. Note that, contrary to the fracture toughness, the effective Young modulus for the H1 −hierarchical honeycomb is higher than for the basic H0 one for the most values of γ 1 .

Fig. 7. Dependence of the normalized fracture toughness (a) and the effective elastic Young modulus (b) for first order hierarchical honeycomb upon the structural parameter γ1 = b1 /b0 . The honeycomb relative density ρ = 0.115.

4.2. Second order hierarchy For the second order hierarchical honeycombs, we also limit our consideration to the case of equal thickness of the beams η1 = η2 = 1. Then each layout is defined by the two parameters γ 1 and γ 2 , which satisfy the constraints

γ2 ≤ γ1 ,

γ1 + γ2 ≤ 1/2

(23)

in accordance with non-overlapping conditions (2). As previously, the uniform thickness of beams is adjusted in order to meet the stipulation of fixed relative density of the honeycomb. The results for the case ρ = 0.115 are presented in Fig. 8a,b for both

fracture modes. As in the H1 case, the fracture toughness values for the basic H0 −honeycomb (γ1 = γ2 = 0) are used for the normalization. The data obtained for the Mode I show that if we consider a honeycomb with H1 hierarchy (the region γ2 = 0 in Fig. 8a) and introduce H2 hierarchy by defining γ 2 increment, then the fracture toughness decreases everywhere besides a small region in the vicinity of the H1 -optimal value γ1 = 0.35. This value changes slightly for the optimal parameter combination (γ1 )opt = 0.32, (γ2 )opt = 0.125 that provides the maximal fracture toughness

(KICH2 /KICH0 )opt = 1.092. Interestingly, all the beam elements for the

optimal configuration have approximately the same length as in

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

157

Fig. 8. Dependence of the Mode I (a), Mode II (b) fracture toughness and effective elastic Young modulus (c) of the second order hierarchical honeycomb upon the structural parameters. The honeycomb relative density ρ = 0.115.

first order hierarchy case. The optimal layouts for the first and second order hierarchical honeycombs are presented schematically at the inset in Fig. 11a. For the Mode II fracture mode the results are similar: there is an improvement and it is more substantial than for the H1 hierarH H chy, (KIIC2 /KIIC0 )opt = 1.2. The optimal parameter values are close to the values for the Mode I case: (γ1 )opt = 0.33 and (γ2 )opt = 0.15. Note, that as for the H1 -case, the Mode II fracture toughness is found to be more sensitive to the hierarchical change of the microstructure. It is interesting to note that, as for the first order hierarchy case, the dependencies of the fracture toughness and effective Young modulus (see Fig. 8a,b,c) of hierarchical honeycomb on the microstructure parameters are quite similar here. However, the improvement of the Young modulus for the optimal parameter combination is considerably more significant than for the fracture toughness.

in the previous examples, ρ = 0.0289, and calculate Mode I fracture toughness for the first and second hierarchy orders. The results are presented in Fig. 9a,b. It is seen that the hierarchy effect on the toughness increase is markedly higher for the considered density than for ρ = 0.115. In particular, for the second order hierarchy, it is about 40 percent. The optimal values of layout parameters are slightly affected, (γ1 )opt = 0.35 for the first order hierarchy and (γ1 )opt = 0.33, (γ2 )opt = 0.15 for the second one. The results for the Mode I fracture toughness improvement are summarized in Fig. 10. Another relative density aspect to be addressed is its influence on the fracture toughness of hierarchical honeycomb for the given layout. For two-dimensional honeycombs, relative density is proportional to the beams thickness t0 /b0 (see (3)–(5)), and, therefore, the graphs presented in the verifying subsection for the basic hexagonal H0 -honeycomb illustrate the dependence of toughness upon the density as well. The type of this dependence is well known

4.3. Influence of the relative density

KIC = Aρ d ,

The influence of the relative density of cellular materials on their mechanical properties is of particular interest (Gibson and Ashby, 1997) Therefore, it is natural to examine the influence of the relative density variation on the considerations presented above. To this end, we consider honeycombs with the density lower than

where A is a parameter independent of ρ and d = 2. For nonhexagonal honeycombs power d may be different (see Fleck and Qiu, 2007) and it is of interest to establish its value for the hexagonal honeycombs of higher hierarchical orders being considered. The results obtained are shown in Fig. 11a,b.

(24)

158

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

Fig. 9. Dependence of the Mode I fracture toughness of first order (a) and second order (b) hierarchical honeycombs upon the structural parameters defined in (1). The honeycomb relative density ρ = 0.0289. Fig. 11. Dependence of the fracture toughness upon the relative density for the first and second order optimal layouts(a), and for the several layouts for the first order honeycomb(b).

crostructure geometries, all hierarchical honeycombs remain bending dominated. For this type of honeycombs, the quadratic dependence of the fracture toughness upon the relative density follows from the general considerations (Maiti et al., 1984). 5. Concluding remarks

Fig. 10. Impact of hierarchy upon the honeycomb fracture toughness for optimal combinations of the structural parameters.

In the Fig. 11a the dependencies are presented for the honeycombs of the first and second order hierarchies with geometric parameters close to the optimal values determined in the previous sections, while in Fig. 11b the first order hierarchy with several combinations γ 1 , γ 2 is addressed. In addition, the results for the basic H0 -honeycomb are shown for the reference in both figures. It is readily seen that, in all cases, the index d of the power law is the same, i.e., formula (24) with d = 2 remains valid. This is not surprising and points to the fact that, in spite of different mi-

High order hierarchical self-similar hexagonal honeycombs can be produced by sequential replacement of each node where three edges join by a hexagon of a smaller scale. The study of their brittle fracture behavior is an expensive computational task due to a large number of degrees of freedom to be incorporated into the analysis. This obstacle is circumvented by employing the nonstandard analysis method based on the discrete Fourier transform, which is found to be a convenient analysis tool for the modeling of periodic hierarchical systems. This approach allows one to replace the analysis of sizeable periodic domain, including many repetitive cells, by the multiple analysis of the single cell. The parametric study of the fracture toughness of the first and second order hierarchical honeycombs with uniform beam element thickness for both Mode I and Mode II fracture modes is performed. It is found that, for the case of fixed relative density, introducing hierarchy with arbitrarily chosen geometric parameters usually leads to the toughness decrease. However, for the optimal

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

parameter combination, the fracture toughness of a honeycomb with a higher hierarchical order is enhanced with reference to the lower order one. The positive effect of hierarchical layout becomes more pronounced for honeycombs with lower relative densities. For example, the comparison of the Mode I fracture toughness of the second order hierarchical honeycomb and the non-hierarchical (zero order) one shows that for relative density ρ = 0.115 the increase is 5.4%, while for ρ = 0.0289 it is 39%. Interestingly, the optimal values of the parameters that give the maximal fracture toughness and the maximal effective Young modulus are approximately the same. Note that the above values are very close to the values reported by Ajdari et al. (2012), which maximize the Young modulus of hierarchical honeycombs with different relative densities. The analysis of the magnitudes of the fracture toughness enhancement for different hierarchical orders demonstrates that the difference between the second and first order hierarchies is more significant than the difference between hierarchies of the first and zero orders. This result, however, probably cannot be extrapolated to the higher order hierarchies, and a kind of saturation for the maximal possible toughness value is expected. This saturation phenomenon regarding the stiffness of self-similar hexagonal honeycombs was observed by Oftadeh et al. (2014), the inefficiency of higher than second hierarchical orders was mentioned by Meza et al. (2015) who studied strength and stiffness of different hierarchical honeycombs. It should be noted that for the higher order hierarchies the microstructure elements may become so small that it is required to include the surface elasticity effects into consideration (see Zhu et al., 2012; Zhu et al., 2014). Such investigation seems to be an essential next step in the topic. Another possible research direction for the future studies of the influence of hierarchy on the fracture behavior is extension of the design space by consideration of heterogeneous honeycombs with non-uniform thickness of beam elements (parameters ηi defined in (1) become the design variables). This will enable further increase of the positive hierarchy effect for the optimal parameter combinations. Employing such heterogeneous layouts was found to be an effective tool for design of stiff (Taylor et al., 2011) and fault-tolerant structures (Cherkaev and Ryvkin, 2017), increase of buckling stress for heterogeneous honeycombs was mentioned by Jang and Kyriakides (2015). The obtained dependence of the fracture toughness for hierarchical honeycombs with the given layout upon their relative density is characterized by the same power scaling law (24) as for a hexagonal non-hierarchical honeycomb. This result shows that, in spite of complicated microstructure, the considered honeycombs remain bending dominated. The length of elements in a high order hierarchical honeycomb is small with respect to the length of elements in the basic zero order one, however, the length scale parameters defining the periodicity of the zero order and high order hierarchies are the same. Therefore, since the fracture toughness scales with square root of the characteristic length parameter, it can be higher for the hierarchical honeycomb than for the zero order one built of identical small elements.

Acknowledgement The authors gratefully acknowledge support by Israel Science Foundation, Grant No. 1494/16.

Appendix Homogeneous eigensolution (K− field) for Mode I crack in isotropic material with Young modulus E and Poisson ratio ν has

the form (e.g., Lawn, 1993)

σ (X1 , X2 ) = K 11

K σ22 (X1 , X2 ) = K σ12 (X1 , X2 ) =

uK1 (X1 , X2 ) = uK2 (X1 , X2 ) =

159





K θ θ 3θ cos 1 − sin sin , 2 2 2 2 π r 1/2   K θ θ 3θ cos 1 + sin sin , 2 2 2 2 π r 1/2 K θ θ 3θ cos sin cos , 2 2  2 2 π r 1/2    1/2 K r θ 3θ , (1 + ν ) (2κ − 1 ) cos − cos 2E 2π 2 2     K r 1/2 θ 3θ , (1 + ν ) (2κ + 1 ) cos − sin 2E 2π 2 2

where

κ=

3−ν for the considered plane stress state, 1+ν

r=



X12 + X22 , and

θ = arccos(X1 /r ).

The rotation angle, which is also required for formulation of the boundary conditions, is derived from the displacements field by a formula

uKθ (X1 , X2 ) =

1 2



∂ uK1 ∂ uK2 − ∂ X2 ∂ X1



Similarly, for Mode II

σ (X1 , X2 ) = K 11

K σ22 (X1 , X2 ) = K σ12 (X1 , X2 ) =

K 2 π r 1/2 K 2 π r 1/2 K 2 π r 1/2

θ



3θ − sin 2 + cos cos 2 2 2 sin

θ

cos

θ

2 2

cos



θ

θ 2

cos

1 − sin

3θ , 2

θ 2

sin

3θ 2

 ,

 ,

  θ 3θ )(1 + ν ) (2κ + 3 ) sin + sin , 2 2    r 1 / 2 K θ 3θ uK2 (X1 , X2 ) = − (1 + ν ) (2κ − 3 ) cos + cos . 2E 2π 2 2 uK1 (X1 , X2 ) =



K r 2E 2π

1 / 2

References Ajdari, A., Jahromi, B.H., Papadopoulos, J., Nayeb-Hashemi, H., Vaziri, A., 2012. Hierarchical honeycombs with tailorable properties. Int. J. Solids Struct. 49 (11–12), 1413–1419. Arabnejad, S., Burnett Johnston, R., Pura, J.A., Singh, B., Tanzer, M., Pasini, D., 2016. High-strength porous biomaterials for bone replacement: a strategy to assess the interplay between cell morphology, mechanical properties, bone ingrowth and manufacturing constraints. Acta Biomater. 30, 345–356. Carpinteri, A., Pugno, N.M., 2008. Mechanics of hierarchical materials. Int. J. Fract. 150, 221–226. Chen, Q., Pugno, N., 2013. Bio-mimetic mechanisms of natural hierarchical materials: a review. J. Mech. Behav. Biomed. Mater. 19, 3–33. Cherkaev, A., Ryvkin, M., 2017. Damage propagation in 2d beam lattices: 2. design of an isotropic fault-tolerant lattice, accepted for archive of applied mechanics. Currey, J.D., 1984. The mechanical adaptations of bones. Princeton University Press, Princeton. Fan, H.L., Jin, F.N., Fang, D.N., 2008. Mechanical properties of hierarchical cellular materials. part i: analysis. Comp. Sci. Technol. 68, 3380–3387. Fleck, N.A., Qiu, X., 2007. The damage tolerance of elastic-brittle, two dimensional isotropic lattices. J. Mech. Phys. Solids 55, 562–588. Gao, H., 2006. Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials. Int. J. Fract. 138, 101–137. Gibson, L., 2012. The hierarchical structure and mechanics of plant materials. J. Royal Soc. Interface 9, 2749–2766. Gibson, L.J., Ashby, M.F., 1997. Cellular Solids. Structure and Properties. Cambridge University Press. Haghpanah, B., Oftadeh, R., Papadopoulos, J., Vaziri, A., 2013. Self-similar hierarchical honeycombs. Proceed. Royal Soc. A 469, 20130022. Jang, W., Kyriakides, S., 2015. On the buckling and crushing of expanded honeycomb. Int. J. Solids Struct. 91, 81–90. Karpov, E.G., Stephen, N.G., Dorofeev, D.L., 2002. On static analysis of finite repetitive structures by discrete fourier transform. Int. J. Solids Struct. 39, 4291–4310. Kooistra, G.W., Deshpande, V., Wadley, H.N.G., 2007. Hierarchical corrugated core sandwich panel concepts. J. Appl. Mech. 74, 259–268.

160

M. Ryvkin, R. Shraga / International Journal of Solids and Structures 152–153 (2018) 151–160

Kucherov, L., Ryvkin, M., 2014. Fracture toughness of open-cell kelvin foam. Int. J. Solids Struct. 51, 440–448. Lakes, R., 1993. Materials with structural hierarchy. Nature 361, 511–515. Lawn, B., 1993. Fracture of brittle solids. Cambridge University Press, Cambridge. Li, D., Yin, J., Dong, L., Lakes, R.S., 2017. Numerical analysis on mechanical behaviors of hierarchical cellular structures with negative poisson ratio. Smart Mater. Struct. 26, 025014. Lipperman, F., Ryvkin, M., Fuchs, M.B., 2007. Fracture toughness of two-dimensional cellular material with periodic microstructure. Int. J. Fract. 146, 279–290. Maiti, S.K., Ashby, M.F., Gibson, L.J., 1984. Fracture toughness of brittle cellular solids. Scr. Metall. 18, 213–217. Meza, L.R., Zelhofer, A.J., Clarke, N., Mateos, A.J., Kochmann, D.M., Greer, J.R., 2015. Resilient 3d hierarchical architected metamaterials. PNAS 112 (3), 11502–11507. Milton, G.W., 2002. Theory of composites. Cambridge University Press, Cambridge. Mousanezhad, D., Ebrahimi, H., Haghpanah, B., Ghosh, R., Ajdari, A., Hamouda, A.M.S., A. Vaziri, A., 2015a. Spiderweb honeycombs. Int. J. Solids Struct. 66, 218–227. Mousanezhad, D., S., B., Ghosh, R., Mahdi, E., Bertoldi, K., Vaziri, A., 2015b. Honeycomb phononic crystals with self-similar hierarchy. Phys. Rev. B 92, 104304. Oftadeh, R., Haghpanah, B., Vella, D., Boudaoud, A., Vaziri, A., 2014. Optimal fractal– like hierarchical honeycombs. Phys. Rev. Lett. 113, 104301. Phani, A.S., Fleck, N.A., 2008. Elastic boundary layers in two-dimensional isotropic lattices. J. Appl. Mech. Trans. ASME 75 (2). Ryvkin, M., 2008. Employing the discrete fourier transform in the analysis of multiscale problems. Intern. J. Multiscale Comput.Eng. 6, 435–449.

Ryvkin, M., Hadar, O., 2015. Employing of the discrete fourier transform for evaluation of crack-tip field in periodic materials. Int. J. Eng. Sci. 86, 10–19. Ryvkin, M., Hadar, O., Kucherov, L., 2017. Multiscale analysis of non-periodic stress state in composites with periodic microstructure. Int. J. Eng. Sci. 121, 167–181. Ryvkin, M., Nuller, B., 1997. Solution of quasieriodic fracture problems by representative cell method. Comput. Mech. 20, 145–149. Schmidt, I., Fleck, N.A., 2001. Ductile fracture of two dimensional cellular structures. Int. J. Fract. 111, 327–342. Sen, D., Buehler, M.J., 2011. Structural hierarchies define toughness and defecttolerance despite simple and mechanically inferior brittle building blocks. Scientific Reports (Nature) 35. doi:10.1038/srep0 0 035. Taylor, C.M., Smith, C.W., Miller, W., Evans, K.E., 2011. The effects of hierarchy on the in-plane elastic properties of honeycombs. Int. J. Solids Struct. 48, 1330–1339. Vigliotti, A., Pasini, D., 2012. Linear multiscale analysis and finite element validation of stretching and bending dominated lattice materials. Mech. Mater. 46, 57–68. Vigliotti, A., Pasini, D., 2013. Mechanical properties of hierarchical lattices. Mech. Mater. 62, 32–43. Zhu, H.X., Wang, Z.B., 2013. Size-dependent and tunable elastic and geometric properties of hierarchical nano-porous materials. Sci. Adv. Mater. 5, 1–10. Zhu, H.X., Yan, L., Zhang, R., Qiu, X.M., 2012. Size-dependent and tunable elastic properties of hierarchical honeycombs with square and equilateral triangular cells. Acta Mater. 60, 4927–4939. Zhu, H.X., Zhang, H.C., You, J.F., Kennedy, D., Wang, Z.B., Fan, T.X., Zhang, D., 2014. The elastic and geometrical properties of micro- and nano-structured hierarchical random irregular honeycombs. J. Mater. Sci. 49 (16), 5690–5702.