## 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/MoS2 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 103 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 [35]. 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.

Fig. 1
Fig. 1
Close modal

## 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 Ei, Poisson ratio νi, and thickness ti, 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 HNtr + (N − 1)tm. The beam of overall length l = 1.2L 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¯≡Em/Er$, the beam slenderness ratio $L¯≡L/H$, the ply thickness ratio $t¯≡tm/tr$ and the number of reinforcement plies N, such that the reinforcement volume fraction $ϕr=[1+t¯(1−1/N)]−1$. 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 er.

### 2.3 Experimental Study

#### 2.3.1 Materials.

The reinforcement plies were sectioned from polycrystalline silicon wafers (PV Crystalox, Abingdon, UK), measuring tr = 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 tm = 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 Er = 160 ± 5 GPa, failure strength sr = 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 ersr/Er = 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 $∼5%$ strain, with failure strains of $>40%$. Table 1 provides the elastic properties for all of the materials used in the study.

Fig. 2
Fig. 2
Close modal
Table 1

Ply thickness ti and Young's modulus Ei, where the subscript “i” takes the value of “r” and “m” for the silicon reinforcement and polymer matrices, respectively

ConstituentMaterialti (mm)Ei (GPa)$E¯=Em/Er$$νi$
ReinforcementSilicon0.231600.22
MatrixLDPE0.050.138.1 × 10−40.45
TF-3100.050.053.1 × 10−40.45
TFL-2EA0.050.021.5 × 10−40.45
ConstituentMaterialti (mm)Ei (GPa)$E¯=Em/Er$$νi$
ReinforcementSilicon0.231600.22
MatrixLDPE0.050.138.1 × 10−40.45
TF-3100.050.053.1 × 10−40.45
TFL-2EA0.050.021.5 × 10−40.45

Note: The modulus ratio is defined as $E¯≡Em/Er$ with the Poisson ratio νi values taken from Ref. [17].

#### 2.3.2 Laminated Beam Fabrication.

We used a CO2 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.

Table 2

Sample dimensions and non-dimensional geometric groups for the composite beam samples tested in this study

MatrixL (mm)H (mm)$E¯$N$t¯$$L¯$
LDPE501.358.1 × 10−450.2137
TF-310501.353.1 × 10−450.2137
TFL-2EA501.351.5 × 10−450.2137
TFL-2EA261.351.5 × 10−450.2119
TFL-2EA121.351.5 × 10−450.219
TF-310503.043.1 × 10−4110.2117
TF-310504.723.1 × 10−4170.2111
MatrixL (mm)H (mm)$E¯$N$t¯$$L¯$
LDPE501.358.1 × 10−450.2137
TF-310501.353.1 × 10−450.2137
TFL-2EA501.351.5 × 10−450.2137
TFL-2EA261.351.5 × 10−450.2119
TFL-2EA121.351.5 × 10−450.219
TF-310503.043.1 × 10−4110.2117
TF-310504.723.1 × 10−4170.2111

#### 2.3.3 Test Methodology.

The beams were mounted for three-point bending (as illustrated in Fig. 1) onto steel rollers, with radii RL/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 σ/sr for a beam with N = 5 reinforcement plies, a $t¯=0.2$ ply thickness ratio, a $L¯=20$ slenderness ratio, and a modulus ratio $E¯=10−2$. 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 sr = Erer, where er = 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¯=10−4$ 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.

Fig. 3
Fig. 3
Close modal

A consequence of this change in behavior with decreasing $E¯$ was that the location of failure changes with $E¯$. For the $E¯=10−2$ 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 ɛ/er = 1). This is the expected failure mode for homogenous brittle materials in three-point bending. By decreasing the modulus ratio to $E¯=10−4$, 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¯=10−6)$, first failure was restricted to the top ply (Fig. 3(c)). Hence, decreasing $E¯$ promoted a change in failure mode from reinforcement rupture of the rear ply to the top ply.

