## Abstract

Composites comprising a high-volume fraction of stiff reinforcements within a compliant matrix are commonly found in natural materials. The disparate properties of the constituent materials endow resilience to the composite, and here we report an investigation into some of the mechanisms at play. We report experiments and simulations of a prototype laminated composite system comprising silicon layers separated by polymer interlayers, where the only failure mechanism is the tensile fracture of the brittle silicon. Two failure modes are observed for such composites loaded in three-point bending: failure under the central roller in (i) the top ply (in contact with the roller) or (ii) the bottom ply (free surface). The former mode is benign with the beam retaining load carrying capacity, whereas the latter leads to catastrophic beam failure. Finite element (FE) simulations confirm this transition in failure mode and inform the development of a reduced order model. Good agreement is shown between measurements, FE simulations, and reduced order predictions, capturing the effects of material and geometric properties on the flexural rigidity, first ply failure mode, and failure load. A failure mechanism map for this system is reported that can be used to inform the design of such laminated composites.

## 1 Introduction

Layered composite beams and plates comprising stiff reinforcements and compliant matrices are widely found in engineering structures (e.g., natural and engineered wood), aerospace vehicles (e.g., fiber-reinforced composites), nano-materials (e.g., graphene/MoS_{2} composites), and biological organisms (e.g., mollusk and crustacean shells). These composites offer desirable combinations of high stiffness, strength, and toughness [1]. For instance, the “brick-and-mortar” bio-composite, nacre, has a toughness 10^{3} times greater than the brittle and stiff aragonite calcium carbonate platelets, which constitute ∼95% of its volume [2]. Studies have revealed that a compliant biopolymer “adhesive” surrounding an architected array of platelets creates this potent toughening mechanism [3–5]. In high performance composites, such as continuous fiber-reinforced laminates and wood, the compliant matrix is instead intercalated between continuous layers of stiff, yet brittle, reinforcement. Here we investigate the flexural rigidity, load bearing capacity, and failure modes of such continuous layered composites under three-point bending, revealing the effects of geometry and material properties on performance.

For layered composites, comprising layers with disparate material properties, homogenization theories, such as classical Euler–Bernoulli, Timoshenko, and higher order beam theories, are unable to accurately capture the stress and strain state within the beam [6]. In contrast, layer-wise theories account for the heterogeneous nature of the beam, but tend to be computationally complicated, with the number of kinematic parameters increasing with the number of layers, therefore limiting their application to beams comprising few plies. For instance, Moussiaux et al. [7] and Adams and Weinstein [8] presented models of a three ply beam, comprising two continuous stiff outer reinforcement layers and a single inner matrix layer. Nguyen et al. [9] developed an exact solution for a three layer Timoshenko beam, while Škec and Jelenić [10], Ranzi [11], and Sousa and da Silva [12] proposed the solutions for multi-layer composite beams. Most of these studies focused on formulations for finite element (FE) implementation, providing little insight into the performance design space.

Simplified layer-wise theories have been presented that are able to capture the discrete nature of laminated composites [13,14]. These theories model the non-monotonic axial displacements through the thickness of the beam as a piecewise linear function to limit the number of kinematic variables. Similarly, Bathe et al. [15] used a variational energy method, assuming a linear variation of shear between layers, to capture the bending of cytoskeletal actin bundles, modeled as stretchable reinforcing fibers connected by cross-linked shear components. Good agreement with measured data was shown, and it resulted in a simple formulation showcasing the dependencies of the key material properties.

The present paper examines the response of a laminated beam comprising compliant matrix plies intercalated between stiff elastic-brittle reinforcement plies subjected to three-point bending (Fig. 1). We use this prototype system to explore the resistance to flexure (i.e., flexural rigidity), the mode of first ply failure, and the failure load. The paper first presents an examination of the material and geometric properties controlling the behavior of the beam using FE analysis. This analysis will reveal a switch in first ply failure location between the bottom (free surface) ply and the top (in contact with the center load) ply. The FE predictions are then validated against a model experimental system. Finally, a reduced order model will reveal the mechanism responsible for this change in failure mode and used to construct a design map.

