A New Biomimetic Composite Structure with Tunable Stiffness and Superior Toughness via Designed Structure Breakage

Mimicking natural structures has been highly pursued recently in composite structure design to break the bottlenecks in the mechanical properties of the traditional structures. Bone has a remarkable comprehensive performance of strength, stiffness and toughness, due to the intricate hierarchical microstructures and the sacrificial bonds within the organic components. Inspired by the strengthening and toughening mechanisms of cortical bone, a new biomimetic composite structure, with a designed progressive breakable internal construction mimicking the sacrificial bond, is proposed in this paper. Combining the bio-composite staggered plate structure with the sacrificial bond-mimicking construction, our new structure can realize tunable stiffness and superior toughness. We established the constitutive model of the representative unit cell of our new structure, and investigated its mechanical properties through theoretical analysis, as well as finite element modeling (FEM) and simulation. Two theoretical relations, respectively describing the elastic modulus decline ratio and the unit cell toughness promotion, are derived as functions of the geometrical parameters and the material parameters, and validated by simulation. We hope that this work can lay the foundation for the stiffness tunable and high toughness biomimetic composite structure design, and provide new ideas for the development of sacrificial bond-mimicking strategies in bio-inspired composite structures.


Introduction
Biological materials, such as bone in mammals and nacre in shells, have attracted constant research interest for their excellent combination of strength, stiffness and toughness [1][2][3][4][5]. For example, bone is a composite material made up of collagen (30%-45% by volume), apatite crystals (30%-50% by volume) and small amounts of non-collagenous proteins [6,7]. The Young's modulus is about 10 GPa, and the tensile strength is about 80-120 MPa [6][7][8][9]. The key reason for its excellent mechanical performance owes to the "staggered lamellae layer", an intricate hierarchical micro-structure inside the cortical bone [10][11][12][13][14]. Besides this, it has also been revealed that the sacrificial bonds within the organic components are among the important factors which account for its excellent energy dissipating mechanisms [15][16][17][18][19]. Sacrificial bonds are defined as relatively weak bonds (and often reversible) that rupture before strong bonds fail under deformation [20]. A load applied to the biological materials would lead to the rupture of sacrificial bonds, which means a huge amount of energy dissipation and a promotion of the toughness of the biological. Meanwhile, due to the retention of strong bonds, the strength of the biological materials would rarely be influenced.
With the established constitutive model of the unit cell of our new "biomimetic staggered structure", we studied the relationship between its mechanical properties and its structural parameters. Although we have not built a laboratory sample of the new "biomimetic staggered structure" for experimental testing, an FEM simulation was carried out to verify the effectiveness of the designed structure. The FEM simulation results agreed well with our theoretical analysis. Ultimately, we wish to put forward a new design and optimization idea for meso/macro-scale biological structures, mimicking "sacrificial bonds" to achieve superior toughness and tunable stiffness mechanical properties.

Mechanical Model and Theoretical Analysis
Considering the periodicity and symmetry of our structure, we take a unit cell, as defined in Figure  1b, to establish the constitutive model, and further analyze the trends in structure parameters and mechanical response. We defined two phases for a unit cell, "phase 0": before the "joint part" breakage; and "phase 1" after the "joint part" breakage. With the designed progressive breakage happening under load, the unit cells would change from "phase 0" to "phase 1" one by one. These will reflect the structure parameters (such as stiffness and toughness) of the whole biomimetic structure changing gradually.

Unit Cell at "Phase 0"
The geometry and material parameters of a unit cell at "phase 0" are shown in Figure 2a. Stress definitions and distribution in a deformed unit cell under uniaxial tensile load, with region partition, are shown in Figure 2b. In a unit cell at "phase 0", there are "hard plate" (in Regions 1,2 and 4), "shear part" (in Region 3), and "joint part" (in Region 5). We make the simplification that "hard plate" and "joint part" can only hold normal stress, while "shear part" only hold shear stress under uniaxial tensile load. We assume that all parts maintain the characteristics of linear elastic deformation. Moreover, under our design the "joint part" is the only breakable region, and the fracture form is brittle fracture. Therefore, the "maximum tensile stress principle" was taken as the failure criterion for the "joint part".

