Analysis and verification of a biomimetic design model based on fish skin

Many biological architectures are Bouligand structures, which comprise uniaxial fiber layers stacked in a periodic helical arrangement and are characterized by high damage resistance. As an effective flexible protective structure, fish skin is a Bouligand structure that protects the body while ensure flexibility during swimming and predation. In this paper, an analytical model inspired by fish skin is established based on previous studies, and the parameters for describing crack growth are determined. Then, mathematical expressions for the local stress intensity factors and plastic zone are used to predict how the helical stacking angle α influences the crack propagation. The results show that crack deflection and twisting improve the fracture toughness of the composite structure greatly, with the optimal fracture toughness being that for α = 60° – 70°. Moreover, biomimetic flexible composite structures inspired by fish skin are produced using silicone and Kevlar fibers. Scanning electron microscopy is used to observe the cross-sectional morphology of the composite structures, showing that the interfaces between the silicone and Kevlar fibers are highly compact. Results from experimental impact tests agree well with the predicted results.


Introduction
Protection of objects against impact loads has gained much attention due to increasing awareness of security issues. According to previous studies [1,2], impact leads to severe damage to protected objects through two ways: firstly, an impactor may penetrate the targets and cause internal damage; secondly, the impact energy and high stress could cause severe deformation and disintegration of the objects. Therefore, the flexible protective structures which can resist the penetration and able to absorb the impact energy are an effective model to resist impact.
To meet such requirements, the biological structures can provide an inspiration for designing anti-impact structure. Bouligand structures are simple and effective biological protective structures that comprise uniaxial fiber layers stacked in a periodic helical arrangement [3]. They are found in different biological tissues, such as fish skin [4][5][6], fish scales [7][8][9], the elytra of Chrysina aurigans and Chrysina limbata [10], the cuticles of Callinectes sapidus and Homarus americanus [11], and the dactyl club of Odontodactylus scyllarus [12,13]. The Bouligand structure gives biological tissue excellent mechanical properties. For example, the dactyl club of O. Scyllarus can endure tremendous impact forces (0.41.5 kN) [14], and arapaima scales are an ideal light and strong biological material, with a single scale bearing a penetration force of 150N [15]. Also, carbon-fiberreinforced polymer composites have been designed based on the Bouligand structure, and the produced material had improved fracture resistance [16,17].
As a flexible material, the mineralization degree of the helically stacked fibers in fish skin is lower than those of the fibers in fish scales and the exoskeleton of beetles, thereby resulting in good deformability. Fish skin is an ideal flexible protective structure that is effective at reducing resistance during swimming and protecting against attack (such as bite). Previous studies have shown that fish skin has impressive mechanical properties in tension and penetration tests and revealed that the tension resistance of fish skin is very close to the tendon tissue [4,18,19]. Also, the helical stacking angle of the fibers in fish skin changes from anterior to posterior, meaning that so do the mechanical properties [4,5]; the likely purpose of this change is to offer different protection functions. The fibers in fish skin are wrapped in a collagen matrix that is mechanically weaker than the fibers [4,20], but the collagen matrix stores energy and provides structural support, thereby making the composite structure more flexible.
Herein, we study the growth of cracks in an elastic composite structure inspired by fish skin. We assumed that the fibers themselves do not fracture and the crack grows along the matrix material. In the composite structure, the fiber diameter and spacing and the helical angle of the layers are constant (see the subsection entitled 'Crack deflection and twist'). The Bouligand structure decreases the anisotropy of the elastic composite structure, meaning that the material can be approximated as being in-plane elastically isotropic [21]. The analysis is simplified greatly by replacing anisotropic materials with isotropic and homogeneous ones. Based on previous studies [22,23], we build an analytical model to investigate crack growth, deflection, and twisting in a flexible Bouligand structure. By analyzing the local stress intensity factors in the local coordinate system of the crack front as a function of the applied stress intensity factors in the global coordinate system of the initial crack (see the subsection entitled 'Model formulation'), the analytical model offers a preliminarily understanding of the crack growth process and toughening mechanisms in a flexible Bouligand structure. We also derive an expression for the plastic zone at the crack tip to study further how the degrees of crack deflection and twisting affect the fracture toughness of the composite structure. Finally, to test and validate the analytical model, we prepare biomimetic flexible protective structures inspired by fish skin by using elastic silicone (Dow Corning Sylgard 184) and Kevlar fibers. Compared with previous studies [22,23], the crack propagation paths in our analytical model are defined as continuous curves, which is closer to the curved characteristic in the Bouligand structure and leads to more accurate analysis results. Furthermore, the effect of the helical angles in the helically stacked fibers on the size of the plastic zone is determined, which can help us to select the optimal helical angles in the preparation of fiber reinforced flexible composite structures. The analytical model is intended to provide valuable guidance for understanding the fracture toughness mechanisms in a flexible Bouligand structure, as well as offering a theoretical support for the design of flexible protective structures against impact.