## 2 Methodology

### 2.1 Problem Definition.

We considered a laminated beam comprising alternating linear-elastic plies of a stiff reinforcement and a compliant matrix, with perfect bonding between the plies (Fig. 1). The outer two plies of the lay-up were reinforcement layers, resulting in a beam with *N* and *N* − 1 reinforcement and matrix plies, respectively. Individual plies were identified by the integer *n*, which had the range −(*N* − 1) ≤ *n* ≤ (*N* − 1) with *n* taking even and odd values for reinforcement and matrix plies, respectively. Each constituent was assumed to be linear elastic with Young's modulus *E*_{i}, Poisson ratio *ν*_{i}, and thickness *t*_{i}, with the subscript “*i*” taking the value “*r*” and “*m*” for “reinforcement” and “matrix,” respectively. The beam was of uniform cross section with a width *b* and thickness *H*, such that *H* ≡ *Nt*_{r} + (*N* − 1)*t*_{m}. The beam of overall length *l* = 1.2*L* was supported at two points separated by a length *L* in the *x* direction and loaded in the negative *z* direction at its mid-span. Application of a load *F* resulted in a mid-span deflection *d*, as measured at the point of the applied load.

The parameter space for this problem was cast into four independent non-dimensional variables, namely, Young's modulus ratio $E\xaf\u2261Em/Er$, the beam slenderness ratio $L\xaf\u2261L/H$, the ply thickness ratio $t\xaf\u2261tm/tr$ and the number of reinforcement plies *N*, such that the reinforcement volume fraction $\varphi r=[1+t\xaf(1\u22121/N)]\u22121$. FE simulations revealed that the influence of the Poisson ratios *ν*_{r} and *ν*_{m} was minor, and hence, the effect of these parameters will not be addressed in detail.

### 2.2 Finite Element Study.

Two-dimensional FE calculations of the composite beam, described above, were performed with the software package abaqus (Dassault Systemes Simulia Corp., Providence, RI) in a small strain setting, using plane stress elements (CPS4R). The materials were modeled as isotropic, linear-elastic solids, with *ν*_{r} = 0.22 and *ν*_{m} = 0.45 (matching the experimental materials detailed in the subsequent section). The beam was mirrored along the *x* axis at the beam mid-span, with one support at *x* = 0 and the mid-span at *x* = *L*/2 (Fig. 1). Each ply had four elements through the thickness except for the outer reinforcement plies, which had 60 elements. The loading rollers were modeled as analytically rigid surfaces, with radii *R* = *L*/20 with frictionless contact between the rollers and the beam. Brittle tensile failure of the reinforcement was specified as occurring when the longitudinal strain *ɛ* at any location within a reinforcement ply reached the tensile rupture strain *e*_{r}.

### 2.3 Experimental Study

#### 2.3.1 Materials.

The reinforcement plies were sectioned from polycrystalline silicon wafers (PV Crystalox, Abingdon, UK), measuring *t*_{r} = 0.23 mm in thickness and dimensions 156 mm × 156 mm in-plane. The following three thermoplastic films, with differing elastic moduli, were used for the matrix plies: a low-density polyethylene (LDPE, Goodfellow Cambridge Ltd., Huntingdon, UK), and two polyurethane (PU) grades (TF-310 and TFL-2EA from Permali Gloucester, Gloucester, UK). The films were of thickness *t*_{m} = 0.05 mm.