Mechanical Model and Theoretical Analysis
Considering the periodicity and symmetry of our structure, we take a unit cell, as defined in Figure 1b, to establish the constitutive model, and further analyze the trends in structure parameters and mechanical response. We defined two phases for a unit cell, "phase 0": before the "joint part" breakage; and "phase 1" after the "joint part" breakage. With the designed progressive breakage happening under load, the unit cells would change from "phase 0" to "phase 1" one by one. These will reflect the structure parameters (such as stiffness and toughness) of the whole biomimetic structure changing gradually. The geometry and material parameters of a unit cell at "phase 0" are shown in Figure 2a. Stress definitions and distribution in a deformed unit cell under uniaxial tensile load, with region partition, are shown in Figure 2b. In a unit cell at "phase 0", there are "hard plate" (in Regions 1,2 and 4), "shear part" (in Region 3), and "joint part" (in Region 5). We make the simplification that "hard plate" and "joint part" can only hold normal stress, while "shear part" only hold shear stress under uniaxial tensile load. We assume that all parts maintain the characteristics of linear elastic deformation. Moreover, under our design the "joint part" is the only breakable region, and the fracture form is brittle fracture. Therefore, the "maximum tensile stress principle" was taken as the failure criterion for the "joint part".
As defined in Figure 2a, dimension b is the half thickness of "hard plate" and "joint part", h is the thickness of the shear part; l a is the length of the shear part (also the uniformly overlapped length of the hard plate); 2l a is the length of the "joint part" (also the non-overlapped length of the hard plate). In this study, a uniformly overlapped structure assumption is due to the more efficient load transfer capability [30,31]. The tensile modulus for "hard plate" and "joint part" are denoted as E m and E e ; while the shear modulus for "region 3" are denoted as G. To easily describine the mechanical model, we also defined several non-dimensional geometrical and material parameters. The geometrical parameters are: the "hard plate" overlapped length to thickness ratio ρ = l a /b; the "Region 3" length to thickness ratio λ = l a /h; the "hard plate" overlapped length to non-overlapped length ratio η = l a /l b ; and the approximate volume fraction of "hard plate" φ = 2b/(2b + h). The material parameters are: the "joint part" tensile modulus to "hard plate" tensile modulus α = E e /E m ; and the "Region 3" shear modulus to "hard plate" tensile modulus β = G/E m .  As defined in Figure 2a, dimension b is the half thickness of "hard plate" and "joint part", h is the thickness of the shear part; a l is the length of the shear part (also the uniformly overlapped length of the hard plate); 2 a l is the length of the "joint part" (also the non-overlapped length of the hard plate). In this study, a uniformly overlapped structure assumption is due to the more efficient load transfer capability [30,31]. The tensile modulus for "hard plate" and "joint part" are denoted as m E and e E ; while the shear modulus for "region 3" are denoted as G . To easily describine the mechanical model, we also defined several non-dimensional geometrical and material parameters.
The geometrical parameters are: the "hard plate" overlapped length to thickness ratio ; and the approximate volume fraction of "hard plate" ( ) At "phase 0", under the coordinate defined in Figure 2b, we established the constitutive model of the unit cell, applying the well-known 'shear-lag' model [26].
where the subscripts stand for the region number.
When considering the boundary conditions, because of the uniformly overlapped staggered configuration, we assumed that stress in Region 4 ( 4 σ ) and Region 5 ( 5 σ ) were both uniform. Under this assumption, a boundary condition can be written as  The region partition 1 to 5 is defined on the deformed unit cell, where Regions 1, 2 and 4 are "hard plate" material, region 3 is "share part" material and Region 5 is "joint part" material. σ i (i = 1, 2, 4, 5) and τ are respectively, the normal and shear stress distributed in every region, under uniaxial tensile stress. Origin "o" and "x-axis" is the coordinate defined while establishing the constitutive model.
At "phase 0", under the coordinate defined in Figure 2b, we established the constitutive model of the unit cell, applying the well-known 'shear-lag' model [26].
where the subscripts stand for the region number. When considering the boundary conditions, because of the uniformly overlapped staggered configuration, we assumed that stress in Region 4 (σ 4 ) and Region 5 (σ 5 ) were both uniform. Under this assumption, a boundary condition can be written as The other boundary condition was obtained by considering the force equilibrium between Region 5 and the left end of Region 2 (or the right end of Region 1).
From the force equilibrium of the whole unit cell, σ 4 and σ 5 should correspond to the following equation φ where φ = 2b/(2b + h) is the approximate volume fraction of "hard plate", and σ is the volume-averaged stress in the unit cell under uniaxial tensile load. Solving Equation (1) with boundary conditions Equations (2) and (3), and simplified with Equation (4); we obtained the stress distribution in Regions 1,2 and 3, at "phase 0": where α = E e /E m , β = G/E m , ρ = l a /b, λ = l a /h, η = l a /l b are the non-dimensional geometrical or material parameters defined before; K is a non-dimensional parameter that represents the complex geometrical and material effects of Regions 1, 2 and 3.
2.1.2. Unit Cell at "Phase 1" Next, we established the changed constitutive model of a unit cell at "phase 1", also under the uniaxial tensile load. In the unit cell at "phase 1", the Region 5 was breaking under deformation and there was only "hard plate" (in Regions 1,2 and 4) and "shear part" (in Region 3), as shown in Figure 3. Moreover, all the geometrical and material parameters of the unit cell can be inherited from "phase 0". Next, we established the changed constitutive model of a unit cell at "phase 1", also under the uniaxial tensile load. In the unit cell at "phase 1", the Region 5 was breaking under deformation and there was only "hard plate" (in Regions 1,2 and 4) and "shear part" (in Region 3), as shown in Figure  3. Moreover, all the geometrical and material parameters of the unit cell can be inherited from "phase 0". For unity of the analysis, we used the "shear-lag" model in Regions 1, 2, 3 and then determined the uniform 4 σ . Therefore, under the coordinate defined in Figure 3, the shear-lag model in Regions 1, 2 and 3 should still be written as Equation (1).
However, the boundary conditions would change as the following: Moreover, from the force equilibrium of the whole unit cell, 4 σ should be derived as Solving Equation (1) with boundary conditions Equation (8) and Equation (9), and simplified with Equation (10), we obtained the stress distribution in Regions 1,2 and 3 at "phase 1". For unity of the analysis, we used the "shear-lag" model in Regions 1, 2, 3 and then determined the uniform σ 4 . Therefore, under the coordinate defined in Figure 3, the shear-lag model in Regions 1, 2 and 3 should still be written as Equation (1).
However, the boundary conditions would change as the following: Moreover, from the force equilibrium of the whole unit cell, σ 4 should be derived as Solving Equation (1) with boundary conditions Equations (8) and (9), and simplified with Equation (10), we obtained the stress distribution in Regions 1, 2 and 3 at "phase 1".

