Toughening in a nacre-like soft-hard layered structure due to weak nonlinearity in the soft layer

TL;DR

This paper demonstrates that weak nonlinearity in the soft layers of nacre-like structures significantly enhances toughness by redistributing stress, providing insights for designing tough bio-inspired composites.

Contribution

It introduces an analytical and numerical model showing how weak nonlinearity in soft layers improves toughness in layered structures, a novel mechanism for bio-inspired material design.

Findings

Weak nonlinearity redistributes stress at crack tips.

Enhanced toughness observed due to nonlinear soft layers.

Model provides a design principle for tough composites.

Abstract

Recently, it has been found experimentally that hydrated nacre exhibits a nonlinear mechanical response. While mechanical nonlinearity has been shown to be important in other biological structures, such as spider webs, the implications of mechanical nonlinearity in nacre have not been explored. Here, we show that the nonlinear mechanical response of nacre can be reproduced by an analytical model, which reflects a nacre-like layered structure, consisting of linear-elastic hard sheets glued together by weakly nonlinear-elastic soft sheets. We develop scaling analysis on this analytical model, and perform numerical simulations using a lattice model, which is a discrete counterpart of the analytical model. Unexpectedly, we find the weak nonlinearity in the soft component significantly contributes to enhancing toughness by redistributing the stress at a crack tip over a wider area. Beyond demonstrating a mechanism that explains the unusual properties of biological nacre, this study points to a general design principle for constructing tough composites using weak nonlinearity, and is useful as a guiding principle to develop artificial layered structures mimicking nacre.

Equations

$$d=d_s+d_h\simeq d_h,$$

$$\sigma$$

$$\epsilon =\frac{E_s}{E_h}{\left(\frac{d_h}{d_s}\right)}^{1/n_s}.$$

$${\sigma }_{yy}\sim E_h{e}_{yy}^{1/n_h}\text{ and }{\sigma }_{yx}\sim \epsilon E_h{e}_{yx}^{1/n_s}$$

$$G$$

$${\sigma }_F$$

$$\lambda \sim \frac{d}{a_0}{\epsilon }^{-\frac{1}{1+1/n_s}}{\left(\frac{G_h/E_h}{a_0}\right)}^{\frac{1}{1+n_h}+\frac{1}{1+1/n_s}-1}$$

$$\frac{{\sigma }_M}{E_h}\sim {\left(\epsilon {\left(\frac{{\sigma }_0}{E_h}\right)}^{\left(1+n_h\right)\left(1+1/n_s\right)}{\left(\frac{L}{d}\right)}^{1+1/n_s}\right)}^{\frac{1}{\left(1+n_h\right)+\left(1+1/n_s\right)}}$$

$$d/a_0=50,\frac{G_h/E_h}{a_0}=1/100,\frac{E_s}{E_h}=\frac{35}{65000},\frac{d_s}{d_h}=1/19.$$

$$F\sim E_h{e}_{yy}^{1/n_h+1}+\epsilon E_h{e}_{yx}^{1/n_s+1}$$

$$(u/Y{)}^{1/n_h+1}\sim \epsilon (u/X{)}^{1/n_s+1}$$

$$E_h(u/Y{)}^{1/n_h+1}aY\sim aG$$

$$\frac{{\sigma }_F}{E_h}\sim {\left(\epsilon {\left(\frac{G/E_h}{X}\right)}^{1+1/n_s}\right)}^{\frac{1}{2+n_h+1/n_s}}\sim {\left(\frac{G/E_h}{Y}\right)}^{\frac{1}{{}^{1+n_h}}}$$

$${\sigma }_F^{MN}\sim E_h{\left(\frac{G_h/E_h}{X_h}\right)}^{\frac{1}{1+n_h}}.$$

$${\sigma }_h\sim E_h{\left(\frac{G_h/E_h}{a_0}\right)}^{\frac{1}{1+n_h}}$$

$$\sigma (r)\sim {\sigma }_0(X/r{)}^{\eta }$$

$$\eta =\frac{1+1/n_s}{2+n_h+1/n_s},$$

$$\sigma (r)\sim E_h{\left(\epsilon {\left(\frac{G/E_h}{r}\right)}^{1+1/n_s}\right)}^{\frac{1}{2+n_h+1/n_s}}\text{.}$$

$$\frac{G}{G_h}=\frac{d}{a_0}{\epsilon }^{-\frac{1}{1+1/n_s}}{\left(\frac{G_h/E_h}{a_0}\right)}^{\frac{1}{1+n_h}+\frac{1}{1+1/n_s}-1}$$

$${\sigma }_0\sim E_h{e}_0^{1/n_h}\text{ and }{e}_0\sim u/Y\text{.}$$

$$X\sim {\left(\epsilon e_0^{-\left(1/n_h-1/n_s\right)}\right)}^{\frac{1}{1+1/n_s}}Y\equiv \tilde{\epsilon }Y$$