The mechanical properties of the silicon wafer were measured by placing a 9 mm wide sample in three-point bending, using an outer roller separation of 50 mm. The silicon was elastic-brittle (Fig. 2(a)), with Young's modulus of *E*_{r} = 160 ± 5 GPa, failure strength *s*_{r} = 143 ± 10 MPa defined as the maximum tensile stress in the beam (the bending stress in the beam is plotted on the right-hand *y*-axis in Fig. 2(a)), and failure strain of *e*_{r} ≡ *s*_{r}/*E*_{r} = 8.8 × 10^{−4} (averaged from ten repeat tests). For the polymeric films, the material properties were obtained from tensile tests, following ASTM standard D882 [16]. Rectangular coupons, measuring 10 mm wide and with a 100 mm gauge length, were loaded in uniaxial tension at a nominal strain rate of 10^{−2} s^{−1}. Figure 2(b) shows representative responses of the tensile tests (five replicate tests were conducted for each material). The polymers were linear elastic to $\u223c5%$ strain, with failure strains of $>40%$. Table 1 provides the elastic properties for all of the materials used in the study.

Constituent | Material | t_{i} (mm) | E_{i} (GPa) | $E\xaf=Em/Er$ | $\nu i$ |
---|---|---|---|---|---|

Reinforcement | Silicon | 0.23 | 160 | – | 0.22 |

Matrix | LDPE | 0.05 | 0.13 | 8.1 × 10^{−4} | 0.45 |

TF-310 | 0.05 | 0.05 | 3.1 × 10^{−4} | 0.45 | |

TFL-2EA | 0.05 | 0.02 | 1.5 × 10^{−4} | 0.45 |

Constituent | Material | t_{i} (mm) | E_{i} (GPa) | $E\xaf=Em/Er$ | $\nu i$ |
---|---|---|---|---|---|

Reinforcement | Silicon | 0.23 | 160 | – | 0.22 |

Matrix | LDPE | 0.05 | 0.13 | 8.1 × 10^{−4} | 0.45 |

TF-310 | 0.05 | 0.05 | 3.1 × 10^{−4} | 0.45 | |

TFL-2EA | 0.05 | 0.02 | 1.5 × 10^{−4} | 0.45 |

Note: The modulus ratio is defined as $E\xaf\u2261Em/Er$ with the Poisson ratio *ν*_{i} values taken from Ref. [17].

#### 2.3.2 Laminated Beam Fabrication.

We used a CO_{2} laser cutter (Universal Laser Systems GmBH, Vienna, Austria) to cut 73 mm × 70 mm sections from the silicon wafers. Polymer plies were then layered between plies of silicon and the ply stack loaded through the thickness to 50 N and placed for 20 min in an oven heated to 145 °C and 175 °C for the LDPE and PU systems, respectively, so as to adhere the polymer to the Si. A circular saw, equipped with a diamond wafering blade, removed an outer 5 mm border from the composite plate to discard any potential heat effected zone from the laser cutting process. The saw was then used to cut *b* = 9 mm wide beams from the composite plate. Finally, the cut faces were wet ground using progressively finer SiC abrasive papers. Table 2 provides the final beam dimensions.

Matrix | L (mm) | H (mm) | $E\xaf$ | N | $t\xaf$ | $L\xaf$ |
---|---|---|---|---|---|---|

LDPE | 50 | 1.35 | 8.1 × 10^{−4} | 5 | 0.21 | 37 |

TF-310 | 50 | 1.35 | 3.1 × 10^{−4} | 5 | 0.21 | 37 |

TFL-2EA | 50 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 37 |

TFL-2EA | 26 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 19 |

TFL-2EA | 12 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 9 |

TF-310 | 50 | 3.04 | 3.1 × 10^{−4} | 11 | 0.21 | 17 |

TF-310 | 50 | 4.72 | 3.1 × 10^{−4} | 17 | 0.21 | 11 |

Matrix | L (mm) | H (mm) | $E\xaf$ | N | $t\xaf$ | $L\xaf$ |
---|---|---|---|---|---|---|

LDPE | 50 | 1.35 | 8.1 × 10^{−4} | 5 | 0.21 | 37 |

TF-310 | 50 | 1.35 | 3.1 × 10^{−4} | 5 | 0.21 | 37 |