Unit Cell Elastic Modulus Analysis
Taking the unit cell as the basic building block of our biomimetic structure, then the structure stiffness can be defined by the material parameters and the geometrical parameters of the unit cells. Therefore, it is necessary to study the elastic modulus of the unit cell.
Following the definition of Hill [32], the volume-averaged stress σ = 1 V V σdV and strain ε = 1 V V dV of a unit cell at both "phase 0" and "phase 1" were calculated separately, where V is the volume of a unit cell. Then, we obtained the effective elastic modulus (E) of the unit cell, at "phase 0" and "phase 1", denoted as E 0 and E 1 , in terms of non-dimensional parameters. At "phase 0" where A = Ktanh(K) is a non-dimensional parameter for written simplicity. At "phase 1" This formula indicates that the effective modulus of the unit cell must decline when the unit cell changes from "phase 0" to "phase 1", and the decline ratio is written as below.
Considering the convenience of the following simulation process and the improvement of the structure, the non-dimensional materials parameters α and β were referred to resin and rubber-like materials, which were widely used in the polyject method 3D-printing. On the other hand, the structural geometric configuration parameters are calculated by the theoretical model after the material parameters are determined, which can clearly reflect the changing trends of the equivalent elastic modulus and structural toughness of the designed structure. We chose the non-dimensional geometrical parameter η and the non-dimensional material parameter α as the independent variables, and studied their influence on E  Table 1. It is shown in Figure 4a that, as η increases, both E 0 and E 1 would increase with an attenuated increase rate, and finally approach their limiting value, derived by Equations (16) and (17) respectively. We particularly defined several critical values of η c upon differtent paramaters. First of all, we denoted the critical value of η for E 1 as η c . As η is the only influence factor on E 1 , we decided on η c as the when η ≥ η c . With our definition, η c can be derived by Equation (18).
However, as also shown in Figure 4a, E 0 would be influenced by not only η but also α; E 0 would increase as α increased from 0.01 to 0.1. Similarly, we denoted the critical value of η for E 0 as η c E 0 , derived as Equation (19), which is influenced by α.
As shown in Figure 4b, the effective modulus decline E  (21).
Usually, in order to realize a large tunable stiffness range in our biomimetic structure, a small  It is shown in Figure 4a that, as η increases, both 0 E and 1 E would increase with an attenuated increase rate, and finally approach their limiting value, derived by Equations (16) and (17) respectively. We particularly defined several critical values of c However, as also shown in Figure 4a, 0 E would be influenced by not only η but also α ; 0 E would increase as α increased from 0.01 to 0.1. Similarly, we denoted the critical value of η for 0 E as ( ) 0 E c η , derived as Equation (19), which is influenced by α .