These FE results can be used to infer the effect of $E¯$ on the flexural rigidity D and the first ply failure load P of the beams. The flexural rigidity D is defined as
$D≡FL348d$
(1)
so that the rigidity of the composite beam equals that of a homogeneous beam of span L and stiffness F/d. Figure 4(a) plots the dependence of the normalized flexural rigidity $D¯≡D/(NDr)$, where $Dr≡Erbtr3/12$ is the flexural rigidity of an individual reinforcement ply, as a function of $E¯$. The flexural rigidity of the previously analyzed N = 5, $t¯=0.2$, and $L¯=20$ beam remained nearly constant as $E¯$ decreased from 100 to 10−2. From $E¯=10−2$ to 10−6, $D¯$ decreased by more than an order of magnitude, with $D¯$ nearing unity at $E¯=10−6$. A value of $D¯=1$ implies that each reinforcement ply bends independently about its own neutral axis and is decoupled from the other plies. Intriguingly, values of $D¯<1$ were also observed and this is due to compliance of the matrix in the through-thickness direction of the beam. With decreasing $L¯$, the $D¯$ versus $E¯$ response shifts such that for a given $E¯$, the normalized rigidity $D¯$ decreases (Fig. 4(a)).
Fig. 4
Fig. 4
Close modal

The normalized first ply failure load $P¯≡P/(NPr)$ versus $E¯$ 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¯=0.2$, and $L¯=20$ beam, $P¯$ decreased with decreasing $E¯$ similar to $D¯$ (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¯$ with increasing $L¯$. In summary, the FE results indicate that the flexural rigidity, failure mode, and failure load depend upon both $L¯$ and $E¯$.

### 3.2 Experimental Measurements and Observations.

The measured load per unit width F/b versus deflection d responses for N = 5 beams with $t¯=0.21$ are plotted in Fig. 5(a). The modulus ratio $E¯$ was controlled by changing the polymer and the slenderness ratio $L¯$ via the span length L. By contrast, Fig. 5(b) plots the F/b versus d responses, with $E¯$ and $t¯$ fixed at 3 × 10−4 and 0.21, respectively, and the slenderness ratio $L¯$ 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¯$ and decreasing $L¯$ as expected from the FE simulations.

Fig. 5
Fig. 5
Close modal

Ultimate failure was catastrophic in all cases in Fig. 5. For instance, for the $E¯=3×10−4$ and $L¯=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¯=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¯=10−4$, N = 5, and $L¯=9$ beam.

Fig. 6
Fig. 6
Close modal

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.

In Euler–Bernoulli beam theory, the beam is assumed to be shear rigid with the beam rotation θ 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 wT of a Timoshenko beam is therefore given by the summation of the bending and shear contributions for a beam in three-point bending
$wT=−Fx2[(3L2−4x2)24DEB+1KAG]$
(2)
where DEB 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, AbH 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 DEB of the composite beam. The Ruess rules of mixtures give the shear modulus as
$G=[ϕrGr+1−ϕrGm]−1$
(3)
where GiEi/[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
$DEB=NDr[1+(1−1N)(E¯t¯3+(1+t¯)2η)]$
(4)
where $η=N[N+1+(N−2)t¯E¯]$ and $Dr=Erbtr3/12$. The effective flexural rigidity of the Timoshenko beam then follows from (1) as
$DT=DEB(1+12DEBKAGL2)−1$
(5)
which can be written in non-dimensional form as
$D¯T≡DTNDr=D¯EB{1+D¯EBϕr2(1−1/N)α×[1+E¯t¯(1−1/N)(1+νr)(1+νm)]}−1$
(6)
where
$α=KE¯2(1+νm)1t¯(L¯Nϕr)2$
(7)
and $D¯EB≡DEB/(NDr)$. In order to calculate the failure load, recall that the longitudinal strain ɛ and moment M are related to the curvature κ/dx via ɛ = −(/dx)z = −(M/DEB)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/(8D) at the bottom of the beam at mid-span. The first ply failure load PT (where the subscript “T” signifies prediction of the Timoshenko beam theory) is obtained by setting ɛmax = er, and we write the normalized failure load as
$P¯T≡PTNPr=D¯EBϕrN$
(8)
Fig. 7
Fig. 7
Close modal

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¯$ matched for large values of $E¯$, the TM under predicted $D¯$ for small $E¯$ (Fig. 4(a)). This was a consequence of the method used to homogenize the shear modulus which resulted in the incorrect limit of $D¯T→0$ as $E¯→0$. 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¯$ and shows no dependence on $L¯$ 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¯→1$ or as a stack of independent reinforcement plies when $E¯→0$. 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.

Here we present an approximate solution for the global bending of the composite beam using a Raleigh–Ritz method, with reinforcement and matrix plies considered to be shear rigid and deformable, respectively. The longitudinal strain within each ply was given by the superposition of a ply bending strain due to bending of the ply about its own neutral axis and a longitudinal strain resulting from the global bending of the beam. We thus write the strain within the nth ply at a location $z~$ measured with respect to the local neutral axis of each ply (Fig. 1) as
$εi=−dθidxz~+du¯(n)dx$
(9)
Here, as before, the subscript i takes values r or m for a reinforcement or matrix ply, respectively, and $u¯(n)$ is the longitudinal displacement field within the nth 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
$0=trθr+tmθm+u¯(n+1)−u¯(n−1)$
(10)
which reduces to
$γm=−(1+1t¯)dwdx−1tm(u¯(n+1)−u¯(n−1))$
(11)
In order to reduce the kinematic variables in the problem, we finally assume a linear variation of $u¯(n)$ from one reinforcement ply to the next such that $u¯(n)=nΔu¯/2$. (This assumption will be shown to be justified by comparing predictions of this model with the FE simulations.) The problem has thus been reduced and cast in terms of two field variables $Δu¯(x)$ and 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 $Δu¯(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
$D¯g≡FL348d1NDr=[96π4∑j=1∞sin2(jπ/2)j4D¯j]−1$
(12)
where
$D¯j=1+(1−1N)[E¯t¯3+(1+t¯)2η1+((jπ)2η/12α)]$
(13)
Numerical experimentation revealed that the first term of the series in (12) approximates the infinite series to within a few percent and hence it suffices to write $D¯g$ as
$D¯g≈1+(1−1N)[E¯t¯3+(1+t¯)2η1+(π2η/12α)]$
(14)
It is worth noting that the difference between (14) and the Euler–Bernoulli model (4) is the factor of [1 + (π2η)/(12α)]−1 in the last term. Similarly, the longitudinal strain at any position within the laminate can be expressed as
$εx(x,z)=−FLtrNDrπ2∑j=1∞[2z~tr+nQj]sin(jπ/2)sin(jπx/L)j2D¯j$
(15)
with
$Qj=12(1+t¯)/(jπ)2+E¯t¯/α12/(jπ)2+η/α+E¯t¯/α$
(16)
Since this analysis only accounts for global bending of the beam, the maximum tensile strain is always located at the bottom of the beam at its mid-span (i.e., x = L/2 and z = −H/2), leading to a normalized first ply failure load
$P¯g≡PgNPr=[8π2∑j=1∞(1+(N−1)Qj)sin2(jπ/2)j2D¯j]−1$
(17)

#### 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¯$. 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.

First consider the top ply resting on a foundation comprising the rest of the composite beam (Fig. 7(b)). The effective modulus of this foundation is given by the Reuss rule of mixtures as
$E¯f≡Ef/Er=1+t¯1+t¯/E¯$
(18)
and Biot [18] then gives the maximum tensile strain ɛl in the top reinforcement ply subjected to a transverse load F as
$εl=Ftr264/331/2DrE¯f1/3$
(19)
The total strain in the top reinforcement ply is the sum of the strain (19) due to local bending and the strain (15) due to global bending. The top ply failure load then follows as
$P¯t≡PtNPr=[433/2ϕrL¯(16E¯f)1/3+8π2∑j=1∞sin2(jπ/2)j2D¯j(1−(N−1)Qj)]−1$
(20)
Recalling that we assume that the reinforcement ply can fail only under tensile strains, the normalized failure load accounting for both local and global bending is
$P¯g+l={P¯gforP¯t≤0min(P¯t,P¯g)forP¯t>0$
(21)
so that top ply failure is operative when $P¯t>0$ and $P¯t.
The flexural rigidity due to combined global and local bending requires an expression for the top ply deflection due to both deformation modes. In a two-dimensional plane stress (or plane strain) setting, line loading of a half-space results in an unbounded displacement of the loading point. However, Biot [18] reported a calibration scheme which enables the use of the Winkler solution to provide a suitable approximation of the maximum displacement. The calibration involves using the continuum solution derived in Ref. [18] to set the Winkler foundation stiffness such that the maximum moment of the top ply in the Winkler and continuum solutions match. Following this procedure, the displacement of the loading point for a reinforcement ply on a foundation comprising the remainder of the composite beam (Fig. 7(b)) is given as
$dl=4311/2Ftr3DrE¯f$
(22)
The total displacement of the loading point is then the sum of the displacement due to this local bending and the global bending solution discussed in Sec. 4.2.1. The normalized effective flexural rigidity ratio accounting for both these contributions therefore follows as
$D¯g+l≡FL348d1NDr=[1D¯g+6439/2ϕr3L¯3N2E¯f]−1$
(23)

#### 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¯$ and $L¯$. 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¯$, $t¯$, and $E¯$ 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¯$.2 This map is interpreted as follows: Consider the case of $t¯=1$ and $E¯=10−4$, which is represented by the dashed-dotted line in the center (green-shaded) region. All designs (i.e., combinations of N and $L¯$) for this choice of $t¯$ and $E¯$ 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¯$ has a relatively small effect on the failure mode, while decreasing $E¯$ caused an appreciable shift whereby top ply failure is dominant over a larger design space.

Fig. 8
Fig. 8
Close modal

The design map in Fig. 8 suggests that in general decreasing $E¯$ 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¯$ 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

2

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≈L¯N(1+t¯)$, and thus, the minimum value of L/tr ≈ 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

A Rayleigh–Ritz minimization scheme was used to approximate the field variables w(x) and $Δu¯(x)$ in the global bending solution of the DPM. For the sake of convenience, we introduce $w^=−w$ 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 Φ = UFd, 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
$Ur=12∫0L[NDr(d2w^dx2)2+ErAr∑nr(nr2dΔu¯dx)2]dx$
(A1)
where nr is the number of the reinforcement ply and thus takes even values in the range −(N − 1) ≤ nr ≤ (N − 1), while Arbtr. On the other hand, the strain energy of the matrix plies is
$Um=12∫0L[EmAm∑nm(nm2dΔu¯dx)2+(N−1)Dm(−1t¯d2w^dx2+1tmdΔu¯dx)2+(N−1)KGmAm{−(1+1t¯)dw^dx+1tmΔu¯}2]dx$
(A2)
where similarly nm is the number of the matrix ply and thus takes odd values in the range −(N − 1) ≤ nr ≤ (N − 1), while Ambtm. The total strain energy is then U = Ur + Um.
In order to employ the Rayleigh–Ritz scheme, we approximate $w^(x)$ and $Δu¯(x)$ via the Fourier series
$w^(x)=∑j=1∞w^jsin(jπxL)$
(A3)
and
$Δu¯(x)=∑j=1∞Δu¯jcos(jπxL)$
(A4)
respectively, and thus reducing the problem to determine the coefficients $w^j$ and $Δu¯j$ via the Rayleigh–Ritz minimization. The expressions (A3) and (A4) are chosen so that the kinematic conditions $w^(x=0)=w^(x=L)=0$ are automatically satisfied, while, without loss of generality, we can set $Δu¯(x=L/2)=0$. The optimum values of $w^j$ and $Δu¯j$ are then calculated by substituting (A3) and (A4) into the expression for Φ and minimizing with respect to $w^j$ and $Δu¯j$.
Setting $∂Φ/∂(Δu¯j)=0$ gives the optimum value of $Δu¯j$ denoted by $Δu¯j*$ as
$Δu¯j*=w^jtrjπLQj$
(A5)
where Qj is given by (16). Substituting $Δu¯j*$ then reduces the expression for Φ to
$Φ=NDrπ44L3∑j=1∞j4w^j2D¯j−F∑j=1∞w^jsin(jπ/2)$
(A6)
where $D¯j$ is given by (13). Now setting $∂Φ/∂w^j=0$ provides the optimum values of $w^j$ denoted by $w^j*$ as
$w^j*=2FL3sin(jπ2)NDrD¯jj4π4$
(A7)

We can thus calculate the deflection of the beam at x = L/2, and the flexural rigidity (12) then follows. Furthermore, substituting the solution for w and $Δu¯$ into (9) gives the longitudinal ply strain (15).

## References

1.
Wegst
,
U. G. K.
,
Bai
,
H.
,
Saiz
,
E.
,
Tomsia
,
A. P.
, and
Ritchie
,
R. O.
,
2015
, “
Bioinspired Structural Materials
,”
Nat. Mater.
,
14
(
1
), pp.
23
36
. 10.1038/nmat4089
2.
Wang
,
J.
,
Cheng
,
Q.
, and
Tang
,
Z.
,
2012
, “
Layered Nanocomposites Inspired by the Structure and Mechanical Properties of Nacre
,”
Chem. Soc. Rev.
,
41
(
3
), pp.
1111
1129
. 10.1039/C1CS15106A
3.
Begley
,
M. R.
,
Philips
,
N. R.
,
Compton
,
B. G.
,
Wilbrink
,
D. V.
,
Ritchie
,
R. O.
, and
Utz
,
M.
,
2012
, “
Micromechanical Models to Guide the Development of Synthetic ‘Brick and Mortar’ Composites
,”
J. Mech. Phys. Solids
,
60
(
8
), pp.
1545
1560
. 10.1016/j.jmps.2012.03.002
4.
Sakhavand
,
N.
, and
Shahsavari
,
R.
,
2015
, “
Universal Composition–Structure–Property Maps for Natural and Biomimetic Platelet–Matrix Composites and Stacked Heterostructures
,”
Nat. Commun.
,
6
(
1
), p.
6523
. 10.1038/ncomms7523
5.
Kim
,
Y.
,
Kim
,
Y.
,
Lee
,
T.-I.
,
Kim
,
T.-S.
, and
Ryu
,
S.
,
2018
, “
An Extended Analytic Model for the Elastic Properties of Platelet-Staggered Composites and Its Application to 3D Printed Structures
,”
Compos. Struct.
,
189
, pp.
27
36
. 10.1016/j.compstruct.2018.01.038
6.
Ghugal
,
Y. M.
, and
Shimpi
,
R. P.
,
2001
, “
A Review of Refined Shear Deformation Theories for Isotropic and Anisotropic Laminated Beams
,”
J. Reinf. Plast. Compos.
,
20
(
3
), pp.
255
272
. 10.1177/073168401772678283
7.
Moussiaux
,
E.
,
Brinson
,
H. F.
, and
Cardon
,
A. H.
,
1987
,
Bending of a Bonded Beam as a Test Method for Adhesive Properties, ADA184813
,
Virginia Tech Center for Adhesion Science
,
Blacksburg, VA
.
8.
,
R. D.
, and
Weinstein
,
A. S.
,
1975
, “
Flexural Stiffness of Sandwich Beams
,”
J. Eng. Mater. Technol.
,
97
(
3
), pp.
264
270
. 10.1115/1.3443294
9.
Nguyen
,
Q.-H.
,
Martinelli
,
E.
, and
Hjiaj
,
M.
,
2011
, “
Derivation of the Exact Stiffness Matrix for a Two-Layer Timoshenko Beam Element With Partial Interaction
,”
Eng. Struct.
,
33
(
2
), pp.
298
307
. 10.1016/j.engstruct.2010.10.006
10.
Škec
,
L.
, and
Jelenić
,
G.
,
2013
, “
Analysis of a Geometrically Exact Multi-Layer Beam With a Rigid Interlayer Connection
,”
Acta Mech.
,
225
(
2
), pp.
523
541
. 10.1007/s00707-013-0972-5
11.
Ranzi
,
G.
,
2008
, “
Locking Problems in the Partial Interaction Analysis of Multi-Layered Composite Beams
,”
Eng. Struct.
,
30
(
10
), pp.
2900
2911
. 10.1016/j.engstruct.2008.04.006
12.
Sousa
,
J. B. M.
, Jr.
, and
da Silva
,
A. R.
,
2010
, “
Analytical and Numerical Analysis of Multilayered Beams With Interlayer Slip
,”
Eng. Struct.
,
32
(
6
), pp.
1671
1680
. 10.1016/j.engstruct.2010.02.015
13.
Averill
,
R. C.
,
1994
, “
Static and Dynamic Response of Moderately Thick Laminated Beams With Damage
,”
Compos. Eng.
,
4
(
4
), pp.
381
395
. 10.1016/S0961-9526(09)80013-0
14.
Tessler
,
A.
,
Di Sciuva
,
M.
, and
Gherlone
,
M.
,
2009
, “
A Refined Zigzag Beam Theory for Composite and Sandwich Beams
,”
J. Compos. Mater.
,
43
(
9
), pp.
1051
1081
. 10.1177/0021998308097730
15.
Bathe
,
M.
,
Heussinger
,
C.
,
Claessens
,
M. M. A. E.
,
Bausch
,
A. R.
, and
Frey
,
E.
,
2008
, “
Cytoskeletal Bundle Mechanics
,”
Biophys. J.
,
94
(
8
), pp.
2955
2964
. 10.1529/biophysj.107.119743
16.
ASTM D882, 2010
,
2010
,
Standard Test Method for Tensile Properties of Thin Plastic Sheeting
,
ASTM International
,
West Conshohocken, PA
.
17.
CES EduPack
,
2015
,
Granta Design Ltd., Cambridge, UK
.
18.
Biot
,
M. A.
,
1937
, “
Bending of an Infinite Beam on an Elastic Foundation
,”
ASME J. Appl. Mech.
,
4
, pp.
A1
A7
.