$$\frac{{\sigma }_F}{{\sigma }_F^{MN}}={\left(\frac{d}{a_0}{\epsilon }^{-\frac{1}{1+1/n_s}}{\left(\frac{G_h/E_h}{a_0}\right)}^{\frac{1}{1+n_h}+\frac{1}{1+1/n_s}-1}\frac{X_h}{Y}\right)}^{\frac{1}{1+n_h}}$$

$$\frac{{\sigma }_F}{{\sigma }_F^{MN}}={\left(\frac{d}{a_0}{\left(\frac{G_h/E_h}{a_0}\right)}^{\frac{1}{1+n_h}+\frac{1}{1+1/n_s}-1}{\left(e_0^{-\left(1/n_h-1/n_s\right)}\right)}^{\frac{1}{1+1/n_s}}\right)}^{\frac{1}{1+n_h}}$$

$${\sigma }_M/E_h\sim e_0^{1/n_h}(a/d{)}^{\eta }$$

Full text

I Introduction

Natural materials often exhibit remarkable hierarchical structures leading to outstanding mechanical characteristics GaoARMR10 ; Fratzl2007Review ; Gao2003PNAS ; OkumuraMRS2015 ; Hellmich2004Bone , as observed in bone, spider silk Vollrath2001NatureReview , and the exoskeletons of crustaceans Club2012ScienceA ; Lobster2010AM . Nacre, which is found in many seashells and protects shellfish from their environment, is composed of soft and hard layers, and has been studied as a prototype material for several decades Currey1977 ; JacksonVincentTurner ; Sarikaya . Researchers have been inspired by nacre’s remarkable structure to develop new materials that demonstrate excellent mechanical performance. KatoNacre2000AM ; Mayer2005Science ; Tomsia2008Science ; BioinspiredPolymerFilmScience2008 ; Science2007LayerPolymer ; bouville2014strong ; sakhavand2015universal .

There have been many studies discussing the mechanisms responsible for the toughness of nacre. Toughening mechanisms have been proposed based on a variety of experimental observations, such as (1) step-wise elongation Smith1999NatureNacre , (2) micro-cracking and crack bridging Kamat2000NatureNacre , (3) thin compressive layers Rao1999Science , (4) rough layer interfaces SuoEvansHuchinson , (5) mineral bridges SongBai2003 , and (6) wavy surface of the plates Barthelat2007JMPS . On the theoretical side, various approaches have been explored, including (1) elastic models Rao1999Science based on analytical solutions OkumuraPGG2001 ; Okumura2002nacre and on scaling arguments DeGennesOkumuraCR ; okumura2005JPCM , (2) viscoelastic models 2003parallel ; Okumura2002nacre , (3) micro-mechanical models Kotha2001Micromechanical , (4) numerical models including finite-element models Katti ; JiGao2004 ; Barthelat2007JMPS , (5) a fuse network model Nukala , and (6) a model with a periodic Young’s modulus fratzl2007periodic .

In general, materials start breaking from the tip of small cracks where stress is concentrated. Therefore, the reduction of the stress concentration at crack tips is key for material toughness Griffith ; Lawn ; Anderson . Using a simple linear model of nacre, and deriving an analytical expression of the stress and strain fields near a crack tip, we have shown that there is a significant reduction in stress near crack tips in nacre OkumuraPGG2001 ; Hamamoto2008PRE . This reduction in stress concentration was numerically confirmed using a simple network model Aoyanagi2009PRE , as well as finite-element calculations Hamamoto2013 . These numerical studies elucidated a physical picture for the mitigation of stress concentration, where enhanced elongation of soft layers effectively suppresses the deformation of the hard layer component, leading to a reduction in stress concentration. This suppression occurs effectively by the fact that the soft layers constitute only a small percentage of the bulk material, and thus stress in the bulk material is governed by the hard layers.

The mechanism of stress reduction above was established using linear models, but the importance of nonlinear response in biological materials has recently attracted greater interest. While we previously showed a simple linear model demonstrates the high mechanical adaptability of spider webs AoyanagiOkumura2010 , the nonlinear response of spider silk was shown to be a key factor in their mechanical superiority Buehler2012Nature . Previously, we have shown the mechanism of stress reduction in nacre using a linear model. Recently, it been shown experimentally that nacre also exhibits a nonlinear mechanical response NanoLett2004Nacre ; Barthelat2007JMPS ; Kakisawa2008Nacre , as spider silk does.

Here we construct a nonlinear model of a nacre-like structure that generalizes the linear model studied in Ref. OkumuraPGG2001 and reflects the recent experimental results in Ref. NanoLett2004Nacre ; Barthelat2007JMPS ; Kakisawa2008Nacre . Using this model, we perform simulations and derive scaling laws for a sample with a line crack. The numerical simulations, and scaling laws derived from this model show that the stress concentration near the crack tips is reduced significantly in the nonlinear case when the stress concentration is significant. The scaling laws demonstrate that toughness and strength are enhanced by a common factor that elucidates simple design principles for developing artificial materials mimicking nacre.