TFL-2EA | 50 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 37 |

TFL-2EA | 26 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 19 |

TFL-2EA | 12 | 1.35 | 1.5 × 10^{−4} | 5 | 0.21 | 9 |

TF-310 | 50 | 3.04 | 3.1 × 10^{−4} | 11 | 0.21 | 17 |

TF-310 | 50 | 4.72 | 3.1 × 10^{−4} | 17 | 0.21 | 11 |

#### 2.3.3 Test Methodology.

The beams were mounted for three-point bending (as illustrated in Fig. 1) onto steel rollers, with radii *R* ≈ *L*/20. A model 5544 Instron, equipped with a 5 or a 500 N load cell, loaded the beams along their mid-span at a constant displacement rate of 2 mm min^{−1}. The mid-span deflection *d* was measured by monitoring the displacement of the central roller using a laser extensometer and microphones were used to pick up acoustic emissions during the bending experiment.

## 3 Results

### 3.1 Finite Element Predictions of the Stress Distribution and Failure Modes.

An FE parametric study was performed to assess the stress distribution and active failure modes within the composite beam. Figure 3(a) shows the distribution of normalized longitudinal stress *σ*/*s*_{r} for a beam with *N* = 5 reinforcement plies, a $t\xaf=0.2$ ply thickness ratio, a $L\xaf=20$ slenderness ratio, and a modulus ratio $E\xaf=10\u22122$. The distribution is shown for an applied three-point bending load *F* = *P* such that the maximum value of the reinforcement ply longitudinal stress *σ*_{max} attains its failure value *s*_{r} = *E*_{r}*e*_{r}, where *e*_{r} = 8.8 × 10^{−4}. A magnified view near the central roller illustrates the non-monotonic variation of the longitudinal stress with *z*. The longitudinal strain *ɛ* profile below the central roller is also included in Fig. 3(a) and reveals a “zig-zag” strain distribution. Such a distribution has been previously reported for laminated fiber-reinforced composites [13,14]. Increasing the relative compliance of the matrix to $E\xaf=10\u22124$ and 10^{−6} (Figs. 3(b) and 3(c)) caused more pronounced zig-zagging due to greater tensile relaxation associated with shearing of the matrix.

A consequence of this change in behavior with decreasing $E\xaf$ was that the location of failure changes with $E\xaf$. For the $E\xaf=10\u22122$ beam, examination of Fig. 3(a) revealed that the bottom ply (ply most distal from the central roller) had the greatest tensile strain and thus failed first (recall that first ply tensile failure occurs when *ɛ*/*e*_{r} = 1). This is the expected failure mode for homogenous brittle materials in three-point bending. By decreasing the modulus ratio to $E\xaf=10\u22124$, the strain became tensile at the bottom of every reinforcement ply, with failure occurring simultaneously at the bottom and top plies of the beam (Fig. 3(b)). On the other hand, for the smallest modulus ratio $(E\xaf=10\u22126)$, first failure was restricted to the top ply (Fig. 3(c)). Hence, decreasing $E\xaf$ promoted a change in failure mode from reinforcement rupture of the rear ply to the top ply.

*D*and the first ply failure load

*P*of the beams. The flexural rigidity

*D*is defined as

*L*and stiffness

*F*/

*d*. Figure 4(a) plots the dependence of the normalized flexural rigidity $D\xaf\u2261D/(NDr)$, where $Dr\u2261Erbtr3/12$ is the flexural rigidity of an individual reinforcement ply, as a function of $E\xaf$. The flexural rigidity of the previously analyzed

*N*= 5, $t\xaf=0.2$, and $L\xaf=20$ beam remained nearly constant as $E\xaf$ decreased from 10

^{0}to 10

^{−2}. From $E\xaf=10\u22122$ to 10