Unit Cell Phase Changing Process Analysis
With the "joint part" breakage, the unit cell will change from "phase 0" to "phase 1" and the mechanical properties of the unit cell will change. Starting with the strength limit of the "joint part", we further studied the relationship between stress and strain during the unit cell phase changing process.
We defined the limiting tensile stress of the "joint part" as σ s e , the limiting tensile stress of the "hard plate" as σ s m , and the limiting shear stress of "Region 3" as τ s . In order to fulfill the assumption that the "joint part" is the only breakable part under our design, the strength parameters within the unit cell (σ s e , σ s m and τ s ) should satisfy the following relationships taking volume-averaged stress σ as the dependent variable and volume-averaged strain ε as the independent variable. Before the "joint part" breakage, the unit cell at "phase 0" should deform under the tensile load with the effective modulus E 0 . Then, we defined the critical average stress σ 0 cr : when σ ≤ σ 0 cr , the breakage will not happen and the unit cell at "phase 0". σ 0 cr could be derived by combining Equation (7) with σ s e : Further, we derived the critical average strain ε cr , which represents the volume-averaged strain ε of the unit cell when the phase-change is just about to happen. With the "generalized Hooke law" at "phase 0", ε cr was derived as the following Considering that, when the "joint part" breaking, the shape of the unit cell must not abruptly change, i.e., saltation is impossible for the strain field, we concluded that ε cr at "phase 0" is the same to that at "phase 1". Therefore, for the equilibrium of the unit cell, the saltation must happen within the stress field, which means the critical average stress at "phase 1" (σ 1 cr ) should not equal σ 0 cr , i.e., σ 1 cr σ 0 cr . Then, σ 1 cr can be derived along with the "generalized Hooke law" at "phase 1", as in the following: After changing into "phase 1", the unit cell will deform with the effective elastic modulus E 1 , until reaching the ultimate strength limit of the unit cell (decided by σ s m and τ s ). With all the characteristic parameters derived above, we obtained the volume-averaged stress-strain relationship of the whole "phase changing" process, as shown in Figure 5.

Structure Toughness Analysis
While the term "toughness" has multiple usages, the "structure toughness" mentioned in this paper is not the material toughness (usually defined as the maximum energy adsorbed per mass before fracture [19]). Instead, we explored hierarchical and structure parameter-dependent toughness-e.g., "structure toughness". This particular toughness is defined as the specific energy dissipation needed for the breakage of the "joint parts", which also reflected the extra toughness increase for our biomimetic structure during "phase changing", compared with the classical platestaggered structure.
We neither made assumptions of "joint part" reformation nor a reversible process. Therefore, the stress-strain response is unidirectional, and no hysteretic behavior needs considering. We denoted the normalized "structure toughness" of a unit cell as cell T , and defined it as the following: Noting that the "joint part" breaking event results in a repeated pattern within every unit cell, the "structure toughness" T for the whole biomimetic structure can be simply summarized as: Therefore, without loss of generality, we took the unit cell as before, and analyzed its cell T .
Substituting Equations (23 ~ 25) into Equation (26), then cell T can be derived as: Figure 5. The volume-averaged stress-strain curve for the unit cell during "phase changing", which is marked with the black line. The effective elastic moduli at "phase 0" and "phase 1" are marked with the red dashed line and blue dashed line, respectively. The critical average strain, ε cr , and the critical average stress, σ 0 cr , as well as σ 1 cr , is marked with black dashed line. The structure stiffness of the unit cell defined in this paper is encircled by the black line and the blue dashed line.