II Results

II.1 Scaling analysis of our model for nacre

In our analytical model, hard and thick layers of thickness $d_{h}$ are glued together with soft and thin layers of thickness $d_{s}$ ($d_{h}>>d_{s}$), with a period

[TABLE]

as illustrated in Fig. 1(a) (The corresponding lattice model for numerical simulation is illustrated in (b), which will be explained below in Sec. II.4.1). In our nonlinear model, the relation between the characteristic magnitude of stress and strain ($\sigma$ and $e$) are given by

[TABLE]

Here, the nonlinear exponent ($n_{s}$ or $n_{h}$) is positive, and the elastic modulus is denoted $E_{h}$ for the hard layers and $E_{s}$ for the soft layers.

As seen below, it is convenient to introduce a small parameter $\varepsilon$ defined as

[TABLE]

In our analytical model, we consider the small $\varepsilon$ limit and we are not interested in the exact values but just the orders of magnitude of the physical quantities.

In the small $\varepsilon$ limit, we show below that the dominant components of the stress and strain fields satisfy the following relations:

[TABLE]

Here, note that the Young’s and shear moduli (e.g., of hard layers) are of the same order of magnitude ($\sim E_{h}$) because they only differ by a numerical coefficient on the order of one. The continuum strain field is defined as $e_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ for $(x_{1,}x_{2},x_{3})=(x,y,z)$. Note that such continuum fields are valid in the present layered structure only on scales larger than the period $d$. Equation (4) is a natural nonlinear extension of the linear relation obtained for the simple model of nacre discussed in Ref. OkumuraPGG2001 .

Now we derive Eq. (4) by considering the elemental modes of deformation appropriate for shear ($\sigma_{yx}$) and stretch ($\sigma_{yy}$) in the $y$ direction. When the composite is sheared in the $y$ direction, only the soft layers are stretched, and $e_{yx}\sim(d_{s}/d)e_{yx}^{(s)}$. From Eq. (2), this results in $\sigma_{yx}\sim\varepsilon E_{h}e_{yx}^{1/n_{s}}$. In contrast, when the composite is stretched in the $y$ direction, both the soft and hard layers are stretched so the stress is dominated by the response from the hard layers, and $\sigma_{yy}\sim E_{h}e_{yy}^{1/n_{h}}$. Here, we have not considered any slip or separation between the interfaces of the layers, because we are interested in deriving the critical condition for failure. The mutual sliding of hexagonal hard platelets observed experimentally in a study on the toughness of nacre NanoLett2004Nacre is regarded as a state that exists after the critical condition for failure is satisfied and is beyond the scope of our description.

We consider a composite governed by Eq. (4) with a crack running in the $x$ direction. The composite is stretched in the $y$ direction as shown in Fig. 1(a), where the plane strain or plane stress condition is satisfied (the sample is thick or thin in the $z$ direction). Thus, the dominant stress components are those pointing in the $y$ direction, or $\sigma_{yy}$ and $\sigma_{yx}$. This is the reason we have considered only these components in the above.

As explained in the Appendix, we can derive the following scaling laws for the fracture toughness and strength of the composite in the presence of a large crack when $\varepsilon$ is small:

[TABLE]

where the common enhancement factor $\lambda$ is given by

[TABLE]

Here, $G_{h}$ and $\sigma_{F}^{MN}$ are the fracture toughness (fracture surface energy) and strength of a brittle monolithic hard material, respectively, and $a_{0}$ is the size of the so-called Griffith cavity, as discussed below. These scaling laws can be viewed as natural nonlinear extensions of results obtained from an analytical solution in Ref. OkumuraPGG2001 (see Appendix).

As explained in the Appendix, from Eqs. (5) and (6), we can derive the following expression for the maximum stress that appears near a tip of a large crack in a network model of system size $L$ when the remote stress $\sigma_{0}$ is applied at the top and bottom edges of the sample:

[TABLE]

The physical implications of Eqs. (7) and (8) will be discussed later.

II.2 Relevance of our analytical model to previous experiments

In this section, we show the relevance of our analytical model developed above to biological nacre by showing that Eq. (4) with appropriate parameters well reproduces previous experiments. Experimental studies have shown the tensile elastic response of the hard layers in nacre is linear with a modulus, typically around 65 GPa (before microscopic internal failure), as confirmed by nano-indentation NanoLett2004Nacre and tension tests Barthelat2007JMPS . This means that Eq. (4) with $n_{h}=1$ and $E_{h}\simeq 65$ GPa reproduces experiments well. Additionally, the shear response of the wet composite has recently been obtained Barthelat2007JMPS , and an experimental stress-strain curve is shown by the star symbols in Fig. 2. This curve can be reproduced at a semi-quantitative level by Eq. (4) with $\varepsilon E_{h}\simeq 72$ MPa and $n_{s}=4$, as shown by the solid curve in Fig. 2. This means that, in biological nacre, the soft layer is weakly nonlinear with $n_{s}\simeq 4$, and $\varepsilon$ is significantly small with $\varepsilon$ $\simeq 0.001$ ($\varepsilon E_{h}\simeq 72$ MPa and $E_{h}\sim 65$ GPa result in $\varepsilon\sim 0.001$).