^{−6}, $D\xaf$ decreased by more than an order of magnitude, with $D\xaf$ nearing unity at $E\xaf=10\u22126$. A value of $D\xaf=1$ implies that each reinforcement ply bends independently about its own neutral axis and is decoupled from the other plies. Intriguingly, values of $D\xaf<1$ were also observed and this is due to compliance of the matrix in the through-thickness direction of the beam. With decreasing $L\xaf$, the $D\xaf$ versus $E\xaf$ response shifts such that for a given $E\xaf$, the normalized rigidity $D\xaf$ decreases (Fig. 4(a)).

The normalized first ply failure load $P\xaf\u2261P/(NPr)$ versus $E\xaf$ relation is included in Fig. 4(b), where $Pr=2srtr2b/(3L)$ is the three-point bending failure load of a single reinforcement ply. For the *N* = 5, $t\xaf=0.2$, and $L\xaf=20$ beam, $P\xaf$ decreased with decreasing $E\xaf$ similar to $D\xaf$ (Fig. 4(b)). The FE results in Fig. 4(b) are marked to indicate the location of first ply failure and show that the transition from bottom to top ply failure shifted to smaller $E\xaf$ with increasing $L\xaf$. In summary, the FE results indicate that the flexural rigidity, failure mode, and failure load depend upon both $L\xaf$ and $E\xaf$.

### 3.2 Experimental Measurements and Observations.

The measured load per unit width *F*/*b* versus deflection *d* responses for *N* = 5 beams with $t\xaf=0.21$ are plotted in Fig. 5(a). The modulus ratio $E\xaf$ was controlled by changing the polymer and the slenderness ratio $L\xaf$ via the span length *L*. By contrast, Fig. 5(b) plots the *F*/*b* versus *d* responses, with $E\xaf$ and $t\xaf$ fixed at 3 × 10^{−4} and 0.21, respectively, and the slenderness ratio $L\xaf$ is controlled by varying *N*. In all cases, the load increased linearly with deflection, confirming that the systems were linear-elastic up to the point of failure. The stiffness (proportional to the slope) increased with increasing $E\xaf$ and decreasing $L\xaf$ as expected from the FE simulations.

Ultimate failure was catastrophic in all cases in Fig. 5. For instance, for the $E\xaf=3\xd710\u22124$ and $L\xaf=37$ beam (Fig. 5(a)), the load decreased by more than 90% after reaching its peak value. An in situ photograph (Fig. 6(a)) taken of this beam after failure showed that four of the five reinforcement plies had fractured under the central roller and only the top Si ply had survived. Conversely, all the matrix plies remained intact, indicating that ultimate failure was limited by the rupture of the Si plies. For the lowest slenderness ratio beams, an acoustic emission (AE) event occurred prior to catastrophic failure (the AE event is identified by a star in Figs. 5(a) and 5(b)). To better understand the origins of this AE, in a repeat test of the $L\xaf=11$ sample in Fig. 5(b), we unloaded the beam immediately after the AE event. Optical inspection of the outer surfaces revealed a crack in the top ply under the location of the central roller (Fig. 6(b)). A similar examination of all other samples revealed that top ply failure also occurred for the $E\xaf=10\u22124$, *N* = 5, and $L\xaf=9$ beam.

These results confirmed the FE prediction of a switch in the first ply failure mode from the bottom to the top ply. While both failure modes occurred under the central roller, the rear ply failure mode was catastrophic, leading to a cascade of ply failures with the load falling by up to 90%. On the other hand, the top ply failure mode was benign with the beam barely registering a change in stiffness after that first failure event.

Figure 5 also compares the FE results with the experimental results. The FE predictions of the failure load *P* and corresponding roller displacement (the load–displacement response is linear and therefore omitted for clarity) are shown by a symbol, where the open and cross-hatched circles represent failure by the bottom and top ply modes, respectively. Good agreement between FE predictions and measurements of the flexural rigidity, failure mode, and first ply failure load were observed for all configurations tested in this study.

## 4 Reduced Order Model for Laminated Beams

### 4.1 Homogenized Models.