Materials and methods
2.1. Crack deflection and twisting Figure1(a) shows the structural characteristics of the helically stacked fibers in the Bouligand structure, which are described by the fiber spacing (d 1 , d 2 ) and the helical angle (α) of adjacent layers. In the present analytical model, crack growth has the following features: (i) the initial crack is straight and parallel to the fibers; (ii) the crack propagates in a finite field; (iii) the fiber structure impedes the crack growth. Consequently, the helically stacked fibers deflect the crack propagation, and the initially straight crack front twists during propagation, as shown in figure1(b).
The origin of the coordinate system is defined and located at the front of the notch ( figure 1(b)). X-axis is the direction of crack propagation, the Y-axis represents the deflection height of the crack, and the Z-axis is parallel to the initial crack front and serves as the width of the crack. Therefore, the twisting crack is represented mathematically as: ( )/ denotes the fiber spacing and α is the helical angle of adjacent layers.

Model formulation
Analyzing the crack deflection, twisting, and propagation on specific paths offers better understanding of how the crack growth mode influences the fracture toughness of the structure. Faber and Evans [22] proposed a theoretical model with which to derive the local stress intensity factors at the twisting crack front according to the applied stress intensity factors (or K-field) by using LEFM. They assumed that the crack propagates in only the linear elastic isotropic matrix material under uniaxial applied mode-I loading. In the theoretical model, the paths of crack propagation are defined as a straight line, as shown in figure2. The model defines two coordinate systems, namely (i) the global coordinate system (x, y, z) at the initial crack front and (ii) the local coordinate system (x′, y′, z′) at the twisting crack front. The z and z′ axes are parallel to the initial crack front and the twisting crack front, respectively. Figure2(b) shows the deflection angle β and the twist angle f, where β denotes the deviation of the twisting crack from the initial crack front and f represents the angle between the z-and z′-axes. Because the twisting crack front is a straight line, f is always constant at any position of the crack front. Under K I loading and plane strain conditions, the local stress intensity factors ¢ k , II and ¢ k III at the twisting crack front can be expressed as [22,23]:    sin  2  cos  2  cos 1 3 sin  2  sin  2  cos  2 cos , 2 cos  2  sin cos cos  2  2  sin  2  cos  2  sin cos 3 cos  2  2 , 2 where K I is the applied stress intensity factor under mode-I loading. Furthermore, Suksangpanya et al built an analytical model based on the Faber-Evans model and defined the path of crack deflection as a continuous curve [23]. However, the twisting crack front in the analytical model is also defined as a straight line, resulting in a constant twist angle f. It is assumed that the crack only grows in the isotropic elastic material under mode I loading. By using the asymptotic stress field [24], the local stress intensity factors at the twisting crack front are derived under plane strain conditions, as shown in Equation II and ¢ k III denote the local stress intensity factors at the twisting crack front, K I is the applied stress intensity factor under mode-I loading. β and f represent the deflection angle and twist angle, respectively. The model is used to analyze the local stress intensity factors and to simulate the local energy release rate at the twisting crack front by ABAQUS. The results show that the accuracy of the model is higher than that of Faber-Evans model. b However, although the Faber-Evans model and the Suksangpanya's model provide a good theoretical framework for understanding how crack twisting affects the fracture toughness of composite structures, the defined straight-line characteristic at the twisting crack front in their models differs from the curved characteristic in the Bouligand structure [13]. Therefore, we develop an analytical model based on their works, as shown in figure 3. To explain accurately the local stress intensity factors at the twisting crack front, four coordinate systems are defined (figure 3(a)), namely the global coordinate system (x, y, z) and its equivalent cylindrical system (r, θ, z), and the local coordinate system (x′, y′, z′) and its equivalent cylindrical system (r′, θ′, z′), where r and r′ are in the xy-and x′y′-planes, respectively. In the coordinate systems, the z-axis is parallel to the initial crack front, the z′-axis is tangent to the curve of the twisting crack front, and the x′-axis is defined as the crack growth direction under mode-I loading. Consequently, the local coordinate system (x′, y′, z′) changes with the position of the twisting crack front (figure3(c)).
To analyze the relationship between the local coordinate system (x′, y′, z′) and the global coordinate system (x, y, z), we also define the deflection angle β′ and the twist angle f ( figure 3(b)). Because the z′-axis can be obtained by rotating the z-axis around the x-axis by angle f, z′ can be expressed mathematically as: where R x is the rotation matrix around the x-axis. Because f is the rotation angle of the twisting crack front relative to the initial crack, according to the definition proposed by Suksangpanya et al [23], f d dx / is the rate of change of f along the x-axis. According to the Bouligand structure in figure 1 / Therefore, f can be expressed as a linear function of the position along the z-axis: Meanwhile, β′ is the angle between the y and y′′ axes, and the y′′-axis is the projection vector of the y′-axis onto the xy-plane ( figure 3(b)). Here, β′ is expressed as: The main difference between the present analytical model and those by Faber-Evans and Suksangpanya et al is the paths of crack deflection and twisting. In this paper, the paths of crack deflection and twisting are defined as continuous variable curves, which is more consistent with the crack propagation paths in the Bouligand structure.