The estimates, $E_{h}\sim 65$ GPa, $\varepsilon E_{h}\simeq 72$ MPa, and $n_{s}\simeq 4$, imply that the soft layer is very soft with $E_{s}\simeq 35$ MPa. This is because the overall volume of the soft layer is typically 5 per cent, implying $d_{s}/(d_{s}+d_{h})\simeq 1/20$ with $d_{h}$ typically 0.5 $\mu$m, and the parameter set $E_{s}\simeq 35$ MPa, $E_{h}\sim 65$ GPa, $d_{s}/d_{h}=1/19$, and $n_{s}=4$ in Eq. (3) gives $\varepsilon E_{h}\sim 72$ MPa. Note that, although the measurement of the force response of the soft component is difficult because of extremely small sample sizes of the soft component, it has been known that the soft component behaves like a soft gel Smith1999NatureNacre ; Kakisawa2008Nacre ; Meyers2008NacreSoft .

II.3 Physical implications of the scaling laws

Here, we discuss the rough orders of magnitude of important parameters in Eq. (7), to better understand its physical implications. The characteristic size of cavities in the hard layers, or the Griffith cavity $a_{0}$, is about $50$ nm, comparable to the thickness of the soft layers $d_{s}$, or $a_{0}\cong d_{s}$. The cavities or defects in the hard layers may be created from protein molecules in the soft layers intruding into, and being trapped within, the mineral crystals during biomineralization process in nacre Gao2003PNAS . The fracture toughness $G_{h}$ is on the order of 10 J/m2, which leads to $G_{h}/E_{h}\sim 0.1$ nm. In summary, typical values of the important factors in the common enhancement factor $\lambda$ can be given as

[TABLE]

Note that this parameter set means $\varepsilon E_{h}\sim 72$ MPa at $n_{s}=4$.

In Fig. 3(a), we show Eq. (7) with the parameter set in Eq. (9). The enhancement factor $\lambda$ as a function of $n_{s}$ sharply rises up until $n_{s}\sim 5-10$ and then tends to plateau. This implies that a relatively weak nonlinearity is sufficient for toughening. Note that in Fig. 3(a), the order of the toughness enhancement $\lambda$ is predicted in the thousands, which is consistent with previous experiments JacksonVincentTurner .

The same physical conclusion can be obtained from Eq. (8) when the ratio $L/d$ and the remote strain $e_{0}$ are relatively large, where $e_{0}=\sigma_{0}/E_{h}$ in the present case of $n_{h}=1$. In Fig. 3(b), Eq. (8) is plotted at $E_{s}/E_{h}=35/65000$ and $d_{s}/d_{h}=1/19$ for different values of $L/d$ and $e_{0}$. Equation (8) represents the maximum stress at the crack tip, and is a measure of strength because a material is considered to be strong when the value of this quantity is small. As seen in Fig. 3(b), this quantity drops sharply up until $n_{s}\sim 5$ and then it approaches a plateau when $e_{0}$ and $L/d$ are relatively large, which suggests that weak nonlinearity is sufficient for toughening. (The sharp drop in maximum stress is less pronounced when $e_{0}$ and $L/d$ are smaller, and this behavior is physically interpreted in the Discussion in terms of the degree of stress concentration.)

II.4 Numerical simulation

II.4.1 Lattice model

We performed numerical simulation using a nonlinear model extended from a linear two-dimensional network model of nacre studied in Ref. Aoyanagi2009PRE . The lattice model for this simulation is illustrated in Fig. 1(b). We consider a grid system in which nodal points are connected by springs that reflect the nonlinearity specified in Eq. (2). The values of $E_{s}$ and $E_{h}$ in our simulation are set as $E_{h}=65$ GPa and $E_{h}=35$ MPa, with $d_{h}/d_{s}=19$, as in Eq. (9), and with $n_{h}=1$ and $n_{s}=4$ to mimic nacre. (See Appendix for further details).

To study the strength of the system against failure, we introduce a line crack in the $x$ direction in a stretched system, and quantify how the stress is concentrated around the crack tips for a fixed deformation at the top and bottom edges of the sample at mechanical equilibrium. The two crack tips are placed at soft layers. This is because the critical condition for failure on a scale larger than the layer period $d$ should be the condition whether a crack develops further when the crack tip is stopped at a soft layer. We calculate the positions of the beads when the elastic energy is minimized by the conjugate gradient method NumericalRecipes .