*θ*related to the deflection

*w*at a longitudinal position

*x*along the beam via the relation

*θ*=

*dw*/

*dx*. Timoshenko extended this classical beam theory to account for shear deformation such as

*θ*=

*dw*/

*dx*+

*γ*, where

*γ*is the engineering shear strain across the cross section of the beam (Fig. 7(a)). The resulting displacement

*w*

_{T}of a Timoshenko beam is therefore given by the summation of the bending and shear contributions for a beam in three-point bending

*D*

_{EB}is the flexural rigidity for a Euler–Bernoulli beam,

*K*is the shear correction factor normally assumed to have a value of 5/6 for rectangular cross sections,

*A*≡

*bH*is the cross-sectional area of the beam, and

*G*is the shear modulus of the beam. In order to employ the Timoshenko beam model, we first need to determine the homogenized modulus

*G*and rigidity

*D*

_{EB}of the composite beam. The Ruess rules of mixtures give the shear modulus as

*G*

_{i}≡

*E*

_{i}/[2(1 +

*ν*

_{i})] (for

*i*=

*r*or

*m*) assuming isotropic elastic plies. Similarly, the Euler–Bernoulli flexural rigidity for the global bending of the composite beam can be estimated by summing the individual contributions of each ply as

*ɛ*and moment

*M*are related to the curvature

*κ*≡

*dθ*/

*dx*via

*ɛ*= −(

*dθ*/

*dx*)

*z*= −(

*M*/

*D*

_{EB})

*z*. For three-point bending, the maximum bending moment is at

*x*=

*L*/2 and thus the maximum longitudinal tensile strain within the beam follows as

*ɛ*

_{max}=

*FLH*/(8

*D*) at the bottom of the beam at mid-span. The first ply failure load

*P*

_{T}(where the subscript “

*T*” signifies prediction of the Timoshenko beam theory) is obtained by setting

*ɛ*

_{max}=

*e*

_{r}, and we write the normalized failure load as

Of course, in this homogenized model, failure is always predicted to occur on the rear face of the composite beam.

#### 4.1.1 Results.

Figure 4 compares the predictions using the Timoshenko model (TM) to the FE results. While the predictions for $D\xaf$ matched for large values of $E\xaf$, the TM under predicted $D\xaf$ for small $E\xaf$ (Fig. 4(a)). This was a consequence of the method used to homogenize the shear modulus which resulted in the incorrect limit of $D\xafT\u21920$ as $E\xaf\u21920$. Furthermore, Fig. 4(b) shows that the failure load prediction of the TM is also poor: the model predicts a negligible dependence of the failure load on $E\xaf$ and shows no dependence on $L\xaf$ contrary to the FE simulations (and measurements). This discrepancy is a consequence of the fact that the Timoshenko beam model a priori assumes a monotonic strain dependence in the *z* direction and thus cannot account for the discrete ply effects seen in Fig. 3 that are critical in setting failure.

### 4.2 Discrete Ply Model.

The inaccuracy of the Timoshenko beam model for a composite beam with a large phase contrast between plies suggests that bending and shear contributions must be considered at the ply level. In particular, the model needs to capture the two limits where the beam behaves either as a homogenous material when $E\xaf\u21921$ or as a stack of independent reinforcement plies when $E\xaf\u21920$. We shall first develop a model for the global bending of the composite beam (Fig. 7(a)) that captures these limits and then proceed to enhance this model to capture the switch in the failure mode from the bottom ply to the top ply by accounting for the local bending of the top ply.

#### 4.2.1 Global Beam Bending.

*n*th ply at a location $z~$ measured with respect to the local neutral axis of each ply (Fig. 1) as

*i*takes values

*r*or

*m*for a reinforcement or matrix ply, respectively, and $u\xaf(n)$ is the longitudinal displacement field within the

*n*th ply due to the global bending field. Moreover,

*θ*

_{r}=

*dw*/