Model derivation
In the process of crack growth, the local stress field changes with the position of the twisting crack front. Therefore, we establish an element at a point in the twisting crack front, assuming that the size of the element is very small and the stress on each surface of the element is uniform, as shown in figure4. σ and τ denote the normal stress and the shear stress of the element, respectively. The analysis of stress state of the element can help us to study the local stress field at a point in the crack front. Furthermore, the element is simplified to a twodimensional plane ( figure 4(b)), s ¢ x and t ¢ ¢ x y denote the normal stress and shear stress on the element plane perpendicular to the x′-axis, respectively; s ¢ y and t ¢ ¢ y x denote the normal stress and shear stress on the element plane perpendicular to the y′-axis, respectively;s ¢ r and t q ¢ ¢ r denote the normal stress and shear stress on the element plane perpendicular to the r′-axis, respectively; s q¢ and t q¢ ¢ r denote the normal stress and shear stress on the element plane parallel to the r′-axis, respectively. The value of t ¢ ¢ x y is equal to that of t ¢ ¢ , y x and the value of t q ¢ ¢ r is equal to that of t q¢ ¢ . r The two-dimensional plane is divided into two parts by section b′c′ and the section a′′ b′c′d (figure 4(c)) is analyzed by using the principle of force balance:  where dS is the cross-sectional area of b′c′. Equation (7) can be simplified as: 2   cos  2  1  sin  2  sin  3  2  2  sin  2  2  cos  2  cos  3  2 , 9 x  By introducing equation (9) into equation (8), the local stress field s t ¢ ¢ [ ] / at a point in the twisting crack front can be expressed: 2   cos  2  1  sin 2 sin  2  cos  3  2  sin  2  1  cos2 sin  3  2   2  sin  2  1  cos  2  cos  3 2   sin  2  2 sin  3  2  cos  2  cos  3  2  tan  cot   2  2 sin  2  1  cos  2  cos  3  2  cos  2  1  sin  2  sin  3  2 tan cot , 10 where K , I K II and K III are the applied stress intensity factors. In this analytical model, only the opening fracture mode (mode I) is considered, so K I ≠ 0, K II =K III =0. According to the relationship between the cylindrical coordinate systems (r′, θ′, z′) and (r, θ, z), the local stress field t s¢ ¢ [ ] / at q¢ =0°is equal to the global stress field t s [ ] / at q =b¢ rotated around the x-axis by angle f, which can be written mathematically as: The functional relationships between the local and applied stress intensity factors can be expressed as:

Specimen preparation
We prepared biomimetic flexible composite structures using silicone (Dow Corning Sylgard 184) and Kevlar fibers. The silicone comprised its liquid silicone and a curing agent. In the curing process, the mixing ratio of the liquid silicone and the curing agent was 10:1, and the mass of the mixture was 64.89g. The process of sample preparation was as follows. First, we put the mixture in a blender and stirred it for 25 min, then we used a vacuum pump for 30 min to remove any bubbles in the mixture. Next, we divided the mixture into (usually) seven parts, poured one part into a rectangular mold, and repeated the bubble-removal process (30 min). We then put the rectangular mold into an electro-thermostatic blast oven at a constant temperature of 100°C for a curing time of 35 min. Next, we laid Kevlar fibers on the solidified silicone, covered them with another part of the mixture, and repeated the bubble-removal and curing processes under the same aforementioned conditions. Each sample contained five layers of Kevlar fibers and had a constant helical stacking angle. We selected the six helical stacking angles of α′=50°, 60°, 65°, 70°, 80°and 90°according to the results obtained from the analytical model (see the subsections entitled 'Features of crack deflection and twisting' and 'K-field and plastic zone in analytical model'). Figure5 shows the arrangement of the fibrous layers for α′=50°sample. The initial fibrous layer aligned with the y-axis was denoted as 0°, and the subsequent fibrous layers were rotated in the counter-clockwise direction. The angle between the 5th layer and the y-axis is 20°( figure 5(c)). Also, the angles between the 5th layer and the y-axis for α′=60°, 65°, 70°, 80°and 90°specimens are 60°, 80°, 100°, 140°and 0°respectively ( figure 6). Figure6 shows the prepared composite samples, the size of which was 120mm×40mm×13mm. There were three specimens for each helical stacking angle α′.

Experiments
Before impact tests, we chose a composite structure (α′=60°) to cut open and to observe the cross-sectional microstructure by means of scanning electron microscopy (SEM). The objective here was to study the bonding state of the two materials (silicone and Kevlar fibers). The SEM (QUANTA 650, FEI: Portland, USA) test was performed at 20 kV.
To verify the analytical model that predicts how the helical stacking angle α′ affects the mechanical properties of the flexible Bouligand structure, we performed impact tests on these composite samples. The main parts of the impact equipment (DIT302A-TS; Wance; China) are its drop-hammer tower, drop hammer, load sensor, velocity measuring device, and pneumatic fixture, as shown in figure7. The mass of the drop hammer was 5kg, and the impact head was a hemisphere with a diameter of 20mm. The impact energy and speed were controlled by adjusting the hammer's drop height. The maximum impact energy and impact height were 300J and 6m, respectively.

Features of crack deflection and twisting
The deflection angle b¢ and twist angle f can be expressed using the defined Cartesian coordinate systems and geometric structure parameters, where b¢ is the deviation of the projection-vector y′′-axis from the initial crack front, and f is the deviation of the twisting crack from the z-axis. Therefore, b¢ and f can be used to describe the shape characteristics of the crack propagation. The relationships between + y x d d ( ) / characterized byf and b¢, respectively, for a =50°, 60°, 70°, 80°and 90°. These curves show the degree of crack deflection and twisting in the propagation process for different values of the stacking angle α. They show that the curve depends on α, indicating that the deflection and twisting state of the propagating crack differs with α. With increasing α, the degree of crack deflection and twisting increases, therefore crack deflection and twisting play a decisive role in the mechanical properties of composite structures by transforming the paths of crack growth.

Microstructural characterization
The cross-sectional microstructure of the composite structure (α′=60°) is shown in figure10. As shown, the interface between the silicone and Kevlar fibers is compact, without pores and cracks (figures10(a) and (c)). Moreover, the spaces between the fibrils are fully filled with silicone (figure10(d)), thereby improving the mechanical properties of the composite structure. The distance between fiber layers is approximately 1.35mm, the thickness of a fiber layer is approximately 0.75mm (figure10(a)), and the diameter of the fibrils is approximately 10μm (figure10(d)). According to the literature [20], the fibrils in fish skin are wrapped in a collagen matrix, and this type of biological structure is effective against applied forces such as bites and impacts through the excellent toughness of the fibers and the large deformation of the collagen matrix. The biomimetic flexible composite structures herein are close to the Bouligand structure in fish skin, thereby offering further understanding of the fracture mechanisms of the flexible Bouligand structure through mechanical tests.