For a systematic comparison, we consider three models, (1) a nonlinear layered model which mimics natural nacre with $\varepsilon\ll 1$ and $(n_{h},n_{s})=(1,4)$, (2) a linear layered model, which is different from the natural nonlinear model only in that $n_{s}$ is set to one, and (3) a monolithic model ($\varepsilon=1$), which is different from the linear model only in that the value of $E_{s}$ is raised to match that of $E_{h}$. We set the remote strain $e_{0}$ to be $0.008$ in each simulation, which is one of the values of $e_{0}$ used in Fig. 3(b).

II.4.2 Stress distribution in the presence of a crack

In Fig. 4(left), a comparison of the stress distribution around a line crack is given for the three models with a crack of the same size, where the maximum of the color scale (red) is set to the maximum stress that appears at the crack tips in the linear model. By comparing the three distribution maps on the left, we can clearly see the stress concentration near the crack tips is significantly reduced and delocalized in the layered cases (a and b), and reduced and delocalized most in the nonlinear case in (a).

To quantify the stress near the crack tips, we plot how the stress changes around one of the crack tips along the $x$ axis in Fig. 4(right) to more clearly see how the stress concentration is minimized in the nonlinear case. As expected, in all the three models, the maximum stress appears at the crack tip, or more precisely, at $r/d_{0}=1$, where $r$ is the distance from the crack tip and $d_{0}$ is the mesh size of the network model. In the nonlinear model, the enhancement factor for the crack tip stress compared with the remote value $\sigma_{0}$ is rather small, and is approximately $2$. The enhancement factors for the linear model and the monolithic model are approximately $4$ and $14$, respectively.

The enhancement of the stress at the crack tip significantly increases as we move from nonlinear to linear to monolithic models. This is in accord with the physical implications demonstrated in Fig. 3, that weak nonlinearity in the soft layer enhances the strength of the layered structure.

II.5 Consistency of our scaling analysis and numerical simulation

In this section, we demonstrate that the results of the numerical simulation semi-quantitatively agree with the scaling analysis. Note that the lattice model for numerical simulation should reduce to the analytical model in the continuum limit.

Equation (8) predicts the correct order of the maximum stress at the crack tip in the three numerical models in Fig. 4(right). The ratios of the maximum stress for nonlinear, linear composite and monolithic models are calculated as $1:2:6$, respectively, from Eq. (8), while the corresponding ratios are given as $1:2:7$ from Fig. 4(right). This level of agreement is remarkable, considering that the assumptions made in the derivation of the scaling laws are only reasonably well satisfied in our simulation for numerical limitations (see Appendix for the details). In this way, we conclude that Eq. (8), and thus Eqs. (5) and (6), are consistent with and supported by the results of our numerical calculation in Fig. 4.

To further confirm semi-quantitative agreement between our simulation and scaling analysis, we plot the maximum stress that appears at a crack tip as a function of the nonlinear index $n_{s}$ in Fig. 5. The three curves in Fig. 5 may be compared with the three solid curves in Fig. 3(b), because they were obtained at the same strain $e_{0}=0.008$ and for the same three values of $L/d$. As expected, the curves in each figure semi-qualitatively agree, further supporting the scaling analysis.

III Discussion

A crucial factor for toughening is stress concentration, and the degree of stress concentration is governed clearly by the remote stress $e_{0}=\sigma_{0}/E_{h}$ and less clearly by the quantity $L/d$. An increase in $L/d$ results in an increase in stress concentration. This is because the stress at a distance $r$ from the crack tip is given as $\sigma(r)\simeq\sigma_{0}(L/r)^{1/(1+n)}$ for a non-linear system governed by the stress-strain relation $\sigma\simeq Ee^{1/n}$, and the maximum stress that appears at the crack tip is given by a cutoff value $\sigma_{0}(L/d)^{1/(1+n)}$ for a system with a scale $d$ below which the continuum description breaks down aoyanagi2017 .

The fact that both $\sigma_{0}/E_{h}$ and $L/d$ control the stress concentration is reflected in Eq. (8). Analytically, the functional dependences of Eq. (8) on $e_{0}$ and $L/d$ are very similar. In fact, as demonstrated numerically in Fig. 3(b), the initial drop in the maximum stress as a function of $n_{s}$ becomes significant as $L/d$ increases, and the same is true as $e_{0}$ increases.

This change with $e_{0}$ of the initial drop in the maximum stress with $n_{s}$ is physically understood from the fact that, if we simply plot $\sigma=Ee$ and $\sigma=Ee^{1/n}$ with $n$ larger than unity, the nonlinear stress ($Ee^{1/n}$) is smaller than the linear stress ($Ee$) when the strain $e$ is larger than unity (the nonlinear stress is larger when $e$ is smaller than unity). This suggests that nonlinearity in the soft layer contributes to stress reduction effectively when strain acting on soft layers near the crack tip is large. This happens even if the remote strain acting on the composite is considerably small. Note that the strain in soft layers is significantly larger than the strain in the soft-hard composite near a crack tip where shear deformation is dominant. Even if the remote strain acting on the composite is small, shear strain of soft layers near a crack tip can be large, as suggested in the paragraph below Eq. (4). (This significantly enhanced deformation in soft layers near a crack tip is demonstrated in previous numerical calculations Aoyanagi2009PRE ; Hamamoto2013 .) In summary, when the stress concentration is high, its reduction due to nonlinearity becomes greater.