Structure Toughness Analysis
While the term "toughness" has multiple usages, the "structure toughness" mentioned in this paper is not the material toughness (usually defined as the maximum energy adsorbed per mass before fracture [19]). Instead, we explored hierarchical and structure parameter-dependent toughness-e.g., "structure toughness". This particular toughness is defined as the specific energy dissipation needed for the breakage of the "joint parts", which also reflected the extra toughness increase for our biomimetic structure during "phase changing", compared with the classical plate-staggered structure.
We neither made assumptions of "joint part" reformation nor a reversible process. Therefore, the stress-strain response is unidirectional, and no hysteretic behavior needs considering. We denoted the normalized "structure toughness" of a unit cell as T cell , and defined it as the following: Noting that the "joint part" breaking event results in a repeated pattern within every unit cell, the "structure toughness" T for the whole biomimetic structure can be simply summarized as: Therefore, without loss of generality, we took the unit cell as before, and analyzed its T cell . Substituting Equations (23)- (25) into Equation (26), then T cell can be derived as: As is shown in Figure 6 with η increasing, the normalized structure toughness of the unit cell T cell would decrease with an attenuated decrease rate, and finally approach the limiting value, written as Equation (29). As α increases from 0.01 to 0.1, T cell would decrease, and the critical value of η for T cell i.e., η c (T cell ), derived by Equation (30), is also influenced by α. (30) Moreover, to gain a high T cell for the unit cell, both α and η should be set as small as possible, under the compromise of E 1 /E 0 and the design constraints. Therefore, with certain material components (i.e., α is defined), the upper limit of η, i.e., η c (T cell ), for the unit cell will be determined by T cell requirement.
As is shown in Figure 6 with η increasing, the normalized structure toughness of the unit cell cell T would decrease with an attenuated decrease rate, and finally approach the limiting value, written as Equation (29). As α increases from 0.01 to 0.1, cell T    Figure 6. The structure toughness change with respct to the non-dimensional geometrical parameter η and the non-dimensional material parameter α, in the logarithmic coordinate. α changed from 0.01 to 0.1 discretely, shown with a black dashed line. η changed from 5 to 500 continuously, and the η c is marked with a black dash-dotted line. Other parameters used for making this figure were inherited from Table 1.

Finite Element Modeling and Simulation
We validated the generality and accuracy of our mechanical model by finite element analysis (FEA) with ABAQUS/CAE, steady-state static, direct solver. Then, we compared the simulation results and the theoretical results of the stiffness decline E 1 /E 0 , as well as the structure toughness T cell of the unit cell.