Impact tests
The results of the impact tests are presented in figure11 and table1. In figure 11, the meaning of different color curves denotes that there are three specimens for each helical stacking angle α′. Table1 shows that the samples with α′=65°had the highest peak load and impact resistance, whereas those with α′=50°had the lowest peak load and impact resistance. The impact resistance of the samples decreased gradually upon increasing α′ from 65°to 90°, but with a decreasing rate of decrease. The maximum impact energy and specific absorption energy (M) of the biomimetic structures with α′ =60°, 65°, and 70°are ordered as M 65°> M 60°> M 70°. The optimal helical stacking angle obtained by the impact tests is consistent with the results from the present analytical model, thereby indicating that the model is valuable and instructive for analyzing flexible Bouligand structures with isotropic matrix materials. The fracture behavior of the specimens in figure11 shows that they failed mainly in delamination fracture under impact loading, with the crack propagating along the silicone matrix. The results show that helically stacked fibers are effective at improving the fracture toughness of composite structures. The reason why the load-displacement curves in figure11 show significant fluctuations is that the helically stacked fibers are effective at dispersing the impact force and impeding the drop hammer from advancing along the thickness direction of the specimens, resulting in fluctuating energy transfer and thus fluctuating loading. Under the impact force, most of the impact energy is absorbed by the silicone structure through deformation, and the stress concentration is prevented by the helically stacked fibers changing the crack propagation paths in the silicone matrix. This collaborative protection mechanism provides unique and superior fracture resistance to the biomimetic composite structures. The impact energy propagates mainly in two directions in the composite structures, namely (i) in the thickness direction of the sample and (ii) along the  silicone matrix to the ends of the sample. Furthermore, during an impact, the fibrous layers of the upper surface of the sample experience compressive stress, those in the middle of the sample experience mainly shear stress, and those in the lower surface of the sample experience mainly tensile stress [26]. Consequently, the fibrous layers at different positions in the composite structures have different stress-mode responses in the impact process.

Conclusion
In this study, we developed an analytical model based on previous studies with which to study crack growth in a Bouligand structure inspired by fish skin. Moreover, we compared the present model with the theoretical model from Faber and Evans [22] and verified the former by impact tests on biomimetic composite structures. Mechanical evaluation of the impact energy, specific absorption energy, and fracture toughness of the biomimetic composite structures showed that the helically stacked fibers make a critical contribution to the damage resistance of the structures, and the test results basically agree with the present model. We draw the following specific conclusions: 1. The deflection angle b¢ and twist angle f were established as parameters in this analysis model (figure3) and were used to describe the crack propagation characteristics. We assumed that the crack propagates only in the isotropic elastic matrix material, and then we used micro-mechanics to analyze our model. The results show that the helically stacked fibers force the crack to take on a more curved shape, and the crack deflection and twisting increase with the helical stacking angle. Furthermore, we used the mathematical expressions for the local stress intensity factors and plastic zone to predict how the crack deflection and twisting affects the flexible Bouligand structure. Finally, the analytical model gave the optimal helical stacking angle to be in the range of 60°-70°.
2. Biomimetic flexible composite samples inspired by the Bouligand structure in fish skin were prepared using silicone and Kevlar fibers, and the cross-sectional microstructure of the samples was observed using SEM. The results showed that the interface between the silicone and Kevlar fibers was tight and without pores. Moreover, the fibrous layers were arranged regularly, and the spacing between them was basically identical. This superior bonding helps greatly to enhance the mechanical properties of the biomimetic composite structures.
3. The analytical model in this study was validated with impact tests on the biomimetic samples, and the tests showed clearly that the biomimetic composite structures were characterized by excellent impact resistance. Delamination fracture was the main failure mode of the composite structures. Moreover, the samples with α′=65°had the best impact resistance, which is consistent with the prediction of the analytical model.
In this work, we have provided guidelines for analyzing and understanding the fracture mechanism of a flexible Bouligand structure, as well as a possible novel and valuable idea for designing flexible protective structures. This is an ongoing project, and next work is to apply the bio-inspired composite structure to the field of engineering, then to validate its defense capability, energy absorption efficiency and penetration resistance by comparing with other anti-impact structures (such as sandwich structure, foamed metal).