The initial drop in the maximum stress as a function of $n_{s}$ is one advantage of weak nonlinearity, and this advantage is more effective as the stress concentration becomes higher, or as $e_{0}$ and/or $L/d$ increases. This behavior of the maximum stress was demonstrated in Fig. 3(b), and can physically understood as explained above.

Another advantage of weak nonlinearity is the increase in the enhancement of $\lambda$ with an increase in $n_{s}$. This is predicted from Eq. (7), as demonstrated in Fig. 3(a), and this expression is obtained in the limit of large $L/d$, or in the case where the stress concentration is significant. This implies that if the stress concentration at a crack tip is low, this second advantage of weak nonlinearity predicted in Eq. (7) may be diminished. This is expected because stress-concentration reduction due to nonlinearity in soft layers (the first advantage) is more efficient when the stress concentration at a crack tip is high, and stress-concentration reduction and the enhancement factor $\lambda$ should be positively correlated. In summary, the two advantages of weak nonlinearity, originating in Eqs. (7) and (8), are physically the same, and the strength of the advantages enhance as $e_{0}$ and $L/d$ increase, or stress concentration increases.

In the case of biological nacre, it has been reported that hard sheets tend to be separated into hexagonal tablets at a significantly low strain of around 0.001 Barthelat2007JMPS . If this is always the case, then the two advantages of weak nonlinearity may not be so significant in real, biological nacre. This is because as seen in Fig. 3(b) the first advantage of weak nonlinearity is visible at $e_{0}=0.008$, and the two advantages tend to decrease as $e_{0}$ decreases (the curve for $e_{0}=0.002$ in Fig. 3(b) shows only a slight initial drop). Note that toughness and strength are governed by the strength of the two advantages at the critical state of failure when $e_{0}$ reaches the threshold strain for separation of tablets.

However, the value of around 0.001 reported in Barthelat2007JMPS is just one sample of nacre, and other natural specimens may vary. If other samples of nacre possess threshold strains which are fairly larger than 0.001, and their mechanical responses are explained by weak nonlinearity in the soft layers, then the weak nonlinearity could be interpreted as a mechanical optimization attained by biological nacre. This is because, if the threshold strain is fairly large, at the critical state of failure, the initial drop in the maximum stress at the crack tip as a function of $n_{s}$ is rather steep, and the maximum stress approaches to a plateau at a small $n_{s}$, and thus weak nonlinearity is sufficient to obtain a high strength of the two advantages.

At any rate, since the present study demonstrates that the advantages of weak nonlinearity in the soft layer is pronounced for a reasonably small strain of 0.008, this mechanism could be useful for developing artificially tough layered structures mimicking nacre, where hard plates do not break into tablets. Note that reducing the strain below 0.008 in numerical simulation is technically difficult due to precision in calculations.

IV Conclusion

We demonstrated that in natural nacre the hard layers are linear and the soft layers are weakly nonlinear ($n_{s}\sim 4$) on the basis of recent experimental results, and studied mechanical properties on the basis of nonlinear models reflecting this nonlinearity both analytically and numerically. The combination of numerical simulation and scaling arguments led to the following conclusions for a reasonably small strain: (1) a nonlinear model significantly contributes to reduce the stress concentration at a crack tip, compared with a linear model, and (2) a weak nonlinearity, as selected in biological nacre, is sufficient for this reduction.

We analytically derived the enhancement factor $\lambda$ common to both the fracture strength and toughness. This factor, given in Eq. (7), shows that the mechanical superiority of nacre-like layered structures is controlled by the factors $d/a_{0}\sim d_{h}/d_{s}$ and $\varepsilon^{-1}=(E_{h}/E_{s})(d_{s}/d_{h})^{1/n_{s}}$, which are larger than one. These factors elucidate two design principles for the enhancement: (1) the soft-hard combination ($E_{s}<E_{h}$) and (2) the thin-thick combination with weak soft nonlinearity ($d_{s}<d_{h}$ and $n_{s}\gtrsim 1$). These principles imply that the soft-hard combination, hierarchical structure, and nonlinearity are all important for developing tough structures mimicking nacre.

Acknowledgement

This research was partly supported by Grant-in-Aid for Scientific Research (A) (No. 24244066) of JSPS, Japan, and by ImPACT Program of Council for Science, Technology and Innovation (Cabinet Office, Government of Japan). We thank Professor Edward Foley (Ochanomizu University) for editing the manuscript.

Appendix

IV.1 Simulation