*dx*is the rotation within a reinforcement ply, while

*θ*

_{m}=

*dw*/

*dx*+

*γ*

_{m}is that for a matrix ply, where

*γ*

_{m}is the shear strain in the ply. Displacement between plies implies

*w*(

*x*), which we determine using a Rayleigh–Ritz minimization scheme (see Appendix A for details). The outcome of this minimization is the full field solution $\Delta u\xaf(x)$ and

*w*(

*x*) from which we can determine quantities of interest. For example, the normalized flexural rigidity of the beam is given by the infinite series

*π*

^{2}

*η*)/(12

*α*)]

^{−1}in the last term. Similarly, the longitudinal strain at any position within the laminate can be expressed as

*x*=

*L*/2 and

*z*= −

*H*/2), leading to a normalized first ply failure load

#### 4.2.2 Combined Global and Local Bending.

The above solution does not account for local bending of the top reinforcement ply immediately under the central loading roller and hence cannot capture the transition of failure from the bottom ply to the top ply with decreasing $E\xaf$. In order to capture the local bending of the top ply, we employed the solution developed by Biot [18] for the bending of a beam resting on a semi-infinite elastic foundation.

*ɛ*

_{l}in the top reinforcement ply subjected to a transverse load

*F*as

#### 4.2.3 Results.

Predictions of the normalized flexural rigidity (23) using the discrete ply model (DPM) are included as solid lines in Fig. 4(a). Excellent agreement is observed between the DPM predictions and the FE simulations across all values of $E\xaf$ and $L\xaf$. Including the effect of local bending due to the through-thickness straining of the matrix plies allows the DPM model to be able to capture the reduction in the normalized flexural rigidity below unity. Equivalent comparisons between the FE and DPM model predictions (21) of the normalized failure load are included in Fig. 4(b). The results show that the DPM model was in good agreement with the FE results both in terms of the failure load and failure mode. DPM predictions of the longitudinal strain profiles below the central roller are also included in Fig. 3 along with the corresponding FE predictions at the same load level. The DPM captured the zig-zag strain distribution to a high level of fidelity. A more extensive investigation over a wide range of parameters demonstrated the fidelity of the DPM model when compared with FE simulations: details of this investigation are omitted here for the sake of brevity.

Through this analysis we showed that local ply bending can create a state of tensile strain at the bottom of the top ply immediately under the central roller. For sufficiently compliant matrices, this tensile strain can surpass the maximum tensile strain in the bottom ply due to global bending, and hence, the failure mode can switch from bottom ply failure to top ply failure. This top ply failure mode is benign, and hence, compliant matrices can give rise to the development of damage tolerance in laminated beams.

### 4.3 Failure Mechanism Map.

The DPM model is used to develop a failure mechanism map to provide insight into designing damage tolerant laminated beams. Four non-dimensional groups *N*, $L\xaf$, $t\xaf$, and $E\xaf$ parametrize the design space with specific combinations of these groups setting whether first failure occurs by failure of the top ply (the benign failure mode) or by failure of the bottom ply that results in catastrophic brittle failure. In Fig. 8, we include such a map with axes of number of plies *N* and beam slenderness ratio $L\xaf$.^{2} This map is interpreted as follows: Consider the case of $t\xaf=1$ and $E\xaf=10\u22124$, which is represented by the dashed-dotted line in the center (green-shaded) region. All designs (i.e., combinations of *N* and $L\xaf$) for this choice of $t\xaf$ and $E\xaf$ that lie to the right of the line will undergo catastrophic failure as the bottom ply fails first, while all designs to the left of the line will undergo the benign failure mode, where top ply fracture is the first failure event. We observe from the map in Fig. 8 that $t\xaf$ has a relatively small effect on the failure mode, while decreasing $E\xaf$ caused an appreciable shift whereby top ply failure is dominant over a larger design space.