FEM Simulation Preparation
The cohesion FEM model can simulate the failure process of materials by introducing the cohesion model within the materials. In this paper, zero-thickness cohesion elements COH2D4 were inserted between every solid elements CPE4R element at Region 5 to simulate the designed progressive damage and breaking process of the "joint part" under external load. The zero thickness COH2D4 element can be regard as two connecting faces of adjacent solid elements. The cohesive elements describe the damage and failure, which have different forms, by the "traction separation law". To remain consistent with the theoretical analysis, we chose the "bilinear cohesive damage model" in the ABAQUS platform, as shown in Figure 7. Under external load, the cohesive element would firstly go through the elastic deformation stage, and when meeting the damage initiation criteria, the damage evolution stage would begin, and, after the cohesive element completely failed, it would be deleted.
progressive damage and breaking process of the "joint part" under external load. The zero thickness COH2D4 element can be regard as two connecting faces of adjacent solid elements. The cohesive elements describe the damage and failure, which have different forms, by the "traction separation law". To remain consistent with the theoretical analysis, we chose the "bilinear cohesive damage model" in the ABAQUS platform, as shown in Figure 7. Under external load, the cohesive element would firstly go through the elastic deformation stage, and when meeting the damage initiation criteria, the damage evolution stage would begin, and, after the cohesive element completely failed, it would be deleted.
In Figure 7, the abscissa and ordinate are, respectively, the cohesive element separation displacement and the traction stress; the slope of the elastic deformation stage is the cohesive element stiffness k ; i δ is the separation displacement at damage initiation; and s e σ is the traction stress at damage initiation, which is equivalent to the "maximum tensile stress" when the "Maxs" principle is selected as the cohesive element failure criterion in ABAQUS; f δ is the failure displacement; and the triangle area enclosed by the elastic deformation stage curve, the damage evolution stage curve and the abscissa is the fracture energy c G .  Figure 8. In order to keep consistent with theoretical analysis, isotropic linear elastic materials were chosen for all parts in the unit cell FEM model, and the "Maxs" principle was selected as the COH2D4 failure criterion. Moreover, the elastic modulus of the COH2D4 elements should be the same as the adjacent CPE4R elements. The geometrical and material properties used for the simulation are listed in Table 2 and  Table 3, respectively. Considering the periodic property of the structure, the periodic boundary In Figure 7, the abscissa and ordinate are, respectively, the cohesive element separation displacement and the traction stress; the slope of the elastic deformation stage is the cohesive element stiffness k; δ i is the separation displacement at damage initiation; and σ s e is the traction stress at damage initiation, which is equivalent to the "maximum tensile stress" when the "Maxs" principle is selected as the cohesive element failure criterion in ABAQUS; δ f is the failure displacement; and the triangle area enclosed by the elastic deformation stage curve, the damage evolution stage curve and the abscissa is the fracture energy G c .
The FEM model with mesh properties of the unit cell is shown in Figure 8. In order to keep consistent with theoretical analysis, isotropic linear elastic materials were chosen for all parts in the unit cell FEM model, and the "Maxs" principle was selected as the COH2D4 failure criterion. Moreover, the elastic modulus of the COH2D4 elements should be the same as the adjacent CPE4R elements. The geometrical and material properties used for the simulation are listed in Tables 2 and 3, respectively. Considering the periodic property of the structure, the periodic boundary conditions were set to this unit cell FEM model, and the uniform-speed displacement load was applied to the left and right edges, as shown in Figure 8.  conditions were set to this unit cell FEM model, and the uniform-speed displacement load was applied to the left and right edges, as shown in Figure 8.

Simulation Results and Discussion
Under the displacement load, the Region 5 break, as designed, and the unit cell FEM model, went through the phase changing process. Figure 9a-f displays the strain and stress contour figures of the deformed FEM model, changing from "phase 0" to "phase 1" during simulation. Especially, the gradient change of stress and strain along y direction owe to the stress concentration that occurred at the intersection of Region 3 and Region 1 (2), which is relatively small compared with the gradient change along x direction, and is neglected while theoretical modeling and analyzing. Figure 9b shows the stress distribution within the unit cell at "phase 0", under uniaxial tensile, which basically coincides with the trend theoretically derived by Equation (5). Figure 9f shows the stress distribution of the unit cell at "phase 1", which basically coincides with Equation (11).

