- Email: [email protected]
Ordered deformation localization in cellular mechanical metamaterials
Accepted Manuscript
Ordered deformation localization in cellular mechanical metamaterials Yafei Zhang , Yujue Wang , C.Q. Chen PII: DOI: Reference:
S0022-5096(18)30618-5 https://doi.org/10.1016/j.jmps.2018.08.025 MPS 3430
To appear in:
Journal of the Mechanics and Physics of Solids
Received date: Revised date: Accepted date:
16 July 2018 25 August 2018 27 August 2018
Please cite this article as: Yafei Zhang , Yujue Wang , C.Q. Chen , Ordered deformation localization in cellular mechanical metamaterials, Journal of the Mechanics and Physics of Solids (2018), doi: https://doi.org/10.1016/j.jmps.2018.08.025
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Ordered deformation localization in cellular mechanical metamaterials Yafei Zhang, Yujue Wang, C. Q. Chen 1 Department of Engineering Mechanics, CNMM and AML, Tsinghua University,
CR IP T
Beijing, 100084, P.R. China ABSTRACT Localized deformations in materials and structures are usually due to instability and are sensitive to various stochastic imperfections introduced during either the manufacturing processes or in-service life. As a result, they always occur in an uncertain manner. In this study, a two dimensional bistable cellular mechanical
AN US
metamaterial (CMM) is designed, showing highly ordered localized deformation under static loading. The metamaterial consists of periodically distributed elliptical holes of dual sizes. By varying the geometrical parameters of the holes, the CMM can be monostable or bistable. The responses of the CMM to uniaxial compression are systematically explored by finite element (FE) simulations. It is shown that the
M
monostable CMM can have either positive or negative Poisson’s ratio and undergo homogeneous deformation up to moderate uniaxial compression. However, the
ED
bistable CMM has negligible Poisson’s ratio and, most interestingly, shows highly ordered localized deformation zones, similar to the multi-soliton feature in fluids and
PT
physics. Based upon the unit cell of the CMM, a simplified spring chain model is developed. In the continuum limit, the model reduces to the Klein-Gordon equation.
CE
Analytical trivial solution and kink-antikink solution of the spring chain model are obtained. They are found to faithfully capture the FE simulated homogeneous deformation of the monostable CMM and the ordered localized deformations of the
AC
bistable CMM. The findings can provide insight into the ordered localized deformation, phase transformation, and static topological soliton of mechanical metamaterials. Keywords:
1
Cellular
mechanical
metamaterial;
Corresponding author, [email protected].
1
bistable;
ordered
localized
ACCEPTED MANUSCRIPT deformation; static topological soliton
1. Introduction Cellular materials and structures, e.g., cellular metals and lattice materials (Ashby et al., 2000; Deshpande et al., 2001), can have either stochastic or periodic microstructures. They have been shown to have attractive mechanical, thermal,
CR IP T
acoustic, and/or physical properties and are promising for various applications in load bearing structural components, energy absorption in impact protection packages, compact heat exchangers, sound and vibration attenuation, among others (Gibson and Ashby, 1997; Banhart, 2001; Lu et al., 2005; Gibson et al., 2010).
Inspired by photonic and acoustic metamaterials with tunable unusual bulk properties,
AN US
such as negative refraction, negative effective mass density and modulus (Liu et al., 2000; Smith et al., 2004; Liu et al., 2011; Kadic et al., 2013; Maznev et al., 2013), periodic cellular mechanical metamaterials (CMMs) have attracted increasing interest from scientists and engineers to achieve desirable mechanical properties. One of the key advantages of CMMs is that their properties depend mainly on their geometrical
M
architecture, rather than their constituent solids (Kochmann and Bertoldi, 2017). For instance, lattice mechanical metamaterials with superior stiffness, strength and
ED
ductility were produced at nanoscales for metals and ceramics (Schaedler et al., 2011; Meza et al., 2014). CMMs have also been designed to have various appealing
PT
properties, e.g., tunable auxetic behaviors (Kolken and Zadpoor, 2017), pattern transitions (Mullin et al., 2007; Overvelde and Bertoldi, 2014; Rafsanjani et al., 2015;
CE
Yang et al., 2016), controllable wave propagation (Bertoldi and Boyce, 2008; Chen et al., 2017), and topological nontrivial properties (Kane and Lubensky, 2014; Chen et
AC
al., 2018; Li et al., 2018). Of particular interest is the instability induced pattern transition in CMMs, because it provides added opportunities to fabricate complex microstructures with desired bulk properties (An et al., 2018; Jiang et al., 2018). A typical example is that, under uniform loading, a two dimensional cellular metamaterial (i.e., an elastic plate with periodic circular holes) firstly undergoes affine deformations. With the load approaching a critical value, their cell walls buckle and the circular holes suddenly transform into a new pattern with alternating orthogonal ellipses. Such pattern 2
ACCEPTED MANUSCRIPT transformations can be triggered by either shortwave length mode or longwave length mode, corresponding to homogeneous and localized deformations, respectively, and can be tuned by porosity and hole shape (Mullin et al., 2007; Overvelde and Bertoldi, 2014). Dynamic localization in the form of solitary waves (single soliton) has also been recently investigated in cellular metamaterials (Deng et al., 2017). In addition to mechanical stimuli, pH/temperature dual-responsive cellular membranes have been reported. Multiple nucleation sites induced by swelling load are formed and distinct
CR IP T
chiral buckling patterns separated by highly localized deformation boundaries are developed (Kang et al., 2013; Wu et al., 2014; Cao et al., 2016). The aforementioned configuration switching from one stable state to another corresponds to two different phases analogous to domain structures in phase transformation materials like ferromagnetic and ferroelectric ceramics (Fatuzzo and Merz, 1967; Remoissenet,
AN US
2010).
As can be seen, localized deformation is a salient feature accompanying the CMMs. In fact, deformation localization is ubiquitous in mechanical systems, e.g., necking (Audoly and Hutchinson, 2016) and shear band (Tvergaard, 1981) in solids, crushing
M
bands in foams and honeycomb structures (Papka and Kyriakides, 1999), wrinkle-tofold transitions of confined thin membranes (Pocivavsek et al., 2008; Diamant and
ED
Witten, 2011). It is widely recognized that deformation localization is sensitive to stochastic imperfections (Bažant and Cedolin, 2010). It is thus not surprising that the
PT
aforementioned localized deformations occur in a somewhat uncertain manner. However, we will show numerically and theoretically in the following that highly
CE
ordered localized deformation can be achieved for a carefully designed CMM under static uniaxial compression.
AC
This paper is organized as follows. First, the CMM and its corresponding finite element (FE) modeling are introduced in Section 2. The FE simulated geometrical parameters dependent Poisson’s ratio and deformation pattern under compression are shown in Section 3. A simplified rod-spring model and its resulted theoretical spring chain model for the simulated deformations are developed in Section 4. Comparisons between the theory predicted and FE simulated homogeneous and ordered localized deformations of the CMM are given in Section 5. The paper is concluded with a few remarks.
3
ACCEPTED MANUSCRIPT
2. Finite element modeling A schematic of the two dimensional (2D) CMM considered in this study is shown in Fig. 1, consisting of periodically distributed orthogonally oriented elliptical holes of dual sizes. A unit cell is marked with ai and bi (i=1 and 2) being the semi-axes in the x and y directions of the corresponding elliptical holes. For big holes (denoted by subscript ―1‖), the semi-major axis is in the y direction while that of the small holes
CR IP T
(denoted by subscript ―2‖) is in the x direction. All elliptical holes have the same ellipse ratio (i.e., b1 a1 a2 b2 ) but two different sizes. A dimensionless parameter a2 a1 is introduced to differentiate the sizes of the holes. All necks between adjacent holes have the same non-dimensional width t t a1 . The
AN US
dimensions of the unit cell are Lx 2a1 (1 t ) and Ly 2a1 ( t ) . The CMM shown in Fig. 1 comprises N x N y 20 5 unit cells, with N x and N y referring to the numbers of unit cells in the x and y directions, respectively. The sizes of the CMM are then N x Lx and N y Ly . It will be shown in the following that,
M
among the three non-dimensional parameters ( , , and t ), is of particular importance in determining the unusually ordered deformation localization in the
ED
CMMs under compression.
t
y
t
o
b1 a1
x
Ly
b2 a2 Lx
AC
CE
PT
,
Fig. 1 A cellular mechanical metamaterial with periodically distributed elliptical holes and its unit cell. The responses of the CMM to static uniaxial compression imposed in the y direction are simulated using the commercial FE software ABAQUS (version 6.1). Meshing using 4-node plane strain elements with reduced integration is adopted. Mesh sensitivity shows that about 10,000 elements for each unit cell are adequate to ensure
4
ACCEPTED MANUSCRIPT convergent numerical results. In the FE simulations, the CMM is sandwiched between two rigid plates and frictionless contact between the CMM and the plates is modelled by the analytical rigid surface/nodes interaction technique. Because we are only interested in moderate compression prior to full densification, self-contact of the holes is not considered in the modeling. Static loading is applied to the rigid plates by prescribed displacement .
An
CR IP T
artificial inter-dissipation factor of 10-4 is used to stabilize the simulations. This particular choice of dissipation is obtained by the trial and error method: Much bigger dissipation would give non-realistic deformed configurations while much less dissipation would not be able to stabilize simulation. The reaction force F and the relative displacement u y between the two rigid plates are used to calculate the
AN US
macroscopic compressive nominal stress y and engineering strain y via
y
uy
(1)
N y Ly
M
y
F w( N x Lx )
ED
where w is the out-of-plane width of CMMs. The solid material constituting the CMM is an incompressible elastomer and is
W C01 ( I1 3) C10 ( I 2 3)
(2)
CE
PT
assumed to following the Mooney–Rivlin hyperelastic model (Rivlin, 1948), i.e.,
where W is the strain energy density, C01 and C10 are material constants related to
AC
shear modulus by 2(C01 C10 ) , and I1 and I 2 are the first and second invariants of B (det B)1/3B . Here, the left Cauchy-Green tensor B is related to the deformation gradient tensor F by B F FT , with F x X and X and x referring to the initial and current configurations, respectively. The Cauchy stress can then be expressed as 2 σ (C10 I1 C01I 2 )I 2C10B 2C01B1 3
5
(3)
ACCEPTED MANUSCRIPT In the simulations, C10 0.284MPa and C01 0.432MPa are adopted. These values are typical for soft polyurethane Hei-Cast 8000.
3. Numerical results 3.1 Poisson’s ratio First, the effects of various are considered, with and t being fixed. A close
CR IP T
inspection of the CMM in Fig. 1 reveals that when is greater than 1 the unit cell displays a re-entrant feature and is thus expected to give rise to an auxetic behavior, i.e., negative Poisson’s ratio. Bertoldi et al. (2010) showed that the buckled phase of a metamaterial with circular holes is very similar to the CMM in Fig. 2(a) with 1 and is indeed auxetic. When 1 , on the contrary, the Poisson’s ratio of the CMM is
AN US
positive. Between them (i.e., 1 ), the Poisson’s ratio becomes negligible. The sign dependence of the Poisson’s ratio is illustrated in the phase diagram in Fig. 2(a). Figure 2(b) shows the FE simulated progressive Poisson’s ratio (i.e., v x y ) as a function of the applied compressive strain y for three selected values of 0.5 , 1,
M
and 1.5, corresponding to the two regimes and the boundary between them in Fig. 2(a). It is clear that the simulated Poisson’s ratios are positive, close to zero, and
ED
negative, consistent with the phase diagram in Fig. 2(a).
(b)
0
1 ,
CE
1
=0.5 =1.0 =1.5
0
-1
1, 0
AC
2
-x/y
PT
a2
(a) 1, 0
0.00
a1
0.05
y
0.10
Fig. 2 (a) Phase diagram for the Poisson’s ratios of CMMs and (b) FE simulated strain trajectories for CMMs under uniaxial compression with a1 3 mm, 2 , t 0.5 and selected values of 0.5, 1.0 and 1.5 .
6
ACCEPTED MANUSCRIPT
y / Ei
Abcdefghijklm 1234567890 0.03 0.09 =0.5 =1.0 =1.5 0.06
0.03
0.03
0.06
y
0.09
CR IP T
0.00 0.00
Fig. 3 FE simulated uniaxial compressive stress-strain curves of CMMs with
2.0, t 0.5 and three different values of 0.5 , 1.0 and 1.5 .
AN US
3.2 Deformation pattern
The associated macroscopic stress-strain curves for 0.5 , 1.0, and 1.5 are presented in Fig. 3, where the stress is normalized by the corresponding effective Young’s modulus Ei (subscript i referring to different models with 0.5 , 1.0, and 1.5) in the vertical direction. By doing so, the initial elastic responses of the three
M
curves collapse into a single master curve. One can see from Fig. 3 that, within the shown strain range, the stress varies monotonically with strain for 1 while
ED
softening is evident for 1. To further explore the deformation mechanisms responsible for the difference in the curves in Fig. 3, we show the corresponding
PT
unreformed and deformed ( y 10% ) configurations in Fig. 4. It can be seen that the deformed patterns for both 0.5 and 1.5 in Fig. 4(a-b) are homogenous. For 1
CE
(Fig. 4(c)), however, three buckling induced localized deformation zones are present. It has been generally accepted that deformation localization in materials and
AC
structures is usually sensitive to various stochastic imperfections either inherently or externally induced. As a result, deformation localization usually occurs in a somewhat random manner and prediction of its precise location is difficult. However, Fig. 4(c) shows a striking feature that the localized zones are highly ordered instead of randomly distributed. We have varied the ratio , all showing uniform deformation for 1 up to moderate compression (e.g., 20% in strain). For 1, highly ordered localized deformation is achieved if other geometrical parameters are properly chosen. Moreover, FE results show that the feature of the transition from uniform to
7
ACCEPTED MANUSCRIPT localized deformation in Fig. 4 persists, providing that N y is no more than a fraction of N x . Results are not shown for brevity. The condition of N y
N x implies that the
CMMs be long strips which are of particular interest in this study.
Undeformed
Deformed
0
(a)
CR IP T
0.5 0
(b)
1.5
0
(c)
AN US
1.0
Fig. 4 FE results of the undeformed and deformed configurations of CMMs with (a) deformed configurations is 10%.
ED
4. Theoretical analysis
M
0.5 , (b) 1.5 , and (c) 1.0 . The macroscopic compressive strain for the
In this section, a simplified theoretical model is developed to uncover the underlying
PT
mechanisms responsible for the simulated uniform and ordered localized deformations shown in Fig. 4. Day et al. (1992) studied the elastic moduli of thin
CE
plates with circular holes and found that, when the size of the necks between adjacent holes is small, deformation and elastic strain energy are concentrated in the neck
AC
region and local deformation at the necks dominates the entire plate. In the deformed configurations of the CMMs, when the neck width t compared to other geometrical parameters is small, localized neck deformation prevails. Accordingly, a simplified unit cell model with rods connected by torsional springs can be derived from the dotted lines superimposed on the original unit cell shown in Fig. 2(a). A rod-spring model subject to vertical compression of displacement u y 2 is given in Fig. 5. The vertical rods shown in blue have an initial inclined anticlockwise angle 0 and initial length L0 . Note that the sign of 0 corresponds to whether is less than, equal to or 8
ACCEPTED MANUSCRIPT greater than 1. The axial stiffness of the elastomer rods can be approximated by a linear spring with stiffness K , the bending stiffness of the yellow and blue necks are denoted by the stiffness C1 and C2 of the torsional springs, while the shearing stiffness of blue necks is represented by Cs .
u y 2
= rod
tension
=
C1
neck
bending stiffness
L0
=
C2 bending stiffness
AN US neck
stiffness
CR IP T
0
K
shear
Cs stiffness
Fig. 5 Schematic of the rod-spring model and the corresponding stiffness for the unit
M
cell of the CMM. The compressive displacement is denoted by u y 2 .
ED
With the imposed vertical compression u y 2 and ignoring asymmetric modes, the strain energy of the unit cell rod-spring model can be expressed as (4)
PT
U cell , L 4C1 ( 0 )2 2C2 ( 0 )2 2K ( L L0 )2
CE
where and L are the current inclined angle and length of the elastic rods, respectively, and are related to each other by L (L0 cos 0 )sec . Defining the
AC
rotation angle as
0 and taking approximation of sec( 0 ) in terms of
Taylor’s series for small , Eq. (4) can be rewritten up to O( 5 ) as U cell g0 g1 g2 2 g3 3 g4 4
with
9
(5)
ACCEPTED MANUSCRIPT g 0 2 f 0 K H 0 f 0 H 0 2 L0 g1 4 f1 K H 0 f 0 H 0 L0
g 2 2 2C1 C2 2 H 0 K f12 2 f 0 f 2 H 0 2 f 2 L0
(6)
g3 2 K 2 f3 L0 2 f1 f 2 f 0 f 3 H 0 H 0
g 4 2 K 2 f 4 L0 f 22 2 f1 f 3 2 f 0 f 4 H 0 H 0
CR IP T
where H 0 L0 cos 0 and fi (i 0,.., 4) are the coefficients of Taylor’s series expansion of sec( 0 ) (e.g., f0 sec0 , f1 sec0 tan 0 , etc.).
Recall that the deformations in Fig. 4 are independent of N y and remains so when
N y is no more than a fraction of N x . A quasi one dimensional (1D) rod-spring model
AN US
with N x unit cells in the x direction can thus be constructed for the CMM (see, Fig. 6(a)). Under uniform vertical compression u y 2 , the characteristic rotation angle of the n-th unit cell can be denoted by n . Based on the assumption of symmetric deformations within each unit cell, the strain energy of the 1D system can be obtained
M
as
ED
U system Cs d 02 (sin n 1 sin n ) 2 2C2 n2 U cell n n
= Cs d 02 (sin n 1 sin n )2 U eff n ,
(7)
PT
n
where Cs is the shearing stiffness of the torsional spring between adjacent unit cells
CE
and U eff n g0 g1 n (2C2 g2 ) n2 g3 n3 g4 n4
is
the
effective
on-site
potential of the n-th unit cell. Equation (7) indicates that the actual pattern of the
AC
deformed metamaterials is the result of the competition between the interaction energy Cs d02 (sin n1 sin n )2 and the on-site potential U eff . A spring chain model corresponding to Eq. (7) is shown in Fig. 6(b). Applying the principle of minimum potential energy to Eq. (7), i.e., U system n 0 , one gets the following equilibrium equation for the n-th unit cell of the 1D system in Fig. 6(b), i.e.,
10
ACCEPTED MANUSCRIPT
2Cs d02 (2 n n1 n1 +h2 )
U eff n n
0
(8)
with h2 sin 20 n21 2 n1 n 6 n2 2 n1n 12n / 4 . 4.1 Macroscopic homogenous deformation The FE results in Section 3 show that, under the vertical compression, the CMMs with
CR IP T
1 (i.e., the initial inclined angle 0 0 ) are bending-dominated and undergo macroscopically homogenous deformation. With the developed theory (8), we will uncover the deformation mechanisms associated with 0 0 and 0 0 . Without loss of generality, we first consider 0 0 . In this case, during the initial loading, the on-
satisfying U cell
AN US
site potential of unit cell U cell is monostable with a single well at cell 0 , indicating that there is only one stable state for the cell
unit cells. Therefore, all unit cells in the metamaterial operate at the bottom of the stable wells during deformation. Although it is hard to get the general solution of Eq.
M
(8), its trivial solution, i.e., constant rotation with n =eff (n 1,..., N ) can be obtained through U eff 0 , yielding
ED
eff
PT
eff
g3 8 g 2 g 4 3g32 Q , 5/3 4 g4 2 g 4Q 3 27/3 g 4
(9)
CE
where Q is given by
AC
Q 216 g 2 s g3 g 4 54 g33 432 g1 g 42 13
4(24 g 2 s g 4 9 g32 )3 (216 g 2 s g3 g 4 54 g33 432 g1 g 42 ) 2
(10)
with g2 s (2C2 g2 ) . The trivial constant rotation solution (9) indicates a macroscopically homogenous deformation. 4.2 Ordered localized deformation When 0 0 , on the contrary, the metamaterial under vertical compression shows a
11
ACCEPTED MANUSCRIPT nontrivial feature of clusters of ordered localized deformation zones (Fig. 4(c)). Most importantly, this feature cannot be represented by the response of a unit cell under the same loading condition. It is in fact a type of low-energy excitation above the lowest energy state analyzed in Section 4.1. In view of the 1D spring chain model (8), we will show in the following that the excited ordered localization zones are analogous to the domain walls in phase transformational ferromagnetic and ferroelectric ceramics (Fatuzzo and Merz, 1967; Remoissenet, 2010) and can be regarded as static
With 0 0 , Eq.(8) can be simplified as
U eff n n
0
AN US
2Cs d02 (2 n n 1 n 1 )
CR IP T
topological solitons (Deng et al., 2017; Zhang et al., 2018).
with
(11)
2
(12)
M
2 U eff n =U 0 n2 2 K 2 U0 U0
where U 0 K (3L0 8 )( L0 ) 6 0 and 6 K ( L0 ) (2C1 C2 ) . It is
ED
noted that the last two terms on the right-hand side of Eq. (12) are independent of n and will be omitted in the following derivation. Clearly, the effective on-site potential
AC
and
CE
PT
U eff n is a double-well potential (i.e., bi-stable) under the condition of 0 , i.e.,
KL20 8(C1 C2 )
(13)
c m
(14)
with c L0 L20 8(C1 C2 ) / K
2 and
m
8L0 3 . If the conditions (13) and
(14) are not satisfied, U eff n is monostable, as illustrated in the phase diagram in Fig. 6(c) with c .
12
ACCEPTED MANUSCRIPT (c)
(a)
(b) unit cell
interaction
KL02
…
bi-stable
on-site potential
1
…
8 monostable
CR IP T
C1 C2
Fig. 6 (a) The simplified rod-spring model of the CMM, (b) its equivalent 1D spring chain model, and (c) the phase diagram of the on-site potential with single and double wells for c m .
Consider the case that the metamaterial in Fig. 6 satisfies Eq. (13). Under small
AN US
loading (i.e., c ), U eff n is monostable at n 0 and the metamaterial undergoes uniform deformation. Increasing the compression beyond c , U eff n transits from monostable to bistable and each constituent unit cells in the metamaterial have two stable states. Under such a circumstance, instability and localization may
M
occur and the deformed pattern is governed by the competition between the double
ED
well on-site potential and the interaction energy between adjacent unit cells. When N x in Fig. 6, a continuum field to describe the evolution of the
PT
characteristic angle n can be derived with the following quantities:
CE
n (x n ) 2 2 (n ) (n ) ..., x 2 x 2
AC
n 1 (t ) [(n 1) ] (n )
where
is the lattice parameter in the x direction and
(15)
2(1 t ) . Based on Eq.
(15), the continuum limit of the equilibrium equation (11) can be obtained as
2 V 0. x 2 with
13
(16)
ACCEPTED MANUSCRIPT 1 2U 0 2 V 4 Cs d02 2 U0
2
2 1 2 2 . 4
(17)
which is the well-known Klein-Gordon equation (Remoissenet, 2010). Without loss of generality, taking x 0 0 , multiplying Eq. (16) by d u d x and integrating with respect to x yield
2 2 2 ( ) 0 x 2
/ 2a2
(18)
0 x 2V ( ) x 0
.
/ 2 2 and
/ 2 2 . Then, x and in Eq.
2
/ 2b2
Define
two
AN US
where
CR IP T
2
quantities
(18) can be related to each other through the elliptic integral
d
0
x
(a 2 2 )(b2 2 )
/ 2 d x / 2x 0
(19)
ED
(Remoissenet, 2010), i.e.,
M
the solution of which can be expressed by the elliptic integral of the first kind
a sn b / 2 x; m
CE
PT
with the period
AC
where m
P
4 b /2
EllipticK m /
(20)
(21)
a 2 b2 is the modulus of the elliptic sine function sn x; m , EllipticK m
is the elliptic integral of modulus m , a 2m (1 m) , and b a
m . Once the
period P is determined, m can be evaluated for a given problem. Formula (20) is the kink-antikink solution to the 1D Klein-Gordon system (16) and represents a periodic topological solitary excitation (Zhang et al., 2018). It is noted that Eq. (20) depends upon the stiffness of the rods and necks K , C1 , C2 and Cs . In turn, these stiffness are determined by the geometrical parameters and the constituent 14
ACCEPTED MANUSCRIPT solid of the CMM in Fig. 1. As will be shown later, the non-trivial solution (20) actually describes the ordered localized deformations of the metamaterial in Fig. 4(c). 4.3 Spring stiffness of the simplified rod-spring model Day et al. (1992) developed simplified rod-spring models to calculate the effective tension stiffness K , shear and bending stiffness ( Cs , C1 and C2 ) of the wall necks in
CR IP T
the plates with circular holes. They focused on thin necks (i.e., small t ). Coulais (2016) extended the model for elliptical holes and straight beams were replaced by curved beams when is larger than 1. Both neck width t and ellipse ratio were discussed. In order to get analytical results, two limiting cases such as extremely small and large t were considered. It should be pointed out that the model considered
AN US
by Coulais (2016) is different from the CMM shown in Fig. 1.
When the neck width t or the ellipse ratio is finite, however, it is difficult to get the analytical expressions of K , C1 , C2 and Cs for the considered CMMs. Here, we conduct extensive FE simulations to quantify the dependence of the stiffness upon t
M
and with 1 . A semi-empirical expression of the dependence is assumed to have
ED
the following form of
y , t G B t
(22)
PT
where y represents either K , C1 , C2 or Cs , G is the parameter for
CE
nondimensionalization of y , and B and are the functions to be determined by curve fitting the corresponding FE results. Therefore, the dimensionless stiffness
y , t
G . With dimensional analysis, G can be
AC
can be defined as y , t
expressed as G = ( 1 t ) , and a12 for tension, shear and bending stiffness, respectively. Based upon the unit cell model in Fig. 1 and the selected material constants in Eq. (3), linear perturbation analysis is used to calculate the tension stiffness K and shear stiffness Cs by applying uniaxial tension and simple shear, respectively. The periodic boundary conditions suggested by Chen et al. (1999) are adopted to ensure deformation compatibility. Figure 7(a) shows the uniaxial tension model and the 15
ACCEPTED MANUSCRIPT definition of K F (2 w) with F and
being the associated load and
displacement. The FE calculated normalized tension stiffness K as a function of nondimensional neck width t is shown in Fig. 7(b) for various ellipse ratios . From the results in Fig. 7(b), one can extract B and in Eq. (22) for K (denoted by BK and K with subscript K referring to tension stiffness K ). They are
governed by two asymptotes for 1 and
1 , i.e., K ( ) 1 0.121 0.349
as 1 and K ( ) 0.00257 0.557 as
1 . Guided by the observation,
K can be assumed to have the form of
ED
3.3
1.2
0.1
(d)
BK()
1
0.50
3
4.0
1.9
PT
CE
K()
0.1
FE Eq.(23)
0.52
AC
(23)
1
t
FE Eq.(24)
0.6
0.54
2
1 M
2.6
1
K
M
K F 2 w
0.56
(b) 1
(c)
M
AN US
K ( ) tanh A 1 1 1 1 1 2 (a)
CR IP T
plotted in Fig. 7(c) and 7(d) as symbols. It can be seen from Fig. 7(c) that BK is
0.4 0.2 1
4
2
3
4
Fig. 7 Dependence of the tension stiffness K upon the geometrical parameters t and
with 1 . (a) FE model to calculate K , (b) FE calculated K as a function of t for various , (c) and (d) curve fittings of the semi-empirical parameters K and
16
ACCEPTED MANUSCRIPT BK for K . Fitting Eq. (23) against the FE results in Fig. 7(c) gives A 2.11 and M 32.2 . The fitted curve is given in Fig. 7(c) as solid line, agreeing well with the FE results. From Fig. 7(d), it is evident that curve fitting of BK for K gives BK 1.08 0.395 0.0425 2 .
CR IP T
(a)
(24)
(b)
4.0 3.3
Cs
0.1
2.6
0.01
0.1
(c)
FE Eq.(25)1
1.2
1
t
(d)
FE Eq.(25)2
0.54
M
Bs()
1.3 1.2 1.1 2
3
0.53 4
1
PT
1
ED
s()
1.9
AN US
Cs F 2 w
2
3
4
CE
Fig. 8 Dependence of stiffness Cs upon the geometrical parameters t and with
1 . (a) FE model for calculating Cs , (b) FE calculated Cs as a function of t for
AC
various , (c) and (d) curve fittings of the semi-empirical parameters s and
Bs for Cs .
The FE model for calculating the shear stiffness Cs for various geometrical parameters is illustrated in Fig. 8(a). The obtained FE results are given in Fig. 8(b), showing that the shear stiffness Cs can also be approximated by Eq. (22) with curve fitted parameters
17
ACCEPTED MANUSCRIPT s ( ) 1.45 0.0957 ,
(25)
Bs 0.552 0.0119 0.00163 .
where subscript s denotes shear stiffness. At 1 , s ( ) 1.35 , very close to the results of circular hole obtained by Day et al. (1992).
(a)
(b)
C1 M w
3.3
CR IP T
C1
0.1
2.6
0.01
1.9
0.001
0.1 (d)
FE Eq.(27)1
t
2.3
2
3
4
ED
1
M
2.2
1.2
FE Eq.(27)2
0.5
B1()
1()
2.4
2.1
1
AN US
(c) 2.5
4.0
0.4 0.3
1
2
3
4
Fig. 9 Dependence of stiffness C1 upon the geometrical parameters t and with
PT
1 . (a) FE model to calculate C1 ; (b) FE calculated C1 as a function of t for
CE
various ; (c) and (d) curve fittings of the semi-empirical parameters 1 and
AC
B1 for C1 .
To calculate the bending stiffness C1 and C2 , the FE models given in Figs. 9(a) and 10(a) are used. The models are analogous to that of Day et al. (1992). The FE calculated dependence of C1 and C2 upon t with selected values are shown in Figs. 9(b) and 10(b), respectively. Again, fitting the numerical results against formula (22) gives
18
ACCEPTED MANUSCRIPT C1 G B1 t 1
(26)
C2 G B2 t 2 with
1 2.63 0.154
(27)
and
2 2.58 0.102 B2 0.592 0.0729
CR IP T
B1 0.608 0.0747
(28)
AN US
for the considered range of t [0.05, 0.7] and [1.2, 4] . As expected, when 1 ,
1 and 2 approach the theoretical results given by Day et al. (1992) for circular holes, i.e., 2.5.
3.3
0.01
2.6
0.001
1.9 0.1
(d)
FE Eq.(28)1
CE
(c) 2.5
AC
1
t
2.1
1.2
FE Eq.(28)2
0.5
B2()
2.3
4.0
0.1
C2
ED
C2 M w
PT
2()
(b)
M
(a)
0.4 0.3
1
2
3
4
1
2
3
4
Fig. 10 Dependence of stiffness C2 upon the geometrical parameters t and with
1 . (a) FE model to calculate C2 ; (b) FE calculated C2 as a function of t for various ; (c) and (d) curve fittings of the semi-empirical parameter 2 and 19
ACCEPTED MANUSCRIPT B2 for C2 .
5. Comparison between FE and theoretical results We now compare the developed theoretical predictions in Section 4 with the FE simulated results. From the FE calculated deformed configurations in Fig. 4, the characteristic rotation for each unit cell located at the middle row can be extracted
CR IP T
using the definition indicated in the simplified model in Fig. 5. The corresponding FE results on the variation of along the x direction are shown in Fig. 11 as symbols for 0.5 , 1, and 1.5. As can be seen, the characteristic rotations for 0.5 and 1.5 are almost constant, a direct consequence of the homogeneous deformations in Fig. 4(a-b). When 1 , on the contrary, is not a constant and its sign alternatively
AN US
changes along x . Note that is the relative rod rotation. A negative with 1 indicates clockwise rotation of the originally vertical rods and a localized zone consisting of re-entrant microstructures. The extracted constant for 0.5 and 1.5 and regular variation of with 1 correlate reasonably well the FE deformation patterns in Fig. 4, showing that the defined characteristic rotation is a
M
suitable parameter to quantify the deformation of the CMM under compression.
ED
For the purpose of comparison, the corresponding theoretical predictions obtained in Section 4 are also included in Fig. 11 as lines. When 0.5 and 1.5, the effective
PT
on-site potential U eff in Eq. (8) is monostable and the trivial homogeneous deformation solution (9) applies. Specifically, Eq. (9) gives eff 0.159 for
CE
0.5 , with the stiffness K 0.934MPa , C1 0.621N , and C2 0.590N obtained from FE simulations. Similarly, Eq. (9) gives eff 0.183 for 1.5 , with
AC
K 1.47MPa , C1 1.17N , and C2 1.13N .
For the case of 1 , the effective on-site potential U eff in Eq. (8) becomes bistable when a moderate compression (e.g., 10%) is applied. Then, the kink-antikink solution (20) is adopted, with the stiffness K 1.34MPa , C1 1.19N , C2 1.10N ,
Cs 0.279MPa , and m 0.35 . As can be seen from Fig. 11, the theoretical predictions for 0.5 and 1.5 agree well with the FE results. For 1 , the overall agreement between the theoretical solution and the FE results is also very good: both 20
ACCEPTED MANUSCRIPT the amplitude and the spacing of the localized deformation zones are captured. Note that the kink-antikink solution is for periodic systems while the system in Fig. 4c is of finite size. It is thus not surprise that, for 1 , the theoretical kink-antikink solution (20) slightly deviates from the FE results at the two ends. =0.5
Theory
=1.0
=1.5
0.4 0.0 -0.4
=0.5
FE -5
-0.8 -10
=1.0 0 x/l
(a)
CR IP T
Rotation angle
0.8
=1.5 5
10
Fig. 11 Comparison of the FE and theoretical predictions of the characteristic rotation
AN US
angle for the deformed CMMs shown in Fig. 4(c). (a) (b)
M
0.50
Theory
0.25
ED
0.00 -0.25 -10
PT
Rotation angle
(c)
-5
CE
Fig. 12 A quasi 1D CMM
FE
0 x/l
5
10
with the geometry parameters a1 3.0 mm, 1 ,
t 0.55 and 2.33 . (a) The initial configuraiton, (b) the ordered localized
AC
deformation under compressive strain 9.7% , and (c) the theoretical kink-antikink solution and FE predicted variation of the characteristic rotation angle. To further validate the theoretical analysis in conjunction with the determined semiempirical results in Section 4, a quasi-one dimensional CMM with N x N y 20 1 and different geometrical parameters 1 , t 0.55 and 2.33 is considered. The particular choice of the geometrical parameters is not used in determining the semiempirical results in Eqs. (23-28). The undeformed and deformed ( 9.7% )
21
ACCEPTED MANUSCRIPT configurations are shown in Fig. 12(a-b), respectively. Regularly spaced localized deformation zones are present in the deformed shape. From the semi-empirical formulae (Eqs. (23-28)), one has K 1.30MPa , Cs 0.340MPa , C1 1.40N ,
C2 1.40N , and m 0.45 . Accordingly, Eq. (20) predicts an ordered variation of characteristic rotation angle and is plotted in Fig. 12(c) as the solid line, showing good
(b) B C
0.04
Rotation angle
A
0.00 0.00
0.04
y
D
0.08
Theory
FE -5
AN US
A
0.8 0.4 0.0 -0.4 -0.8 -10
D
B C
0.02
(c)
0
A
A
B
B
M
y / E
(a)
CR IP T
agreement with the FE result.
C
0 x/l
C
D
D
5
10
ED
Fig. 13 Deformation evolution of a CMM ( a1 3.0 mm, 1 , t 0.3 and 1.5 ) under uniaxial compression. (a) The FE predicted normalized compressive stress-
PT
strain curve, (b) the original and FE simulated deformed configurations at selected loading levels marked in the stress-strain curve, and (c) the theoretical and FE
CE
predicted characteristic rotation angle. It is noted that the periodic kink-antikink solution (20) depends not only upon the
AC
geometrical parameters but also the magnitude of compression (i.e., ). To further validate the theoretical analysis in Section 4, a quasi-one dimensional CMM with
N x N y 20 1 and geometrical parameters 1 , t 0.3 and 1.5 is
considered. Again, the particular choice of the geometrical parameters is not used in determining the semi-empirical results (Eqs. (23-28)). The FE simulated compressive stress-strain curve and the deformation configurations at selected load levels (A–D) are shown in Fig. 13(a-b), respectively. For small deformations (e.g., Point A in the linear region), the CMM undergoes uniform deformation. With the load increasing, 22
ACCEPTED MANUSCRIPT the stress-strain curve gradually deviates from its linear response. Upon reaching Point B (slightly before the peak load), localized deformation zones start to develop. When the loading increases beyond the critical strain at the peak, strain softening accompanied by localized deformation zones can be observed (Point C). Further increasing the strain from C to D, minor strain hardening can be observed. Interestingly, there is little change in the locations of the localized zones, showing that
CR IP T
deformation localized zones are robust. With Eq. (20), the characteristic rotation angles of the quasi-one dimensional CMM at different loading levels (i.e., A–D) can be theoretically predicted ( m 0.17, 0.41, 0.53 and 0.79 for A–D, respectively) and are plotted in Fig. 13(c) as lines, together with the FE results denoted by symbols. It is clear that the FE simulated evolution of the localization zones as characterized by the
AN US
amplitude and spacing of the rotation angle is faithfully captured by the developed theoretical model.
6. Concluding remarks
A two dimensional CMM comprising of elliptical holes of dual sizes is designed.
M
Systematical FE simulations are conducted to explore its responses to uniaxial compression. It is shown that the CMM can have positive, negative, or zero Poisson’s
ED
ratio and can also be monostable or bistable, depending upon its geometrical parameters. Under moderate compression, the bistable CMM show ordered and highly
PT
localized deformation zones in the direction perpendicular to the compression, in
CE
sharp contrast to the homogeneous deformation of the monostable CMM. A simplified rod-spring model is derived from the unit cell of the CMM to uncover the deformation mechanisms and a theoretical spring chain model is developed.
AC
Dependence of the model parameters upon the geometrical parameters of the CMM is quantified by FE simulations. In the continuum limit, the theoretical model reduces to the Klein-Gordon equation. Analytical trivial solution (i.e., constant rotation) and kink-antikink solution are obtained, corresponding to the models with monostable and bistable potentials, respectively. The constant rotation solution indicates homogeneous deformation and correlates well to the FE simulated deformed configurations of the monostable CMM. The kink-antikink solution faithfully captures the features (i.e., periodic spacing and magnitude) of the ordered localization zones of the bistable 23
ACCEPTED MANUSCRIPT CMM under compression. Note that the feature of kink-antikink pairs is ubiquitous (e.g., domain walls in phase transform materials and multi-solitary excitation in waves). The excellent agreement between the FE simulations and theoretical predictions implies that the ordered deformation localization zones are similar to multi-soliton in fluids and physics. This study can thus provide insight into the excitation of ordered localized deformation, phase transformation and static
CR IP T
topological soliton, as well as the design of mechanical metamaterials. Acknowledgements This work is supported by the National Natural Science Foundation of China (Nos. 11732007 and 11472149) and the Science Challenge Project of China (Grant No. TZ2018007).
AN US
References
An, S., Kim, B., Kwon, S., Moon, G., Lee, M., Jhe, W., 2018. Bifurcation-enhanced ultrahigh sensitivity of a buckled cantilever. Proc. Natl. Acad. Sci. 115, 2884– 2889. https://doi.org/10.1073/pnas.1716067115
M
Ashby, M., Evans, A., Fleck, N., Gibson, L., Hutchinson, J., Wadley, H., 2000. Metal Foams: A Design Guide. Butterworth-Heinemann, Oxford.
ED
Audoly, B., Hutchinson, J.W., 2016. Analysis of necking based on a one-dimensional model. J. Mech. Phys. Solids 97, 68–91. https://doi.org/10.1016/j.jmps.2015.12.018
PT
Banhart, J., 2001. Manufacture, characterisation and application of cellular metals and metal foams. Prog. Mater. Sci. 46, 559–632. https://doi.org/10.1016/S00796425(00)00002-5
CE
Bažant, Z.P., Cedolin, L., 2010. Stability of structures: elastic, inelastic, fracture and damage theories. World Scientific Pub, Hackensack, NJ.
AC
Bertoldi, K., Boyce, M., 2008. Mechanically triggered transformations of phononic band gaps in periodic elastomeric structures. Phys. Rev. B 77, 052105. https://doi.org/10.1103/PhysRevB.77.052105 Bertoldi, K., Reis, P.M., Willshaw, S., Mullin, T., 2010. Negative Poisson’s ratio behavior induced by an elastic instability. Adv. Mater. 22, 361–366. https://doi.org/10.1002/adma.200901956 Cao, B., Wu, G., Xia, Y., Yang, S., 2016. Buckling into single-handed chiral structures from pH-sensitive hydrogel membranes. Extreme Mech. Lett. 7, 49–54. https://doi.org/10.1016/j.eml.2015.12.011 Chen, C., Lu, T.J., Fleck, N.A., 1999. Effect of imperfections on the yielding of twodimensional foams. J. Mech. Phys. Solids 47, 2235–2272. 24
ACCEPTED MANUSCRIPT https://doi.org/10.1016/s0022-5096(99)00030-7 Chen, H., Nassar, H., Huang, G.L., 2018. A study of topological effects in 1D and 2D mechanical lattices. J. Mech. Phys. Solids. https://doi.org/10.1016/j.jmps.2018.04.013 Chen, Y., Li, T., Scarpa, F., Wang, L., 2017. Lattice metamaterials with mechanically tunable poisson’s ratio for vibration control. Phys. Rev. Appl. 7, 024012.
CR IP T
Coulais, C., 2016. Periodic cellular materials with nonlinear elastic homogenized stress-strain response at small strains. Int. J. Solids Struct. 97–98, 226–238. https://doi.org/10.1016/j.ijsolstr.2016.07.025 Day, A.R., Snyder, K.A., Garboczi, E.J., Thorpe, M.F., 1992. The elastic moduli of a sheet containing circular holes. J. Mech. Phys. Solids 40, 1031–1051.
AN US
Deng, B., Raney, J.R., Tournat, V., Bertoldi, K., 2017. Elastic vector solitons in soft architected materials. Phys. Rev. Lett. 118, 204102. https://doi.org/10.1103/PhysRevLett.118.204102 Deshpande, V., Ashby, M., Fleck, N., 2001. Foam topology: bending versus stretching dominated architectures. Acta Mater. 49, 1035–1040.
M
Diamant, H., Witten, T.A., 2011. Compression induced folding of a sheet: an integrable system. Phys. Rev. Lett. 107, 164302. https://doi.org/10.1103/PhysRevLett.107.164302 Fatuzzo, E., Merz, W.J., 1967. Ferroelectricity. North-Holland Pub. Co., Amsterdam.
ED
Gibson, L.J., Ashby, M.F., 1997. Cellular solids: Structure and properties. Cambridge University Press, Cambridge.
PT
Gibson, L.J., Ashby, M.F., Harley, B.A., 2010. Cellular materials in nature and medicine. Cambridge University Press, Cambridge.
CE
Jiang, Y., Liu, Z.Y., Matsuhisa, N., Qi, D., Leow, W.R., Yang, H., Yu, J., Chen, G., Liu, Y., Wan, C., Liu, Z.J., Chen, X., 2018. Auxetic mechanical metamaterials to enhance sensitivity of stretchable strain sensors. Adv. Mater. 30, 1706589. https://doi.org/10.1002/adma.201706589
AC
Kadic, M., Bückmann, T., Schittny, R., Wegener, M., 2013. Metamaterials beyond electromagnetism. Rep. Prog. Phys. 76, 126501. https://doi.org/10.1088/00344885/76/12/126501 Kane, C.L., Lubensky, T.C., 2014. Topological boundary modes in isostatic lattices. Nat. Phys. 10, 39–45. https://doi.org/10.1038/nphys2835 Kang, S.H., Shan, S., Noorduin, W.L., Khan, M., Aizenberg, J., Bertoldi, K., 2013. Buckling-induced reversible symmetry breaking and amplification of chirality using supported cellular structures. Adv. Mater. 25, 3380–3385. https://doi.org/10.1002/adma.201300617
25
ACCEPTED MANUSCRIPT Kochmann, D.M., Bertoldi, K., 2017. Exploiting Microstructural Instabilities in Solids and Structures: From Metamaterials to Structural Transitions. Appl. Mech. Rev. 69, 050801-050801–24. https://doi.org/10.1115/1.4037966 Kolken, H.M.A., Zadpoor, A.A., 2017. Auxetic mechanical metamaterials. RSC Adv. 7, 5111–5129. https://doi.org/10.1039/C6RA27333E Li, S., Zhao, D., Niu, H., Zhu, X., Zang, J., 2018. Observation of elastic topological states in soft materials. Nat. Commun. 9, 1370. https://doi.org/10.1038/s41467-018-03830-8
CR IP T
Liu, X.N., Hu, G.K., Huang, G.L., Sun, C.T., 2011. An elastic metamaterial with simultaneously negative mass density and bulk modulus. Appl. Phys. Lett. 98, 251907. https://doi.org/10.1063/1.3597651 Liu, Z., Zhang, X., Mao, Y., Zhu, Y.Y., Yang, Z., Chan, C.T., Sheng, P., 2000. Locally resonant sonic materials. Science 289, 1734–1736.
AN US
Lu, T.J., Valdevit, L., Evans, A.G., 2005. Active cooling by metallic sandwich structures with periodic cores. Prog. Mater. Sci. 50, 789–815. https://doi.org/10.1016/j.pmatsci.2005.03.001 Maznev, A.A., Every, A.G., Wright, O.B., 2013. Reciprocity in reflection and transmission: What is a `phonon diode’? Wave Motion 50, 776–784. https://doi.org/10.1016/j.wavemoti.2013.02.006
ED
M
Meza, L.R., Das, S., Greer, J.R., 2014. Strong, lightweight, and recoverable threedimensional ceramic nanolattices. Science 345, 1322–1326. https://doi.org/10.1126/science.1255908
PT
Mullin, T., Deschanel, S., Bertoldi, K., Boyce, M., 2007. Pattern transformation triggered by deformation. Phys. Rev. Lett. 99, 084301. https://doi.org/10.1103/PhysRevLett.99.084301
CE
Overvelde, J.T.B., Bertoldi, K., 2014. Relating pore shape to the non-linear response of periodic elastomeric structures. J. Mech. Phys. Solids 64, 351–366. https://doi.org/10.1016/j.jmps.2013.11.014
AC
Papka, S.D., Kyriakides, S., 1999. In-plane biaxial crushing of honeycombs—: Part II: Analysis. Int. J. Solids Struct. 36, 4397–4423. https://doi.org/10.1016/S0020-7683(98)00225-X Pocivavsek, L., Dellsy, R., Kern, A., Johnson, S., Lin, B., Lee, K.Y.C., Cerda, E., 2008. Stress and fold localization in thin elastic membranes. Science 320, 912–916. Rafsanjani, A., Akbarzadeh, A., Pasini, D., 2015. Snapping mechanical metamaterials under tension. Adv. Mater. 27, 5931–5935. https://doi.org/10.1002/adma.201502809 Remoissenet, M., 2010. Waves called solitons: concepts and experiments. Springer Science & Business Media., Berlin. 26
ACCEPTED MANUSCRIPT Rivlin, R.S., 1948. Large elastic deformations of isotropic materials IV. Further developments of the general theory. Philos. Trans. R. Soc. Math. Phys. Eng. Sci. 241, 379–397. https://doi.org/10.1098/rsta.1948.0024 Schaedler, T.A., Jacobsen, A.J., Torrents, A., Sorensen, A.E., Lian, J., Greer, J.R., Valdevit, L., Carter, W.B., 2011. Ultralight metallic microlattices. Science 334, 962–965. https://doi.org/10.1126/science.1211649
CR IP T
Smith, D.R., Pendry, J.B., Wiltshire, M.C.K., 2004. Metamaterials and negative refractive index. Science 305, 788–792. https://doi.org/10.1126/science.1096796 Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389–407. https://doi.org/10.1007/BF00036191 Wu, G., Xia, Y., Yang, S., 2014. Buckling, symmetry breaking, and cavitation in periodically micro-structured hydrogel membranes. Soft Matter 10, 1392– 1399. https://doi.org/10.1039/C3SM51640G
AN US
Yang, D., Jin, L., Martinez, R.V., Bertoldi, K., Whitesides, G.M., Suo, Z., 2016. Phase-transforming and switchable metamaterials. Extreme Mech. Lett. 6, 1– 9. https://doi.org/10.1016/j.eml.2015.11.004
AC
CE
PT
ED
M
Zhang, Y., Li, B., Zheng, Q.S., Chen, C.Q., 2018. Static periodic topological soliton in mechanical metamaterials. Under preparation.
27
Ordered deformation localization in cellular mechanical metamaterials Yafei Zhang , Yujue Wang , C.Q. Chen PII: DOI: Reference:
S0022-5096(18)30618-5 https://doi.org/10.1016/j.jmps.2018.08.025 MPS 3430
To appear in:
Journal of the Mechanics and Physics of Solids
Received date: Revised date: Accepted date:
16 July 2018 25 August 2018 27 August 2018
Please cite this article as: Yafei Zhang , Yujue Wang , C.Q. Chen , Ordered deformation localization in cellular mechanical metamaterials, Journal of the Mechanics and Physics of Solids (2018), doi: https://doi.org/10.1016/j.jmps.2018.08.025
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Ordered deformation localization in cellular mechanical metamaterials Yafei Zhang, Yujue Wang, C. Q. Chen 1 Department of Engineering Mechanics, CNMM and AML, Tsinghua University,
CR IP T
Beijing, 100084, P.R. China ABSTRACT Localized deformations in materials and structures are usually due to instability and are sensitive to various stochastic imperfections introduced during either the manufacturing processes or in-service life. As a result, they always occur in an uncertain manner. In this study, a two dimensional bistable cellular mechanical
AN US
metamaterial (CMM) is designed, showing highly ordered localized deformation under static loading. The metamaterial consists of periodically distributed elliptical holes of dual sizes. By varying the geometrical parameters of the holes, the CMM can be monostable or bistable. The responses of the CMM to uniaxial compression are systematically explored by finite element (FE) simulations. It is shown that the
M
monostable CMM can have either positive or negative Poisson’s ratio and undergo homogeneous deformation up to moderate uniaxial compression. However, the
ED
bistable CMM has negligible Poisson’s ratio and, most interestingly, shows highly ordered localized deformation zones, similar to the multi-soliton feature in fluids and
PT
physics. Based upon the unit cell of the CMM, a simplified spring chain model is developed. In the continuum limit, the model reduces to the Klein-Gordon equation.
CE
Analytical trivial solution and kink-antikink solution of the spring chain model are obtained. They are found to faithfully capture the FE simulated homogeneous deformation of the monostable CMM and the ordered localized deformations of the
AC
bistable CMM. The findings can provide insight into the ordered localized deformation, phase transformation, and static topological soliton of mechanical metamaterials. Keywords:
1
Cellular
mechanical
metamaterial;
Corresponding author, [email protected].
1
bistable;
ordered
localized
ACCEPTED MANUSCRIPT deformation; static topological soliton
1. Introduction Cellular materials and structures, e.g., cellular metals and lattice materials (Ashby et al., 2000; Deshpande et al., 2001), can have either stochastic or periodic microstructures. They have been shown to have attractive mechanical, thermal,
CR IP T
acoustic, and/or physical properties and are promising for various applications in load bearing structural components, energy absorption in impact protection packages, compact heat exchangers, sound and vibration attenuation, among others (Gibson and Ashby, 1997; Banhart, 2001; Lu et al., 2005; Gibson et al., 2010).
Inspired by photonic and acoustic metamaterials with tunable unusual bulk properties,
AN US
such as negative refraction, negative effective mass density and modulus (Liu et al., 2000; Smith et al., 2004; Liu et al., 2011; Kadic et al., 2013; Maznev et al., 2013), periodic cellular mechanical metamaterials (CMMs) have attracted increasing interest from scientists and engineers to achieve desirable mechanical properties. One of the key advantages of CMMs is that their properties depend mainly on their geometrical
M
architecture, rather than their constituent solids (Kochmann and Bertoldi, 2017). For instance, lattice mechanical metamaterials with superior stiffness, strength and
ED
ductility were produced at nanoscales for metals and ceramics (Schaedler et al., 2011; Meza et al., 2014). CMMs have also been designed to have various appealing
PT
properties, e.g., tunable auxetic behaviors (Kolken and Zadpoor, 2017), pattern transitions (Mullin et al., 2007; Overvelde and Bertoldi, 2014; Rafsanjani et al., 2015;
CE
Yang et al., 2016), controllable wave propagation (Bertoldi and Boyce, 2008; Chen et al., 2017), and topological nontrivial properties (Kane and Lubensky, 2014; Chen et
AC
al., 2018; Li et al., 2018). Of particular interest is the instability induced pattern transition in CMMs, because it provides added opportunities to fabricate complex microstructures with desired bulk properties (An et al., 2018; Jiang et al., 2018). A typical example is that, under uniform loading, a two dimensional cellular metamaterial (i.e., an elastic plate with periodic circular holes) firstly undergoes affine deformations. With the load approaching a critical value, their cell walls buckle and the circular holes suddenly transform into a new pattern with alternating orthogonal ellipses. Such pattern 2
ACCEPTED MANUSCRIPT transformations can be triggered by either shortwave length mode or longwave length mode, corresponding to homogeneous and localized deformations, respectively, and can be tuned by porosity and hole shape (Mullin et al., 2007; Overvelde and Bertoldi, 2014). Dynamic localization in the form of solitary waves (single soliton) has also been recently investigated in cellular metamaterials (Deng et al., 2017). In addition to mechanical stimuli, pH/temperature dual-responsive cellular membranes have been reported. Multiple nucleation sites induced by swelling load are formed and distinct
CR IP T
chiral buckling patterns separated by highly localized deformation boundaries are developed (Kang et al., 2013; Wu et al., 2014; Cao et al., 2016). The aforementioned configuration switching from one stable state to another corresponds to two different phases analogous to domain structures in phase transformation materials like ferromagnetic and ferroelectric ceramics (Fatuzzo and Merz, 1967; Remoissenet,
AN US
2010).
As can be seen, localized deformation is a salient feature accompanying the CMMs. In fact, deformation localization is ubiquitous in mechanical systems, e.g., necking (Audoly and Hutchinson, 2016) and shear band (Tvergaard, 1981) in solids, crushing
M
bands in foams and honeycomb structures (Papka and Kyriakides, 1999), wrinkle-tofold transitions of confined thin membranes (Pocivavsek et al., 2008; Diamant and
ED
Witten, 2011). It is widely recognized that deformation localization is sensitive to stochastic imperfections (Bažant and Cedolin, 2010). It is thus not surprising that the
PT
aforementioned localized deformations occur in a somewhat uncertain manner. However, we will show numerically and theoretically in the following that highly
CE
ordered localized deformation can be achieved for a carefully designed CMM under static uniaxial compression.
AC
This paper is organized as follows. First, the CMM and its corresponding finite element (FE) modeling are introduced in Section 2. The FE simulated geometrical parameters dependent Poisson’s ratio and deformation pattern under compression are shown in Section 3. A simplified rod-spring model and its resulted theoretical spring chain model for the simulated deformations are developed in Section 4. Comparisons between the theory predicted and FE simulated homogeneous and ordered localized deformations of the CMM are given in Section 5. The paper is concluded with a few remarks.
3
ACCEPTED MANUSCRIPT
2. Finite element modeling A schematic of the two dimensional (2D) CMM considered in this study is shown in Fig. 1, consisting of periodically distributed orthogonally oriented elliptical holes of dual sizes. A unit cell is marked with ai and bi (i=1 and 2) being the semi-axes in the x and y directions of the corresponding elliptical holes. For big holes (denoted by subscript ―1‖), the semi-major axis is in the y direction while that of the small holes
CR IP T
(denoted by subscript ―2‖) is in the x direction. All elliptical holes have the same ellipse ratio (i.e., b1 a1 a2 b2 ) but two different sizes. A dimensionless parameter a2 a1 is introduced to differentiate the sizes of the holes. All necks between adjacent holes have the same non-dimensional width t t a1 . The
AN US
dimensions of the unit cell are Lx 2a1 (1 t ) and Ly 2a1 ( t ) . The CMM shown in Fig. 1 comprises N x N y 20 5 unit cells, with N x and N y referring to the numbers of unit cells in the x and y directions, respectively. The sizes of the CMM are then N x Lx and N y Ly . It will be shown in the following that,
M
among the three non-dimensional parameters ( , , and t ), is of particular importance in determining the unusually ordered deformation localization in the
ED
CMMs under compression.
t
y
t
o
b1 a1
x
Ly
b2 a2 Lx
AC
CE
PT
,
Fig. 1 A cellular mechanical metamaterial with periodically distributed elliptical holes and its unit cell. The responses of the CMM to static uniaxial compression imposed in the y direction are simulated using the commercial FE software ABAQUS (version 6.1). Meshing using 4-node plane strain elements with reduced integration is adopted. Mesh sensitivity shows that about 10,000 elements for each unit cell are adequate to ensure
4
ACCEPTED MANUSCRIPT convergent numerical results. In the FE simulations, the CMM is sandwiched between two rigid plates and frictionless contact between the CMM and the plates is modelled by the analytical rigid surface/nodes interaction technique. Because we are only interested in moderate compression prior to full densification, self-contact of the holes is not considered in the modeling. Static loading is applied to the rigid plates by prescribed displacement .
An
CR IP T
artificial inter-dissipation factor of 10-4 is used to stabilize the simulations. This particular choice of dissipation is obtained by the trial and error method: Much bigger dissipation would give non-realistic deformed configurations while much less dissipation would not be able to stabilize simulation. The reaction force F and the relative displacement u y between the two rigid plates are used to calculate the
AN US
macroscopic compressive nominal stress y and engineering strain y via
y
uy
(1)
N y Ly
M
y
F w( N x Lx )
ED
where w is the out-of-plane width of CMMs. The solid material constituting the CMM is an incompressible elastomer and is
W C01 ( I1 3) C10 ( I 2 3)
(2)
CE
PT
assumed to following the Mooney–Rivlin hyperelastic model (Rivlin, 1948), i.e.,
where W is the strain energy density, C01 and C10 are material constants related to
AC
shear modulus by 2(C01 C10 ) , and I1 and I 2 are the first and second invariants of B (det B)1/3B . Here, the left Cauchy-Green tensor B is related to the deformation gradient tensor F by B F FT , with F x X and X and x referring to the initial and current configurations, respectively. The Cauchy stress can then be expressed as 2 σ (C10 I1 C01I 2 )I 2C10B 2C01B1 3
5
(3)
ACCEPTED MANUSCRIPT In the simulations, C10 0.284MPa and C01 0.432MPa are adopted. These values are typical for soft polyurethane Hei-Cast 8000.
3. Numerical results 3.1 Poisson’s ratio First, the effects of various are considered, with and t being fixed. A close
CR IP T
inspection of the CMM in Fig. 1 reveals that when is greater than 1 the unit cell displays a re-entrant feature and is thus expected to give rise to an auxetic behavior, i.e., negative Poisson’s ratio. Bertoldi et al. (2010) showed that the buckled phase of a metamaterial with circular holes is very similar to the CMM in Fig. 2(a) with 1 and is indeed auxetic. When 1 , on the contrary, the Poisson’s ratio of the CMM is
AN US
positive. Between them (i.e., 1 ), the Poisson’s ratio becomes negligible. The sign dependence of the Poisson’s ratio is illustrated in the phase diagram in Fig. 2(a). Figure 2(b) shows the FE simulated progressive Poisson’s ratio (i.e., v x y ) as a function of the applied compressive strain y for three selected values of 0.5 , 1,
M
and 1.5, corresponding to the two regimes and the boundary between them in Fig. 2(a). It is clear that the simulated Poisson’s ratios are positive, close to zero, and
ED
negative, consistent with the phase diagram in Fig. 2(a).
(b)
0
1 ,
CE
1
=0.5 =1.0 =1.5
0
-1
1, 0
AC
2
-x/y
PT
a2
(a) 1, 0
0.00
a1
0.05
y
0.10
Fig. 2 (a) Phase diagram for the Poisson’s ratios of CMMs and (b) FE simulated strain trajectories for CMMs under uniaxial compression with a1 3 mm, 2 , t 0.5 and selected values of 0.5, 1.0 and 1.5 .
6
ACCEPTED MANUSCRIPT
y / Ei
Abcdefghijklm 1234567890 0.03 0.09 =0.5 =1.0 =1.5 0.06
0.03
0.03
0.06
y
0.09
CR IP T
0.00 0.00
Fig. 3 FE simulated uniaxial compressive stress-strain curves of CMMs with
2.0, t 0.5 and three different values of 0.5 , 1.0 and 1.5 .
AN US
3.2 Deformation pattern
The associated macroscopic stress-strain curves for 0.5 , 1.0, and 1.5 are presented in Fig. 3, where the stress is normalized by the corresponding effective Young’s modulus Ei (subscript i referring to different models with 0.5 , 1.0, and 1.5) in the vertical direction. By doing so, the initial elastic responses of the three
M
curves collapse into a single master curve. One can see from Fig. 3 that, within the shown strain range, the stress varies monotonically with strain for 1 while
ED
softening is evident for 1. To further explore the deformation mechanisms responsible for the difference in the curves in Fig. 3, we show the corresponding
PT
unreformed and deformed ( y 10% ) configurations in Fig. 4. It can be seen that the deformed patterns for both 0.5 and 1.5 in Fig. 4(a-b) are homogenous. For 1
CE
(Fig. 4(c)), however, three buckling induced localized deformation zones are present. It has been generally accepted that deformation localization in materials and
AC
structures is usually sensitive to various stochastic imperfections either inherently or externally induced. As a result, deformation localization usually occurs in a somewhat random manner and prediction of its precise location is difficult. However, Fig. 4(c) shows a striking feature that the localized zones are highly ordered instead of randomly distributed. We have varied the ratio , all showing uniform deformation for 1 up to moderate compression (e.g., 20% in strain). For 1, highly ordered localized deformation is achieved if other geometrical parameters are properly chosen. Moreover, FE results show that the feature of the transition from uniform to
7
ACCEPTED MANUSCRIPT localized deformation in Fig. 4 persists, providing that N y is no more than a fraction of N x . Results are not shown for brevity. The condition of N y
N x implies that the
CMMs be long strips which are of particular interest in this study.
Undeformed
Deformed
0
(a)
CR IP T
0.5 0
(b)
1.5
0
(c)
AN US
1.0
Fig. 4 FE results of the undeformed and deformed configurations of CMMs with (a) deformed configurations is 10%.
ED
4. Theoretical analysis
M
0.5 , (b) 1.5 , and (c) 1.0 . The macroscopic compressive strain for the
In this section, a simplified theoretical model is developed to uncover the underlying
PT
mechanisms responsible for the simulated uniform and ordered localized deformations shown in Fig. 4. Day et al. (1992) studied the elastic moduli of thin
CE
plates with circular holes and found that, when the size of the necks between adjacent holes is small, deformation and elastic strain energy are concentrated in the neck
AC
region and local deformation at the necks dominates the entire plate. In the deformed configurations of the CMMs, when the neck width t compared to other geometrical parameters is small, localized neck deformation prevails. Accordingly, a simplified unit cell model with rods connected by torsional springs can be derived from the dotted lines superimposed on the original unit cell shown in Fig. 2(a). A rod-spring model subject to vertical compression of displacement u y 2 is given in Fig. 5. The vertical rods shown in blue have an initial inclined anticlockwise angle 0 and initial length L0 . Note that the sign of 0 corresponds to whether is less than, equal to or 8
ACCEPTED MANUSCRIPT greater than 1. The axial stiffness of the elastomer rods can be approximated by a linear spring with stiffness K , the bending stiffness of the yellow and blue necks are denoted by the stiffness C1 and C2 of the torsional springs, while the shearing stiffness of blue necks is represented by Cs .
u y 2
= rod
tension
=
C1
neck
bending stiffness
L0
=
C2 bending stiffness
AN US neck
stiffness
CR IP T
0
K
shear
Cs stiffness
Fig. 5 Schematic of the rod-spring model and the corresponding stiffness for the unit
M
cell of the CMM. The compressive displacement is denoted by u y 2 .
ED
With the imposed vertical compression u y 2 and ignoring asymmetric modes, the strain energy of the unit cell rod-spring model can be expressed as (4)
PT
U cell , L 4C1 ( 0 )2 2C2 ( 0 )2 2K ( L L0 )2
CE
where and L are the current inclined angle and length of the elastic rods, respectively, and are related to each other by L (L0 cos 0 )sec . Defining the
AC
rotation angle as
0 and taking approximation of sec( 0 ) in terms of
Taylor’s series for small , Eq. (4) can be rewritten up to O( 5 ) as U cell g0 g1 g2 2 g3 3 g4 4
with
9
(5)
ACCEPTED MANUSCRIPT g 0 2 f 0 K H 0 f 0 H 0 2 L0 g1 4 f1 K H 0 f 0 H 0 L0
g 2 2 2C1 C2 2 H 0 K f12 2 f 0 f 2 H 0 2 f 2 L0
(6)
g3 2 K 2 f3 L0 2 f1 f 2 f 0 f 3 H 0 H 0
g 4 2 K 2 f 4 L0 f 22 2 f1 f 3 2 f 0 f 4 H 0 H 0
CR IP T
where H 0 L0 cos 0 and fi (i 0,.., 4) are the coefficients of Taylor’s series expansion of sec( 0 ) (e.g., f0 sec0 , f1 sec0 tan 0 , etc.).
Recall that the deformations in Fig. 4 are independent of N y and remains so when
N y is no more than a fraction of N x . A quasi one dimensional (1D) rod-spring model
AN US
with N x unit cells in the x direction can thus be constructed for the CMM (see, Fig. 6(a)). Under uniform vertical compression u y 2 , the characteristic rotation angle of the n-th unit cell can be denoted by n . Based on the assumption of symmetric deformations within each unit cell, the strain energy of the 1D system can be obtained
M
as
ED
U system Cs d 02 (sin n 1 sin n ) 2 2C2 n2 U cell n n
= Cs d 02 (sin n 1 sin n )2 U eff n ,
(7)
PT
n
where Cs is the shearing stiffness of the torsional spring between adjacent unit cells
CE
and U eff n g0 g1 n (2C2 g2 ) n2 g3 n3 g4 n4
is
the
effective
on-site
potential of the n-th unit cell. Equation (7) indicates that the actual pattern of the
AC
deformed metamaterials is the result of the competition between the interaction energy Cs d02 (sin n1 sin n )2 and the on-site potential U eff . A spring chain model corresponding to Eq. (7) is shown in Fig. 6(b). Applying the principle of minimum potential energy to Eq. (7), i.e., U system n 0 , one gets the following equilibrium equation for the n-th unit cell of the 1D system in Fig. 6(b), i.e.,
10
ACCEPTED MANUSCRIPT
2Cs d02 (2 n n1 n1 +h2 )
U eff n n
0
(8)
with h2 sin 20 n21 2 n1 n 6 n2 2 n1n 12n / 4 . 4.1 Macroscopic homogenous deformation The FE results in Section 3 show that, under the vertical compression, the CMMs with
CR IP T
1 (i.e., the initial inclined angle 0 0 ) are bending-dominated and undergo macroscopically homogenous deformation. With the developed theory (8), we will uncover the deformation mechanisms associated with 0 0 and 0 0 . Without loss of generality, we first consider 0 0 . In this case, during the initial loading, the on-
satisfying U cell
AN US
site potential of unit cell U cell is monostable with a single well at cell 0 , indicating that there is only one stable state for the cell
unit cells. Therefore, all unit cells in the metamaterial operate at the bottom of the stable wells during deformation. Although it is hard to get the general solution of Eq.
M
(8), its trivial solution, i.e., constant rotation with n =eff (n 1,..., N ) can be obtained through U eff 0 , yielding
ED
eff
PT
eff
g3 8 g 2 g 4 3g32 Q , 5/3 4 g4 2 g 4Q 3 27/3 g 4
(9)
CE
where Q is given by
AC
Q 216 g 2 s g3 g 4 54 g33 432 g1 g 42 13
4(24 g 2 s g 4 9 g32 )3 (216 g 2 s g3 g 4 54 g33 432 g1 g 42 ) 2
(10)
with g2 s (2C2 g2 ) . The trivial constant rotation solution (9) indicates a macroscopically homogenous deformation. 4.2 Ordered localized deformation When 0 0 , on the contrary, the metamaterial under vertical compression shows a
11
ACCEPTED MANUSCRIPT nontrivial feature of clusters of ordered localized deformation zones (Fig. 4(c)). Most importantly, this feature cannot be represented by the response of a unit cell under the same loading condition. It is in fact a type of low-energy excitation above the lowest energy state analyzed in Section 4.1. In view of the 1D spring chain model (8), we will show in the following that the excited ordered localization zones are analogous to the domain walls in phase transformational ferromagnetic and ferroelectric ceramics (Fatuzzo and Merz, 1967; Remoissenet, 2010) and can be regarded as static
With 0 0 , Eq.(8) can be simplified as
U eff n n
0
AN US
2Cs d02 (2 n n 1 n 1 )
CR IP T
topological solitons (Deng et al., 2017; Zhang et al., 2018).
with
(11)
2
(12)
M
2 U eff n =U 0 n2 2 K 2 U0 U0
where U 0 K (3L0 8 )( L0 ) 6 0 and 6 K ( L0 ) (2C1 C2 ) . It is
ED
noted that the last two terms on the right-hand side of Eq. (12) are independent of n and will be omitted in the following derivation. Clearly, the effective on-site potential
AC
and
CE
PT
U eff n is a double-well potential (i.e., bi-stable) under the condition of 0 , i.e.,
KL20 8(C1 C2 )
(13)
c m
(14)
with c L0 L20 8(C1 C2 ) / K
2 and
m
8L0 3 . If the conditions (13) and
(14) are not satisfied, U eff n is monostable, as illustrated in the phase diagram in Fig. 6(c) with c .
12
ACCEPTED MANUSCRIPT (c)
(a)
(b) unit cell
interaction
KL02
…
bi-stable
on-site potential
1
…
8 monostable
CR IP T
C1 C2
Fig. 6 (a) The simplified rod-spring model of the CMM, (b) its equivalent 1D spring chain model, and (c) the phase diagram of the on-site potential with single and double wells for c m .
Consider the case that the metamaterial in Fig. 6 satisfies Eq. (13). Under small
AN US
loading (i.e., c ), U eff n is monostable at n 0 and the metamaterial undergoes uniform deformation. Increasing the compression beyond c , U eff n transits from monostable to bistable and each constituent unit cells in the metamaterial have two stable states. Under such a circumstance, instability and localization may
M
occur and the deformed pattern is governed by the competition between the double
ED
well on-site potential and the interaction energy between adjacent unit cells. When N x in Fig. 6, a continuum field to describe the evolution of the
PT
characteristic angle n can be derived with the following quantities:
CE
n (x n ) 2 2 (n ) (n ) ..., x 2 x 2
AC
n 1 (t ) [(n 1) ] (n )
where
is the lattice parameter in the x direction and
(15)
2(1 t ) . Based on Eq.
(15), the continuum limit of the equilibrium equation (11) can be obtained as
2 V 0. x 2 with
13
(16)
ACCEPTED MANUSCRIPT 1 2U 0 2 V 4 Cs d02 2 U0
2
2 1 2 2 . 4
(17)
which is the well-known Klein-Gordon equation (Remoissenet, 2010). Without loss of generality, taking x 0 0 , multiplying Eq. (16) by d u d x and integrating with respect to x yield
2 2 2 ( ) 0 x 2
/ 2a2
(18)
0 x 2V ( ) x 0
.
/ 2 2 and
/ 2 2 . Then, x and in Eq.
2
/ 2b2
Define
two
AN US
where
CR IP T
2
quantities
(18) can be related to each other through the elliptic integral
d
0
x
(a 2 2 )(b2 2 )
/ 2 d x / 2x 0
(19)
ED
(Remoissenet, 2010), i.e.,
M
the solution of which can be expressed by the elliptic integral of the first kind
a sn b / 2 x; m
CE
PT
with the period
AC
where m
P
4 b /2
EllipticK m /
(20)
(21)
a 2 b2 is the modulus of the elliptic sine function sn x; m , EllipticK m
is the elliptic integral of modulus m , a 2m (1 m) , and b a
m . Once the
period P is determined, m can be evaluated for a given problem. Formula (20) is the kink-antikink solution to the 1D Klein-Gordon system (16) and represents a periodic topological solitary excitation (Zhang et al., 2018). It is noted that Eq. (20) depends upon the stiffness of the rods and necks K , C1 , C2 and Cs . In turn, these stiffness are determined by the geometrical parameters and the constituent 14
ACCEPTED MANUSCRIPT solid of the CMM in Fig. 1. As will be shown later, the non-trivial solution (20) actually describes the ordered localized deformations of the metamaterial in Fig. 4(c). 4.3 Spring stiffness of the simplified rod-spring model Day et al. (1992) developed simplified rod-spring models to calculate the effective tension stiffness K , shear and bending stiffness ( Cs , C1 and C2 ) of the wall necks in
CR IP T
the plates with circular holes. They focused on thin necks (i.e., small t ). Coulais (2016) extended the model for elliptical holes and straight beams were replaced by curved beams when is larger than 1. Both neck width t and ellipse ratio were discussed. In order to get analytical results, two limiting cases such as extremely small and large t were considered. It should be pointed out that the model considered
AN US
by Coulais (2016) is different from the CMM shown in Fig. 1.
When the neck width t or the ellipse ratio is finite, however, it is difficult to get the analytical expressions of K , C1 , C2 and Cs for the considered CMMs. Here, we conduct extensive FE simulations to quantify the dependence of the stiffness upon t
M
and with 1 . A semi-empirical expression of the dependence is assumed to have
ED
the following form of
y , t G B t
(22)
PT
where y represents either K , C1 , C2 or Cs , G is the parameter for
CE
nondimensionalization of y , and B and are the functions to be determined by curve fitting the corresponding FE results. Therefore, the dimensionless stiffness
y , t
G . With dimensional analysis, G can be
AC
can be defined as y , t
expressed as G = ( 1 t ) , and a12 for tension, shear and bending stiffness, respectively. Based upon the unit cell model in Fig. 1 and the selected material constants in Eq. (3), linear perturbation analysis is used to calculate the tension stiffness K and shear stiffness Cs by applying uniaxial tension and simple shear, respectively. The periodic boundary conditions suggested by Chen et al. (1999) are adopted to ensure deformation compatibility. Figure 7(a) shows the uniaxial tension model and the 15
ACCEPTED MANUSCRIPT definition of K F (2 w) with F and
being the associated load and
displacement. The FE calculated normalized tension stiffness K as a function of nondimensional neck width t is shown in Fig. 7(b) for various ellipse ratios . From the results in Fig. 7(b), one can extract B and in Eq. (22) for K (denoted by BK and K with subscript K referring to tension stiffness K ). They are
governed by two asymptotes for 1 and
1 , i.e., K ( ) 1 0.121 0.349
as 1 and K ( ) 0.00257 0.557 as
1 . Guided by the observation,
K can be assumed to have the form of
ED
3.3
1.2
0.1
(d)
BK()
1
0.50
3
4.0
1.9
PT
CE
K()
0.1
FE Eq.(23)
0.52
AC
(23)
1
t
FE Eq.(24)
0.6
0.54
2
1 M
2.6
1
K
M
K F 2 w
0.56
(b) 1
(c)
M
AN US
K ( ) tanh A 1 1 1 1 1 2 (a)
CR IP T
plotted in Fig. 7(c) and 7(d) as symbols. It can be seen from Fig. 7(c) that BK is
0.4 0.2 1
4
2
3
4
Fig. 7 Dependence of the tension stiffness K upon the geometrical parameters t and
with 1 . (a) FE model to calculate K , (b) FE calculated K as a function of t for various , (c) and (d) curve fittings of the semi-empirical parameters K and
16
ACCEPTED MANUSCRIPT BK for K . Fitting Eq. (23) against the FE results in Fig. 7(c) gives A 2.11 and M 32.2 . The fitted curve is given in Fig. 7(c) as solid line, agreeing well with the FE results. From Fig. 7(d), it is evident that curve fitting of BK for K gives BK 1.08 0.395 0.0425 2 .
CR IP T
(a)
(24)
(b)
4.0 3.3
Cs
0.1
2.6
0.01
0.1
(c)
FE Eq.(25)1
1.2
1
t
(d)
FE Eq.(25)2
0.54
M
Bs()
1.3 1.2 1.1 2
3
0.53 4
1
PT
1
ED
s()
1.9
AN US
Cs F 2 w
2
3
4
CE
Fig. 8 Dependence of stiffness Cs upon the geometrical parameters t and with
1 . (a) FE model for calculating Cs , (b) FE calculated Cs as a function of t for
AC
various , (c) and (d) curve fittings of the semi-empirical parameters s and
Bs for Cs .
The FE model for calculating the shear stiffness Cs for various geometrical parameters is illustrated in Fig. 8(a). The obtained FE results are given in Fig. 8(b), showing that the shear stiffness Cs can also be approximated by Eq. (22) with curve fitted parameters
17
ACCEPTED MANUSCRIPT s ( ) 1.45 0.0957 ,
(25)
Bs 0.552 0.0119 0.00163 .
where subscript s denotes shear stiffness. At 1 , s ( ) 1.35 , very close to the results of circular hole obtained by Day et al. (1992).
(a)
(b)
C1 M w
3.3
CR IP T
C1
0.1
2.6
0.01
1.9
0.001
0.1 (d)
FE Eq.(27)1
t
2.3
2
3
4
ED
1
M
2.2
1.2
FE Eq.(27)2
0.5
B1()
1()
2.4
2.1
1
AN US
(c) 2.5
4.0
0.4 0.3
1
2
3
4
Fig. 9 Dependence of stiffness C1 upon the geometrical parameters t and with
PT
1 . (a) FE model to calculate C1 ; (b) FE calculated C1 as a function of t for
CE
various ; (c) and (d) curve fittings of the semi-empirical parameters 1 and
AC
B1 for C1 .
To calculate the bending stiffness C1 and C2 , the FE models given in Figs. 9(a) and 10(a) are used. The models are analogous to that of Day et al. (1992). The FE calculated dependence of C1 and C2 upon t with selected values are shown in Figs. 9(b) and 10(b), respectively. Again, fitting the numerical results against formula (22) gives
18
ACCEPTED MANUSCRIPT C1 G B1 t 1
(26)
C2 G B2 t 2 with
1 2.63 0.154
(27)
and
2 2.58 0.102 B2 0.592 0.0729
CR IP T
B1 0.608 0.0747
(28)
AN US
for the considered range of t [0.05, 0.7] and [1.2, 4] . As expected, when 1 ,
1 and 2 approach the theoretical results given by Day et al. (1992) for circular holes, i.e., 2.5.
3.3
0.01
2.6
0.001
1.9 0.1
(d)
FE Eq.(28)1
CE
(c) 2.5
AC
1
t
2.1
1.2
FE Eq.(28)2
0.5
B2()
2.3
4.0
0.1
C2
ED
C2 M w
PT
2()
(b)
M
(a)
0.4 0.3
1
2
3
4
1
2
3
4
Fig. 10 Dependence of stiffness C2 upon the geometrical parameters t and with
1 . (a) FE model to calculate C2 ; (b) FE calculated C2 as a function of t for various ; (c) and (d) curve fittings of the semi-empirical parameter 2 and 19
ACCEPTED MANUSCRIPT B2 for C2 .
5. Comparison between FE and theoretical results We now compare the developed theoretical predictions in Section 4 with the FE simulated results. From the FE calculated deformed configurations in Fig. 4, the characteristic rotation for each unit cell located at the middle row can be extracted
CR IP T
using the definition indicated in the simplified model in Fig. 5. The corresponding FE results on the variation of along the x direction are shown in Fig. 11 as symbols for 0.5 , 1, and 1.5. As can be seen, the characteristic rotations for 0.5 and 1.5 are almost constant, a direct consequence of the homogeneous deformations in Fig. 4(a-b). When 1 , on the contrary, is not a constant and its sign alternatively
AN US
changes along x . Note that is the relative rod rotation. A negative with 1 indicates clockwise rotation of the originally vertical rods and a localized zone consisting of re-entrant microstructures. The extracted constant for 0.5 and 1.5 and regular variation of with 1 correlate reasonably well the FE deformation patterns in Fig. 4, showing that the defined characteristic rotation is a
M
suitable parameter to quantify the deformation of the CMM under compression.
ED
For the purpose of comparison, the corresponding theoretical predictions obtained in Section 4 are also included in Fig. 11 as lines. When 0.5 and 1.5, the effective
PT
on-site potential U eff in Eq. (8) is monostable and the trivial homogeneous deformation solution (9) applies. Specifically, Eq. (9) gives eff 0.159 for
CE
0.5 , with the stiffness K 0.934MPa , C1 0.621N , and C2 0.590N obtained from FE simulations. Similarly, Eq. (9) gives eff 0.183 for 1.5 , with
AC
K 1.47MPa , C1 1.17N , and C2 1.13N .
For the case of 1 , the effective on-site potential U eff in Eq. (8) becomes bistable when a moderate compression (e.g., 10%) is applied. Then, the kink-antikink solution (20) is adopted, with the stiffness K 1.34MPa , C1 1.19N , C2 1.10N ,
Cs 0.279MPa , and m 0.35 . As can be seen from Fig. 11, the theoretical predictions for 0.5 and 1.5 agree well with the FE results. For 1 , the overall agreement between the theoretical solution and the FE results is also very good: both 20
ACCEPTED MANUSCRIPT the amplitude and the spacing of the localized deformation zones are captured. Note that the kink-antikink solution is for periodic systems while the system in Fig. 4c is of finite size. It is thus not surprise that, for 1 , the theoretical kink-antikink solution (20) slightly deviates from the FE results at the two ends. =0.5
Theory
=1.0
=1.5
0.4 0.0 -0.4
=0.5
FE -5
-0.8 -10
=1.0 0 x/l
(a)
CR IP T
Rotation angle
0.8
=1.5 5
10
Fig. 11 Comparison of the FE and theoretical predictions of the characteristic rotation
AN US
angle for the deformed CMMs shown in Fig. 4(c). (a) (b)
M
0.50
Theory
0.25
ED
0.00 -0.25 -10
PT
Rotation angle
(c)
-5
CE
Fig. 12 A quasi 1D CMM
FE
0 x/l
5
10
with the geometry parameters a1 3.0 mm, 1 ,
t 0.55 and 2.33 . (a) The initial configuraiton, (b) the ordered localized
AC
deformation under compressive strain 9.7% , and (c) the theoretical kink-antikink solution and FE predicted variation of the characteristic rotation angle. To further validate the theoretical analysis in conjunction with the determined semiempirical results in Section 4, a quasi-one dimensional CMM with N x N y 20 1 and different geometrical parameters 1 , t 0.55 and 2.33 is considered. The particular choice of the geometrical parameters is not used in determining the semiempirical results in Eqs. (23-28). The undeformed and deformed ( 9.7% )
21
ACCEPTED MANUSCRIPT configurations are shown in Fig. 12(a-b), respectively. Regularly spaced localized deformation zones are present in the deformed shape. From the semi-empirical formulae (Eqs. (23-28)), one has K 1.30MPa , Cs 0.340MPa , C1 1.40N ,
C2 1.40N , and m 0.45 . Accordingly, Eq. (20) predicts an ordered variation of characteristic rotation angle and is plotted in Fig. 12(c) as the solid line, showing good
(b) B C
0.04
Rotation angle
A
0.00 0.00
0.04
y
D
0.08
Theory
FE -5
AN US
A
0.8 0.4 0.0 -0.4 -0.8 -10
D
B C
0.02
(c)
0
A
A
B
B
M
y / E
(a)
CR IP T
agreement with the FE result.
C
0 x/l
C
D
D
5
10
ED
Fig. 13 Deformation evolution of a CMM ( a1 3.0 mm, 1 , t 0.3 and 1.5 ) under uniaxial compression. (a) The FE predicted normalized compressive stress-
PT
strain curve, (b) the original and FE simulated deformed configurations at selected loading levels marked in the stress-strain curve, and (c) the theoretical and FE
CE
predicted characteristic rotation angle. It is noted that the periodic kink-antikink solution (20) depends not only upon the
AC
geometrical parameters but also the magnitude of compression (i.e., ). To further validate the theoretical analysis in Section 4, a quasi-one dimensional CMM with
N x N y 20 1 and geometrical parameters 1 , t 0.3 and 1.5 is
considered. Again, the particular choice of the geometrical parameters is not used in determining the semi-empirical results (Eqs. (23-28)). The FE simulated compressive stress-strain curve and the deformation configurations at selected load levels (A–D) are shown in Fig. 13(a-b), respectively. For small deformations (e.g., Point A in the linear region), the CMM undergoes uniform deformation. With the load increasing, 22
ACCEPTED MANUSCRIPT the stress-strain curve gradually deviates from its linear response. Upon reaching Point B (slightly before the peak load), localized deformation zones start to develop. When the loading increases beyond the critical strain at the peak, strain softening accompanied by localized deformation zones can be observed (Point C). Further increasing the strain from C to D, minor strain hardening can be observed. Interestingly, there is little change in the locations of the localized zones, showing that
CR IP T
deformation localized zones are robust. With Eq. (20), the characteristic rotation angles of the quasi-one dimensional CMM at different loading levels (i.e., A–D) can be theoretically predicted ( m 0.17, 0.41, 0.53 and 0.79 for A–D, respectively) and are plotted in Fig. 13(c) as lines, together with the FE results denoted by symbols. It is clear that the FE simulated evolution of the localization zones as characterized by the
AN US
amplitude and spacing of the rotation angle is faithfully captured by the developed theoretical model.
6. Concluding remarks
A two dimensional CMM comprising of elliptical holes of dual sizes is designed.
M
Systematical FE simulations are conducted to explore its responses to uniaxial compression. It is shown that the CMM can have positive, negative, or zero Poisson’s
ED
ratio and can also be monostable or bistable, depending upon its geometrical parameters. Under moderate compression, the bistable CMM show ordered and highly
PT
localized deformation zones in the direction perpendicular to the compression, in
CE
sharp contrast to the homogeneous deformation of the monostable CMM. A simplified rod-spring model is derived from the unit cell of the CMM to uncover the deformation mechanisms and a theoretical spring chain model is developed.
AC
Dependence of the model parameters upon the geometrical parameters of the CMM is quantified by FE simulations. In the continuum limit, the theoretical model reduces to the Klein-Gordon equation. Analytical trivial solution (i.e., constant rotation) and kink-antikink solution are obtained, corresponding to the models with monostable and bistable potentials, respectively. The constant rotation solution indicates homogeneous deformation and correlates well to the FE simulated deformed configurations of the monostable CMM. The kink-antikink solution faithfully captures the features (i.e., periodic spacing and magnitude) of the ordered localization zones of the bistable 23
ACCEPTED MANUSCRIPT CMM under compression. Note that the feature of kink-antikink pairs is ubiquitous (e.g., domain walls in phase transform materials and multi-solitary excitation in waves). The excellent agreement between the FE simulations and theoretical predictions implies that the ordered deformation localization zones are similar to multi-soliton in fluids and physics. This study can thus provide insight into the excitation of ordered localized deformation, phase transformation and static
CR IP T
topological soliton, as well as the design of mechanical metamaterials. Acknowledgements This work is supported by the National Natural Science Foundation of China (Nos. 11732007 and 11472149) and the Science Challenge Project of China (Grant No. TZ2018007).
AN US
References
An, S., Kim, B., Kwon, S., Moon, G., Lee, M., Jhe, W., 2018. Bifurcation-enhanced ultrahigh sensitivity of a buckled cantilever. Proc. Natl. Acad. Sci. 115, 2884– 2889. https://doi.org/10.1073/pnas.1716067115
M
Ashby, M., Evans, A., Fleck, N., Gibson, L., Hutchinson, J., Wadley, H., 2000. Metal Foams: A Design Guide. Butterworth-Heinemann, Oxford.
ED
Audoly, B., Hutchinson, J.W., 2016. Analysis of necking based on a one-dimensional model. J. Mech. Phys. Solids 97, 68–91. https://doi.org/10.1016/j.jmps.2015.12.018
PT
Banhart, J., 2001. Manufacture, characterisation and application of cellular metals and metal foams. Prog. Mater. Sci. 46, 559–632. https://doi.org/10.1016/S00796425(00)00002-5
CE
Bažant, Z.P., Cedolin, L., 2010. Stability of structures: elastic, inelastic, fracture and damage theories. World Scientific Pub, Hackensack, NJ.
AC
Bertoldi, K., Boyce, M., 2008. Mechanically triggered transformations of phononic band gaps in periodic elastomeric structures. Phys. Rev. B 77, 052105. https://doi.org/10.1103/PhysRevB.77.052105 Bertoldi, K., Reis, P.M., Willshaw, S., Mullin, T., 2010. Negative Poisson’s ratio behavior induced by an elastic instability. Adv. Mater. 22, 361–366. https://doi.org/10.1002/adma.200901956 Cao, B., Wu, G., Xia, Y., Yang, S., 2016. Buckling into single-handed chiral structures from pH-sensitive hydrogel membranes. Extreme Mech. Lett. 7, 49–54. https://doi.org/10.1016/j.eml.2015.12.011 Chen, C., Lu, T.J., Fleck, N.A., 1999. Effect of imperfections on the yielding of twodimensional foams. J. Mech. Phys. Solids 47, 2235–2272. 24
ACCEPTED MANUSCRIPT https://doi.org/10.1016/s0022-5096(99)00030-7 Chen, H., Nassar, H., Huang, G.L., 2018. A study of topological effects in 1D and 2D mechanical lattices. J. Mech. Phys. Solids. https://doi.org/10.1016/j.jmps.2018.04.013 Chen, Y., Li, T., Scarpa, F., Wang, L., 2017. Lattice metamaterials with mechanically tunable poisson’s ratio for vibration control. Phys. Rev. Appl. 7, 024012.
CR IP T
Coulais, C., 2016. Periodic cellular materials with nonlinear elastic homogenized stress-strain response at small strains. Int. J. Solids Struct. 97–98, 226–238. https://doi.org/10.1016/j.ijsolstr.2016.07.025 Day, A.R., Snyder, K.A., Garboczi, E.J., Thorpe, M.F., 1992. The elastic moduli of a sheet containing circular holes. J. Mech. Phys. Solids 40, 1031–1051.
AN US
Deng, B., Raney, J.R., Tournat, V., Bertoldi, K., 2017. Elastic vector solitons in soft architected materials. Phys. Rev. Lett. 118, 204102. https://doi.org/10.1103/PhysRevLett.118.204102 Deshpande, V., Ashby, M., Fleck, N., 2001. Foam topology: bending versus stretching dominated architectures. Acta Mater. 49, 1035–1040.
M
Diamant, H., Witten, T.A., 2011. Compression induced folding of a sheet: an integrable system. Phys. Rev. Lett. 107, 164302. https://doi.org/10.1103/PhysRevLett.107.164302 Fatuzzo, E., Merz, W.J., 1967. Ferroelectricity. North-Holland Pub. Co., Amsterdam.
ED
Gibson, L.J., Ashby, M.F., 1997. Cellular solids: Structure and properties. Cambridge University Press, Cambridge.
PT
Gibson, L.J., Ashby, M.F., Harley, B.A., 2010. Cellular materials in nature and medicine. Cambridge University Press, Cambridge.
CE
Jiang, Y., Liu, Z.Y., Matsuhisa, N., Qi, D., Leow, W.R., Yang, H., Yu, J., Chen, G., Liu, Y., Wan, C., Liu, Z.J., Chen, X., 2018. Auxetic mechanical metamaterials to enhance sensitivity of stretchable strain sensors. Adv. Mater. 30, 1706589. https://doi.org/10.1002/adma.201706589
AC
Kadic, M., Bückmann, T., Schittny, R., Wegener, M., 2013. Metamaterials beyond electromagnetism. Rep. Prog. Phys. 76, 126501. https://doi.org/10.1088/00344885/76/12/126501 Kane, C.L., Lubensky, T.C., 2014. Topological boundary modes in isostatic lattices. Nat. Phys. 10, 39–45. https://doi.org/10.1038/nphys2835 Kang, S.H., Shan, S., Noorduin, W.L., Khan, M., Aizenberg, J., Bertoldi, K., 2013. Buckling-induced reversible symmetry breaking and amplification of chirality using supported cellular structures. Adv. Mater. 25, 3380–3385. https://doi.org/10.1002/adma.201300617
25
ACCEPTED MANUSCRIPT Kochmann, D.M., Bertoldi, K., 2017. Exploiting Microstructural Instabilities in Solids and Structures: From Metamaterials to Structural Transitions. Appl. Mech. Rev. 69, 050801-050801–24. https://doi.org/10.1115/1.4037966 Kolken, H.M.A., Zadpoor, A.A., 2017. Auxetic mechanical metamaterials. RSC Adv. 7, 5111–5129. https://doi.org/10.1039/C6RA27333E Li, S., Zhao, D., Niu, H., Zhu, X., Zang, J., 2018. Observation of elastic topological states in soft materials. Nat. Commun. 9, 1370. https://doi.org/10.1038/s41467-018-03830-8
CR IP T
Liu, X.N., Hu, G.K., Huang, G.L., Sun, C.T., 2011. An elastic metamaterial with simultaneously negative mass density and bulk modulus. Appl. Phys. Lett. 98, 251907. https://doi.org/10.1063/1.3597651 Liu, Z., Zhang, X., Mao, Y., Zhu, Y.Y., Yang, Z., Chan, C.T., Sheng, P., 2000. Locally resonant sonic materials. Science 289, 1734–1736.
AN US
Lu, T.J., Valdevit, L., Evans, A.G., 2005. Active cooling by metallic sandwich structures with periodic cores. Prog. Mater. Sci. 50, 789–815. https://doi.org/10.1016/j.pmatsci.2005.03.001 Maznev, A.A., Every, A.G., Wright, O.B., 2013. Reciprocity in reflection and transmission: What is a `phonon diode’? Wave Motion 50, 776–784. https://doi.org/10.1016/j.wavemoti.2013.02.006
ED
M
Meza, L.R., Das, S., Greer, J.R., 2014. Strong, lightweight, and recoverable threedimensional ceramic nanolattices. Science 345, 1322–1326. https://doi.org/10.1126/science.1255908
PT
Mullin, T., Deschanel, S., Bertoldi, K., Boyce, M., 2007. Pattern transformation triggered by deformation. Phys. Rev. Lett. 99, 084301. https://doi.org/10.1103/PhysRevLett.99.084301
CE
Overvelde, J.T.B., Bertoldi, K., 2014. Relating pore shape to the non-linear response of periodic elastomeric structures. J. Mech. Phys. Solids 64, 351–366. https://doi.org/10.1016/j.jmps.2013.11.014
AC
Papka, S.D., Kyriakides, S., 1999. In-plane biaxial crushing of honeycombs—: Part II: Analysis. Int. J. Solids Struct. 36, 4397–4423. https://doi.org/10.1016/S0020-7683(98)00225-X Pocivavsek, L., Dellsy, R., Kern, A., Johnson, S., Lin, B., Lee, K.Y.C., Cerda, E., 2008. Stress and fold localization in thin elastic membranes. Science 320, 912–916. Rafsanjani, A., Akbarzadeh, A., Pasini, D., 2015. Snapping mechanical metamaterials under tension. Adv. Mater. 27, 5931–5935. https://doi.org/10.1002/adma.201502809 Remoissenet, M., 2010. Waves called solitons: concepts and experiments. Springer Science & Business Media., Berlin. 26
ACCEPTED MANUSCRIPT Rivlin, R.S., 1948. Large elastic deformations of isotropic materials IV. Further developments of the general theory. Philos. Trans. R. Soc. Math. Phys. Eng. Sci. 241, 379–397. https://doi.org/10.1098/rsta.1948.0024 Schaedler, T.A., Jacobsen, A.J., Torrents, A., Sorensen, A.E., Lian, J., Greer, J.R., Valdevit, L., Carter, W.B., 2011. Ultralight metallic microlattices. Science 334, 962–965. https://doi.org/10.1126/science.1211649
CR IP T
Smith, D.R., Pendry, J.B., Wiltshire, M.C.K., 2004. Metamaterials and negative refractive index. Science 305, 788–792. https://doi.org/10.1126/science.1096796 Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389–407. https://doi.org/10.1007/BF00036191 Wu, G., Xia, Y., Yang, S., 2014. Buckling, symmetry breaking, and cavitation in periodically micro-structured hydrogel membranes. Soft Matter 10, 1392– 1399. https://doi.org/10.1039/C3SM51640G
AN US
Yang, D., Jin, L., Martinez, R.V., Bertoldi, K., Whitesides, G.M., Suo, Z., 2016. Phase-transforming and switchable metamaterials. Extreme Mech. Lett. 6, 1– 9. https://doi.org/10.1016/j.eml.2015.11.004
AC
CE
PT
ED
M
Zhang, Y., Li, B., Zheng, Q.S., Chen, C.Q., 2018. Static periodic topological soliton in mechanical metamaterials. Under preparation.
27