The design map in Fig. 8 suggests that in general decreasing $E\xaf$ will switch the failure mode from the catastrophic bottom ply failure mode to the more benign top ply failure mode. However, a reduction in $E\xaf$ also increases the beam compliance and decreases the failure load. Hence, a trade-off exists in terms of maximizing the failure load and the damage tolerance simultaneously in such laminated beams.

## 5 Concluding Remarks

This study examined the three-point bending response of elastic laminated beam comprising alternating layers of a ductile matrix and a brittle reinforcement. Experiments with a prototype silicon/polymer beam revealed that first failure will occur by the fracture of a Si ply. The location of this failure is always under the central roller either in the top ply (in contact with the roller) or in the bottom ply (free surface). Bottom ply failure was catastrophic, leading to a cascade of ply failures and a substantial drop in the load carrying capacity of the beam. On the other hand, top ply failure was benign with the load continuing to increase after this first failure event. This transition in failure mode was confirmed using two-dimensional FE simulations. A reduced order model was developed to gain insight into the mechanisms and that model demonstrated that top ply failure occurred when the strain from local bending of the top ply due to the through-thickness straining of the compliant matrix plies exceeded the strain in the rear ply from the global bending of the beam. Good agreement was seen between experiments, FE simulations, and reduced order model predictions. Finally, a design map was presented showing the design space where different modes are operative. The map illustrates that the benign top ply failure mode can be activated mainly by (a) decreasing the matrix to reinforcement modulus ratio, (b) increasing the number of plies, or (c) decreasing the beam slenderness ratio. We anticipate that the analysis presented here can be extended to mixed-mode loadings (e.g., combined shear and bending) which will have applicability to additively manufactured biomechanical implants—this remains a topic for future studies.

## Footnote

We emphasize that our analysis assumes that each individual reinforcement ply behaves as an Euler–Bernoulli beam, although of course the shearing effects are included at the composite beam level. Thus, the model is only valid for reinforcement ply aspect ratios >10. The aspect ratio of each reinforcement ply in the composite beam is $L/tr\u2248L\xafN(1+t\xaf)$, and thus, the minimum value of *L*/*t*_{r} ≈ 30 in Fig. 8. This lies within the applicability of the Euler–Bernoulli beam assumptions employed in the model.

## Acknowledgment

The work was supported by the Office of Naval Research Grant N62909-15-1-N058 (Program manager, Dr. Judah Goldwasser). V.S.D. has benefitted greatly from inspiring and enjoyable discussions with R.M.M. and looks forward to that continuing for many years to come.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in this paper.

### Appendix A

#### Rayleigh–Ritz Solution for the Three-Point Bending of a Laminated Beam

*w*(

*x*) and $\Delta u\xaf(x)$ in the global bending solution of the DPM. For the sake of convenience, we introduce $w^=\u2212w$ and write the minimization scheme in terms of $w^$. For a given applied three-point force

*F*(Fig. 1), the potential energy of the system is Φ =

*U*−

*Fd*, where

*U*is the strain energy of the beam and $d=w^(x=L/2)$. Recalling that we are assuming the reinforcement plies to be shear rigid, only longitudinal strains contribute to the strain energy in the reinforcement. Conversely, the strain energy of the matrix plies includes contributions from both shear and longitudinal strains. The strain energy of the reinforcement plies is given by

*n*

_{r}is the number of the reinforcement ply and thus takes even values in the range −(

*N*− 1) ≤

*n*

_{r}≤ (

*N*− 1), while

*A*

_{r}≡

*bt*

_{r}. On the other hand, the strain energy of the matrix plies is

*n*

_{m}is the number of the matrix ply and thus takes odd values in the range −(

*N*− 1) ≤

*n*

_{r}≤ (

*N*− 1), while

*A*

_{m}≡

*bt*

_{m}. The total strain energy is then

*U*=

*U*

_{r}+

*U*

_{m}.

*Q*

_{j}is given by (16). Substituting $\Delta u\xafj*$ then reduces the expression for Φ to