Simulation Results and Discussion
Under the displacement load, the Region 5 break, as designed, and the unit cell FEM model, went through the phase changing process. Figure 9a-f displays the strain and stress contour figures of the deformed FEM model, changing from "phase 0" to "phase 1" during simulation. Especially, the gradient change of stress and strain along y direction owe to the stress concentration that occurred at the intersection of Region 3 and Region 1 (2), which is relatively small compared with the gradient change along x direction, and is neglected while theoretical modeling and analyzing. As shown in Figure 10, the stress-strain relationship during FEM simulation agrees well with the theoretical result. The fluctuation of the simulation curve from point a to point b happened during the phase changing process. This is because the breakage of Region 5 would start at point a (  Figure 9b shows the stress distribution within the unit cell at "phase 0", under uniaxial tensile, which basically coincides with the trend theoretically derived by Equation (5). Figure 9f shows the stress distribution of the unit cell at "phase 1", which basically coincides with Equation (11).
As shown in Figure 10, the stress-strain relationship during FEM simulation agrees well with the theoretical result. The fluctuation of the simulation curve from point a to point b happened during the phase changing process. This is because the breakage of Region 5 would start at point a (ε a cr = 1.031 × 10 −4 ) and gradually develop along y direction, as shown in Figure 9e Table 4.   Figure 11. It is predicted by the theoretical analysis and verified by this simulation that, when η is set, as α increases, 0 E would increase with an attenuated increase rate and finally approach a limiting value. As can be seen in the above analysis, the limited value should be defined by η .  Next, we simulated the influence of α on the unit cell effective modulus. In this simulation, α was set as 0.02, 0.05, 0.1, 0.14 and 0.2, by adjusting the material property E e successively, while η was kept as 100. The simulation result was compared with theoretical results, in which α continuously changed from 0.01 to 0.2, as shown in Figure 11. It is predicted by the theoretical analysis and verified by this simulation that, when η is set, as α increases, E 0 would increase with an attenuated increase rate and finally approach a limiting value. As can be seen in the above analysis, the limited value should be defined by η. e was kept as 100. The simulation result was compared with theoretical results, in which α continuously changed from 0.01 to 0.2, as shown in Figure 11. It is predicted by the theoretical analysis and verified by this simulation that, when η is set, as α increases, 0 E would increase with an attenuated increase rate and finally approach a limiting value. As can be seen in the above analysis, the limited value should be defined by η . Figure 11. Comparison between FEM simulation and theoretical results of the effective modulus changing with α , while 100 η = . In FEM simulation, represented by black circles, α was set as

Conclusions
In summary, we have brought up a new biomimetic composite structure with tunable stiffness and superior structure toughness via a designed progressive breakable constituent. We mainly focused on the periodic unit cell of the structure, established the mechanical model of the unit cell and verified it with FEM simulation. Two theoretical relations describing the elastic modulus decline ratio and the unit cell toughness promotion are developed as functions of the geometrical parameters and the material parameters, respectively. Moreover, we demonstrate a strategy to adjust the unit cell stiffness and structure toughness by typical geometrical parameter η (the "hard plate" overlapped length to non-overlapped length ratio) and typical material parameter α (the "joint part" tensile modulus to "hard plate" tensile modulus ratio).
Based on the above theatrical analysis and FEM simulation, we can draw the following conclusions. Firstly, the breakage of the "joint part" within the unit cell, while phase changing, does not mean biomimetic composite structure failure. The strength property of the structure should be decided at "phase 1", by the tensile strength of "hard plate" and the shear strength of "shear part" together. Therefore, with the proper choice of materials for the "hard plate" and "shear part", when the breakage of Region 5 occurs, the unit cell can realize extra energy dissipation and stiffness changing without loss of strength. Secondly, as η increases, both E 0 (the effective modulus at "phase 0") and E 1 (the effective modulus at "phase 1") would increase with an attenuated increase rate and finally approach a limiting value, derived as Equations (16) and (17). However, as α increasing, only E 0 will be influenced, i.e., increase with an attenuated increase rate and finally approach a limiting value. Therefore, the structure stiffness before phase changing will restrict the material selection. Thirdly, E 1 /E 0 (the effective modulus decline ratio) and T cell (the normalized "structure toughness" of a unit cell) are two unique target parameters for our design. They will respectively decide the lower and upper limiting value of η, derived as Equations (21) and (30), which are two constraint conditions for structure configuration design.
To recapitulate, the stiffness changing and toughness promotion of our biomimetic structure can be precisely achieved with theoretically calculated structure parameters, i.e., by quantitatively tailoring the "joint part" breaking process within each unit cell. The theoretical analysis was verified by the FEM simulation. Moreover, the simulation results provide more detailed insights into the mechanical properties of the unit cell, such as the detailed stress and strain field output and the specific breaking process of Region 5 during phase changing. Therefore, combining the theoretical predication and the FEM verification, it is possible to adjust the properties of our biomimetic composite structure.
With superior toughness and tunable stiffness, our biomimetic composite structure can serve as a guideline in designing novel load-bearing structures.
Future extension of this study can involve investigating and designing the dynamic properties of the unit cell, under dynamic load. Taking the unit cell as the basic building block, a hierarchical structure can be constructed with different unit cell arrangement modes. Then, the influence of the unit cell properties, as well as the arrangement modes on the whole structure dynamic response, can be studied.