The simulation model is composed of $960\times 10000$ grid points ($L=10000$ in the unit $d_{0}$), and the system size in the non-stretched state is 960$d_{0}\times 10000d_{0}$, with $d_{s}=d_{0}$ and $d=20d_{0}$ ($d_{h}=19d_{0}$). We introduced a line crack of size $a=320d_{0}$ (whose width is $d_{0}$) with the tips located at soft layers. We stretched the system in the $y$ direction perpendicular to the line crack. The strain in our simulation is set to 0.008, which corresponds to the experimental value of tensile strain on nacre at failure Barthelat2007JMPS , with $E_{h}/E_{s}=65000/35$ in the nonlinear and linear layered models, and $E_{h}/E_{s}=1$ in the monolithic model. Note that assumptions in the scaling theory, e.g., $L,a\gg d\gg d_{0}$, is relaxed (or only marginally satisfied) in the simulation due to practical limitations on numerical calculation.

IV.2 Analytical model

The scaling structure of the elastic energy (per volume) can be given as

[TABLE]

where the numerical coefficients of the both terms are set to one for simplicity. Equation (10) is obtained from the expression $F=(\sigma_{yy}e_{yy}+\sigma_{yx}e_{yx})/(n+1)$, together with Eq. (4).

Equation (10) is a natural nonlinear extension of the energy, first derived in Ref. OkumuraPGG2001 . On the basis of this expression, we develop scaling arguments to derive Eqs. (5), (6), and (8) in the text. These expressions can be viewed as a natural nonlinear extension of those obtained in Ref. OkumuraPGG2001 , and by setting $n_{s}=n_{h}=1$ in Eqs. (5), (6), and (8), we can reproduce the corresponding expressions obtained in Ref. OkumuraPGG2001 from an exact analytical solution in appropriate limits. The details of the derivation of the scaling laws are given below.

IV.2.1 Details of Derivation

To obtain the scaling laws for the nonlinear model of nacre, we consider a crack problem in a generic way. A crack of size $a$ is located along the $x$ axis in a plate whose size in the $y$ direction is characterized by $L$ (the sample size in the $x$ direction is comparable to $L$). This plate is in tension. The following derivation is inspired by the arguments in Ref. okumura2005fracture .

We consider two characteristic length scales, $X$ and $Y$, in the $x$ and $y$ directions. A scaling relation between $X$ and $Y$ is given by balancing the two terms in the energy in Eq. (10):

[TABLE]

The critical condition for fracture is that the elastic energy that is released as a result of the creation of a crack of size $a$ matches the fracture energy to create the crack:

[TABLE]

Here, the fracture energy per area is denoted $G$. Note that if Eq. (11) holds, the two terms in Eq. (10) are of the same order.

From Eqs. (11) and (12), we can express $u$ and $Y$ in terms of $X$, and from these expressions, we obtain the following expression for the $(y,y)$ component of the stress at the critical state of failure that scales as $\sigma\sim E_{h}(u/Y)^{1/n_{h}}\sim G/u$:

[TABLE]

This is a generalized version of a nonlinear Griffith’s failure formula that was proposed theoretically in Ref. Aoyanagi2009JPSJ and confirmed experimentally in Ref.SoneMoriJSPS . If we set $\varepsilon=1$ and $n_{s}=n_{h}$ in Eq. (13) we obtain a nonlinear Griffith’s formula for a monolithic hard material:

[TABLE]

Here, we have set $G=G_{h}$ and $X=X_{h}$ (in fact, $X_{h}=Y_{h}$ because the system is isotropic), where $G_{h}$ and $X_{h}$ ($Y_{h}$) stand for the fracture energy and characteristic scale in the $x$ ($y$) direction for the monolithic hard material.

The nonlinear Griffith formula in Eq. (14) is a nonlinear version of the classic Griffith formula, and in the linear case ($n_{h}=1$), Eq. (14) reduces to the well-known expressions for failure stress: $\sigma\sim\sqrt{E_{h}G/a}$ for $a<L$ and $\sigma\sim\sqrt{E_{h}G/L}$ for $a>L$.

If we set $X_{h}=a_{0}$ in Eq. (14), we obtain the intrinsic failure stress of hard layers when there are no macroscopic cracks:

[TABLE]

Here, $a_{0}$ is the typical size of defects in the layers that play the role of Griffith cavities.

The scaling structure for the singularity of the stress field at a crack tip should be given in the form,

[TABLE]

where $\sigma_{0}$ and $r$ are the characteristic size of the remote stress and the distance from the crack tip, respectively. The exponent $\eta$ can be determined by a principle that, at the critical state of failure [$\sigma_{0}$ is equal to $\sigma_{F}$ in Eq. (13)], the field $\sigma(r)$ becomes independent of $X$ as we approach the singularity at $r=0$ ($r<<X$) Okumura2004DNG . From this principle, we obtain the desired exponent,

[TABLE]

and the singular field at the critical state of failure,

[TABLE]

Note that this principle allows us to reproduce the well-known crack-tip singularity obtained by the Hutchinson, Rice and Rosengren (HRR) Hutchinson1968 ; RiceRosengren1968 , and if we set $n_{s}=n_{h}=n$, $E_{s}=E_{h}$, and $d_{s}=d_{h}$, Eq. (16) reduces to the well-known expression of the HRR singularity: $\sigma(r)\sim\sigma_{0}(X/r)^{\frac{1}{1+n}}$. This further reduces to the well-known singularity in the linear case ($n=1$): $\sigma(r)\sim\sigma_{0}\sqrt{X/r}$.

The maximum force allowed at the crack tip is given by Eq. (16) at $r\sim d$. This is expected because the continuum expression should be cut-off at this scale, and the continuum theory is valid only beyond length scales larger than the layer period $d$. From this second principle, the maximum stress at the crack tip at the critical state of failure is given by Eq. (18) with $r=d$.

In fact, this second principle is numerically confirmed in Refs. Nakagawa ; Aoyanagi2009JPSJ ; Takahashi2014 . The authors considered a line crack of size $a$ in a two dimensional nonlinear network system whose mesh size is $d$, and showed that, as expected, the maximum stress always appears at the crack tips. Moreover, they showed that the maximum stress follows the scaling law $\sigma\simeq\sigma_{0}\left(a/d\right)^{1/(1+n)}$ in a clear way for the remote stress $\sigma_{0}$, whereas the well-known singularity $\sigma\simeq\sigma_{0}\left(a/r\right)^{1/(1+n)}$ near a crack tip at the distance $r$ is predicted in Refs. Hutchinson1968 ; RiceRosengren1968 .

The critical condition for failure is satisfied when this maximum stress matches $\sigma_{h}$ in Eq. (15), and at the critical state, the tip stress reaches the intrinsic failure stress of hard layers. This matching condition leads to an expression of the fracture energy $G$:

[TABLE]

The enhancement factor for the fracture energy is given in the left-hand side.

Another scaling relation between $X$ and $Y$ is obtained by considering the characteristic size of the remote stress $\sigma_{0}$ and the remote strain $e_{0}$. They satisfy the relationships

[TABLE]

From Eqs. (11) and (20), we obtain a useful relation between $X$ and $Y$:

[TABLE]

This implies that $(X,Y)=(\tilde{\varepsilon}L,L)$ when $a/\tilde{\varepsilon}>L$ and $(X,Y)=(a,a/\tilde{\varepsilon})$ otherwise. This is because the natural characteristic length $X$ is $a$ if $L>>a$, but, if $L$ is not so larger than $a$ and if the length scale $Y\sim a/\tilde{\varepsilon}$ is cutoff by $L$ (i.e., if $a/\tilde{\varepsilon}\gtrsim L$), then $Y$ should scale as $L$, which means $X\sim\tilde{\varepsilon}L$. In the case of the monolithic hard material, $X_{h}\sim Y_{h}$ so that $(X_{h},Y_{h})=(L,L)$ when $a\gtrsim L$ and $(X_{h},Y_{h})=(a,a)$ otherwise.

The enhancement factor for the fracture stress is obtained as

[TABLE]

This can be obtained in the following three ways, which are physically the same: (I) Substitute Eq. (19) into Eq. (13). (II) Match Eq. (16) evaluated at $r=d$ with $\sigma_{h}$ and interpret the remote stress $\sigma_{0}$ in Eq. (16) as the fracture strength $\sigma_{F}$. (III) Calculate $\sigma_{0}u_{0}$ where $u_{0}$ is given by $u$ in Eq. (20) and identify $\sigma_{0}$ as $\sigma_{F}$.

When the crack is large so that $X_{h}=Y=L$, we obtain Eq. (6) from Eq. (22). Equation (8) is obtained for such large cracks by setting $r=d$ and $X=\tilde{\varepsilon}L$ in Eq. (16), with using Eq. (20). While Eqs. (5) and (6) are valid only at the critical state of failure, Eq. (8) is valid as long as $\sigma_{0}$ does not exceed the critical stress $\sigma_{F}$.

When the crack is small so that $X_{h}=a$ and $Y=a/\tilde{\varepsilon}$, we obtain from Eq. (22):

[TABLE]

For such small cracks, by setting $r=d$ and $X=a$ in Eq. (16), with using Eq. (20).

[TABLE]

The large crack assumption may be satisfied in our numerical calculations. Note that for $\varepsilon\simeq 0.001$ and $e_{0}=0.008$, then $\tilde{\varepsilon}\simeq 0.1$ so that $a>k\tilde{\varepsilon}L$ holds if k\,\is slightly smaller $1/3$ (i.e., $k$ is of the order of unity) even if $L/d_{0}=10000$, which is the largest value we used in numerical calculation because $a/d_{0}=320$. (If otherwise, $\sigma_{M}$ would be $L/d$-independent.)