# Analytical Prediction and Numerical Verification of Stress Concentration
Profile Around an In-situ Tow Break in Resin-impregnated Filament-wound
Composites

Jiakun Liu<EMAIL_ADDRESS>Stuart Leigh Phoenix<EMAIL_ADDRESS>Sibley
School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY
14853, USA

###### Abstract

A new empirical analytical approach is developed for predicting the stress
concentration profile around an in-situ tow break in filament-wound
composites. A shear-lag analysis is firstly performed to solve for the
perturbational axial displacement of the broken tow and resultant debonding
lengths. Solution of stress field caused by tangential load on the surface of
a elastic half space is then utilized in combination with superposition
concepts to obtain the overload magnitudes in the neighboring tows.
Subsequently, high-fidelity finite element analysis on a representative uni-
directional laminate model under different stress states is performed, and
excellent overall agreement is observed between analytical and numerical
predictions. The proposed method takes into account essential aspects such as
transversely isotropic material properties, in-situ stress states and their
effect on the interfacial frictional forces in the debonded interfaces, and
thus provides a convenient way to evaluate the stress concentration factors in
damaged filament-wound composites. In addition, this approach can be applied
to yield auxiliary failure evaluation criteria for statistical strength
prediction or finite element modeling of filament-wound composites or similar
structures.

###### keywords:

Filament-wound Composites , Resin-impregnated Tows , Stress Concentration
Factor , Finite Element Modeling

soft open fences, separator symbol=; ,

111Current affiliation: Department of Mechanical Engineering and Applied
Mechanics, University of Pennsylvania, 220 S 33rd St, Philadelphia, PA 19104,
USA

##  1 Introduction

###  1.1 Challenges in damage analysis of FWCs

The filament winding technique is efficient, cost-effective, and has been
applied in fields including civil and environmental, aerospace, and automotive
industries. It is primary used to manufacture hollow or convex-shaped
components such as wind turbines, pipes, gas tanks, and composite overwrapped
pressure vessels (COPVs). Fig.1 illustrates basic procedures of filament
winding, during which combed fiber strands are impregnated in hot resin bath,
and then the resin-impregnated tows are wound onto a rotating mandrel to
fabricate the filament-wound composites (FWCs). Filament tows are usually
impregnated right before being wound onto the mandrel (wet-winding), but can
also be priorly-impregnated (dry-winding), or post-impregnated. During
winding, a moving guide (pay-out eye) wraps filament bands onto a rotating
mandrel following a layer-over-layer routine. The motion and path of the guide
are controlled by computer-aided winding machines to help realize desired
layer angle at specific positions and targeted thickness of each winding
layer. After completion of winding, the entire structure is placed in an oven
or under radiant heaters for curing. Thereafter, the mold core may be removed
or left as an integral component of the final product.

Figure 1: Basic procedures for resin-impregnation of carbon fiber tows and
filament winding process (this plot shows wet-winding).

Since FWC structures are widely used in engineering applications, and their
failure could lead to devastating losses (e.g., burst of COPVs), it is of
critical importance to investigate the damage and load redistribution
mechanisms in FWCs. However, these remain as challenging tasks due to the
complicated material and structural characteristics [1, 2, 3].

The first challenge in studying failure mechanisms of composites is associated
with the length scales. Composite materials typically exhibit a strong
dependence on varying length scales, commonly categorized into: micro, meso,
and component scales. Micro scale is generally in the magnitude of microns due
to the typical fiber diameters (e.g., $3$ to $10\mu\text{m}$). Meso scale
usually ranges from several to hundreds of millimeters depending on the object
of interest, and a typical meso-scopic representative volume element (RVE) is
used to show the geometrical features. For instance, a periodic RVE that
demonstrate the morphology of woven textile composites. Component scale is on
the magnitude of a structure or part (e.g., composite wing blades) for
structural analysis. This multi-scale analysis concept is important because
certain features may only be analyzed or modeled at a specific length scale.
Also, due to the complicated geometries, material heterogeneity, and enormous
number of fibers in composite structures, direct linking and information
sharing between scales is impractical. As a result, computational
homogenization and multi-scale modeling methods are widely applied for
analysis of composites.

In addition to length scale, another challenge is the morphological complexity
of FWCs. Firstly, the local pattern of filament can be highly intricate. For
example, overwrap of a cylindrical COPV usually has both both helical and hoop
layers, each with specific winding angles and layer thicknesses, and
henceforth, the stacking angles of tows vary among different locations and
also between distinct layers. Secondly, the material integrity and
distribution are not homogeneous. For example, raw roving strands (used as
source material of filaments) can be regarded as carbon fiber tows embedded in
dry matrix, and as these tows pass through hot resin flows, additional wet
resin coatings are adhered to the strand surfaces (shown schematically in
Fig.2), resulting in resin-rich areas among tows and between different layers.

In a mechanical perspective, unevenly distributed fiber volumes lead to ’weak
interfaces’ affecting local stress redistribution and damage paths, which, in
some cases, could dominate the strength and failure mechanisms of composites.
Such kinds of heterogeneity and resin-rich areas can be evidently observed in
both intact and damaged filament-wound and/or prepreg-fabricated composites
[4, 5, 6]. Meanwhile, discrete failure patterns in form of split strands/tows
can be observed in widely available photos of burst/destructed FWCs [7, 8, 9,
10, 11, 12, 13], where interfacial cracks among broken tows and ’brush-like’
damage patterns (split and protruding tows) are commonly observed failure
modes. These phenomena are expected considering the fact that initial raw
materials used for filament wounding process are in the formation of such
’tows’ or ’strands’.

###  1.2 Theoretical strength prediction of fibrous composites

To provide theoretical foundation for strength prediction of FWCs (e.g.,
COPVs), or more generally, fiber-reinforced polymeric composites, an extensive
amount of studies on stress redistribution mechanisms between (micro-scopic)
fibers and matrices have been presented in the past decades, such as Refs.
[14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], among many others. These
work typically focus on a fiber-matrix system with assumed load-sharing
mechanism among one or multiple fiber breaks in specific geometrical patterns.
Through which, the stress concentration factor (SCF) regarding a cluster of
fiber breaks (at micro scale) can be quantified, and thereafter, the strength
of large volume composites can be predicted by finding a ’critical cluster
size’ (number of broken fibers) above which the SCF would be large enough to
trigger catastrophic failure. This concept is somehow comparable to the
’critical crack length’ in classical liner elastic fracture mechanics (LEFM)
for strength prediction of brittle materials, above which the resultant stress
intensity factor would become significant enough to trigger unstable crack
propagation.

When applying the above micro-mechanical fiber breakage analysis results to
predict strength of component-scale composites with fiber volume $V$ (here
referring to the number of fibers), the entire composite is typically
conceptualized as an ’aggregate’, or a ’linkage’, of $n$ inter-connected
micro-scale fiber-matrix bundles, where each bundle contains $m$ fibers so
that $V=mn$. Then the statistical strength prediction is performed by finding
the critical number of fiber breaks (and distributions) and the corresponding
probabilities of failure of such a linked fiber bundles [26, 27, 28]. To apply
micro-mechanical results, bundle size has to be limited (e.g., $m<20$) so that
load sharing mechanisms correlated to all possible fiber breakage amount and
patterns in a bundle, and henceforth SCFs in unbroken fibers in the bundle,
can be analytically quantified. And all $n$ bundles are subsequently assumed
to share similar load-sharing mechanisms. However, in reality, load
redistribution among broken and intact fibers may vary among different
locations, and may also depend on multiple factors such as winding patterns,
in-situ stress states. For instance, in a cross-ply ($\pm\theta^{\degree}$)
laminate, assume a cluster of fiber breaks occurred around the top of an inner
$+\theta^{\degree}$ ply below a $-\theta^{\degree}$ ply, then the load
redistribution mechanisms between fibers in the same ply would be different
from the mechanism between fibers across the cross plies.

###  1.3 Numerical approaches for failure analysis of FWCs

Other than analytical methods and experimental studies, finite element
analysis (FEA) is nowadays widely applied for damage analysis and strength
prediction of FWC structures [7, 29, 8, 30, 31, 32, 33, 34, 9], commonly
directly at component scale (e.g., a numerical model of a COPV). In which, the
composite part is typically treated as a ’smeared’ continuum laminate
consisted of many winding plies, and the effective properties of such a
laminate is obtained by classical laminate theory . Then conventional
continuum damage mechanics (CDM) based failure criteria (e.g., max strain,
Tsai-Hill, Hashin) or proposed failure evaluation methods with more elaborate
algorithms are incorporated with these elements as an ’indicator’ of their
damage levels. Henceforth the burst strength prediction of FWCs depends on the
evaluation of selected failure criteria. However, as mentioned earlier,
composites typically exhibit strong scale dependency and characteristics,
whereas CDM-based criteria are inherently homogenized evaluations on continuum
solid elements with ’smeared’ laminate properties. Therefore, despite
component-scale structure and loading features (e.g. COPV shape, effective
layer-level stress states) are inherently retained in these kinds of studies,
advances in understanding of micro- and meso-scale damage mechanisms would be
beneficial in terms of refining in CDM-based criteria, or serving as a failure
evaluation method for multi-scale damage analysis.

###  1.4 Necessity of damage analysis of a broken tow at meso scale

In general, the above discussion indicates that a study that covers essential
damage mechanisms around a broken tow and also realizes accurate prediction of
resultant overload magnitudes is strongly needed. Such a study, on the one
side, could provide extra conceptual foundations for statistical strength
prediction in analytical approaches that commonly merely rely on micro-
mechanical results, and on the other side, could be applied to enhance, or
integrate with, existing failure criteria used for component-scale numerical
models.

Under such a motivation, in this study we develop a new analytical approach to
quantify the overload magnitude around a tow break inside a FWC structure, and
then compare analytical results to numerical predictions realized by high-
fidelity finite element simulations. The remaining of this article is
organized as follows: In Section 2, we start with a shear-lag based analysis
that predicts axial deformation of broken tow and interfacial debonding
lengths. Then in Section 3, we derive algorithms based on classical contact
mechanics (point load analysis) that quantitatively evaluates stress profiles
around a broken tow. Subsequently, in Section 4, a series of finite element
simulations are performed and compared with analytical predictions. Finally, a
discussion with concluding remarks are given in Section 5.

##  2 Shear-lag based analysis of an in-situ tow break

Many studies involving SCF of fiber breaks applied so-called shear-lag
analysis, wherein fibers are considered to be relatively rigid and hence carry
all the tensile load, while the epoxy matrix interfaces are the ’load
redistributor’ that only carry and deliver shear stresses. Earlier work
typically assume perfectly-bonded fiber/matrix interfaces even when a fiber
breaks [14, 15]. Some later studies included debonding into consideration,
wherein the interfacial behavior is commonly described in a ’debond-then-slip’
manner [17, 20, 21, 35, 23].

In this study, a similar shear-lag concept is used. In specific, it is assumed
that the tow breakage, resulted from a critical cluster of fiber breaks inside
the tow, will also cause debonding along tow interfaces. Then as the broken
tow slides against the surrounding neighbors, interfacial cracks keep
propagating until a local equilibrium state has been reached. Such an
equilibrium state is when the accumulated shear force, provided by frictional
forces in debonded interfaces and elastic shear in deformed interfaces, are
sufficient to ’stop’ the broken tow and fully recovers its stress after a
certain distance. In some conditions, a tow break may become unstable, for
example, a broken tow on the outer surface of an COPV may just pump out of the
overwrap due to insufficient boundary constraints. In other situations, the
overload exerted by tow break(s) may become high enough to trigger continued
tow breaks thereby the entire structure will fail catastrophically in a
transient amount of time, essentially a typical failure mode (burst) of FWCs.

###  2.1 Representative tow cross section

A schematic plot of tow cross-section and resin-rich zones is shown in Fig.2.
Unlike in the micro-mechanical analysis where fiber and matrix have comparable
dimensions, the resin-rich zone between tows are considered to be much smaller
than the cross-sectional dimension of composite tows. The shape of the cross-
section is set to be ’almost rectangular’ with height $h$ and width $w$, while
in practice it could be somewhere between oval and rectangular.

Let $t_{\Leftrightarrow}$ and $t_{\Updownarrow}$ be the interface thickness
between neighboring tows and adjacent layers, respectively, where the
directional arrows used as subscripts are intended to described the interface
relations: the left-right arrow $\Leftrightarrow$ indicating the intra-ply
property/interaction between tows in the same ply, while the up-down arrow
$\Updownarrow$ denotes the inter-ply property/interaction across different
plies. These subscripts will be used consistently for several other intra- and
inter-ply variables in the following contents.

Figure 2: Schematic plot of material formation in a resin-impregnated tow and
dimensions of its representative cross section.

In general, the aspect ratio of tow cross section ranges from $1.0<w/h<3.0$,
and the cross-sectional area $w\cdot h$ is typically less than
$1\,\text{mm}^{2}$. The number of fibers within the cross-sectional area of a
tow generally ranges from 800 to 2400. For instance, a high modulus carbon
fiber strands, Toray® T1100G, is available in 12k and 24k (density of fibers
in tow), based on the fiber diameter ($5\mu\text{m}$) and its typical volume
fraction ($60\%$), the cross-sectional area of 12k tow can thus be
approximated to be around $0.393\text{mm}^{2}$.

###  2.2 Coordinate system and notations involving the tow break

Consider a small ’piece’ cut out of a filament-wound composites, it can then
be characterized as a laminate consisted of plies whose stacking layups are
equivalent to the winding angles at the cut location. Use x-y-z to denote the
global Cartesian coordinate system (CSYS), and use $1{-}2{-}3$ to refer to the
local tow material CSYS of each ply, where ’$1$’ is the longitudinal fiber
orientation of tows, ’$2$’ the transverse direction (vertical to ’$1$’), and
’$3$’ the through-thickness orientation, respectively. Let the winding
direction (longitudinal tow orientation ’1’) of the ply containing the broken
tow aligned with x-axis of the global system (Fig.3). The winding angle of the
$i\text{th}$ ply, $\theta^{i}\in\interval{-90^{\circ}}{90^{\circ}}$ denotes
the angle between local fiber orientation of the ply to the global x-axis. Let
$z_{i}$ ($i=1,2...N$) be the coordinate of bottom surfaces of plies in the $z$
direction, and $z_{N+1}$ be those of the top surfaces of the $N$th ply. The
thickness of the $i$th ply is thus $h^{(i)}=z_{i+1}-z_{i}$. These conventions
will be used in the remaining contents, where subscripts $x,y,z$ refer to any
measure/value of a variable in the global coordinates, subscripts $1,2,3$
refer to measure/value in the local material orientation of a specific ply,
and superscript $(i)$ denote those of the $i\text{th}$ ply.

In this study, a primary interest is the stress profile in the neighboring
tows right next to the broken tow, since these regions/locations are expected
to suffer highest overload. Set the center of this ply as the origin $O$
(Fig.3), then following the previously defined superscript $(k)$, let $(k,j)$
be the $j\text{th}$ nearest tow (in the $k$th ply) next to the center tow
along the positive $x$ direction, and $(k,0)$ refers to the broken tow itself.
Similarly, as shown in Fig.3B, denote the tow regions in the above and below
plies in the same manner.

Figure 3: Coordinate system and notations for plies and tows. (A) Intra-ply
slice view ($z=0$ plane) around the tow break, showing notations of the
neighboring tows in the $k$th ply. (B) Through-thickness slice view in the
plane of tow break ($x=0$ plane), showing notations of the neighboring tow
regions in different plies. (C) If the adjacent plies have different winding
angles, the effective tow region may cover a portion of multiple tows. Here as
an example, assume the first ply above the broken tow has a winding angle of
$-30^{\circ}$, then superscript ($k+1,j$) corresponds to an effective region
covering multiple tows.

Now let the $k$th ply be an inner ply ($2\leq k\leq N-1$) with winding angle
$0^{\circ}$, assume that a tow break occurs in the center of this ply. The
dimensions of selected representative laminate are thought to be sufficiently
large to cover the stress recovery lengths, so that the stresses or strains at
the boundaries are unaffected by the broken tow. Assume that the intra-ply
stresses are uniform through the thickness of each ply, whereas the through-
thickness stresses may have a slight variation. Denote, respectively, the
’far-field’ intra-ply stresses and average through-thickness stress of the ply
as

$\sigma^{(k)}_{11,\infty}=\sigma^{(k)}_{xx,\infty},\quad\sigma^{(k)}_{22,\infty}=\sigma^{(k)}_{yy,\infty},\quad\bar{\sigma}^{(k)}_{33}=\bar{\sigma}^{(k)}_{zz}=\int_{z_{k}}^{z_{k+1}}\sigma^{(k)}_{33}dz$
(1)

where the equivalence between local ($1{-}2{-}3$) and global ($x{-}y{-}z$)
stresses is due to the assumed $0^{\circ}$ winding angle this ply. It is also
noted that the intra-ply shear stress $\sigma^{(k)}_{12}$ of this $0^{\circ}$
ply, under the bi-axial loading, is neglected in the following shear analysis.

Fig.4 shows a schematic plot of the broken tow and its intra-ply neighboring
tows. The resultant intra-ply debonding lengths, under the in-situ stress
states, are $L_{\Leftrightarrow}$ on each side. The distance from the origin
to the surfaces of broken tow is $\delta_{cod}$. Since the deformation and
debonding lengths of this ply is assumed to be symmetric about the plane of
tow break, here $\delta_{cod}$ is half of the real ’gap’ between the two
fracture surfaces.

Figure 4: Intra-ply cut-view of an in-situ tow break occurred at center of the
$k$th ply with winding angle $\theta^{k}=0^{\circ}$.

Fig.5 shows a schematic inter-ply cross-sectional cut-view of the broken tow
and its intra- and inter-ply neighboring tows, where the inter-ply debonding
lengths on the top and bottom interfaces of the broken tow are assumed to be
$L_{\Updownarrow}$. This is based on the assumption that $\sigma_{33}^{(k)}$
does not vary significantly through one ply in the through thickness
direction, so that the debonding lengths between the above and below ply
interfaces are comparable if these plies have the same winding angles.

Figure 5: Inter-ply cross-sectional cut-view of an in-situ tow break and
corresponding interface debonding.

Note that both Fig.5 and Fig.4 show a representative laminate as a stacking of
$0^{\degree}$ tows having the same ply thicknesses. In general, since winding
layers can have different angles and thicknesses, the region of interest in
the other layers are not necessarily in the form of individual tow/yarns.
Instead, such a region may cover multiple tows due to the mismatched winding
angles (as illustrated in Fig.3C showing a $-30^{\circ}$ layer).

###  2.3 Elastic properties and stress states

Material properties of composite tows are regarded here as linear elastic,
brittle, transversely isotropic. Thus the lamina/ply stress-strain relation,
in the corresponding local CSYS, is

$\underset{-}{\varepsilon^{(i)}_{local}}=\underset{=}{S_{tow}}:\underset{-}{\sigma^{(i)}_{local}}\quad\rightarrow\quad\underset{-}{\sigma^{(i)}_{local}}=\underset{=}{S_{tow}}^{-1}:\underset{-}{\varepsilon^{(i)}_{local}}=\underset{=}{C_{tow}}:\underset{-}{\varepsilon^{(i)}_{local}}$
(2)

where the single and double underlines are used to denote first and second
order tensors, $\varepsilon^{(i)}_{\,local}$ and $\sigma^{(i)}_{local}$ are
the local strain and stress vectors, respectively, ${S}_{tow}$ is the
compliance matrix of the tow, and its inverse $C_{tow}$ is the stiffness
matrix.

There are five independent variables to be defined: let $E_{l}$ be the Young’s
modulus of composite tows in the longitudinal (fiber) direction, $E_{t}$ be
the modulus in transverse directions; $\nu_{lt}$ and $\nu_{tt}$ be the
corresponding Poisson’s ratios; $G_{lt}$ be the longitudinal-to-transverse
shear modulus. Then Eq.2 takes the form

$\begin{bmatrix}\varepsilon_{11}\\\ \varepsilon_{22}\\\ \varepsilon_{33}\\\
\gamma_{23}\\\ \gamma_{13}\\\
\gamma_{12}\end{bmatrix}=\begin{bmatrix}1/E_{l}&-\nu_{tl}/E_{t}&-\nu_{tl}/E_{t}&0&0&0\\\
-\nu_{lt}/E_{l}&1/E_{t}&-\nu_{tt}/E_{t}&0&0&0\\\
-\nu_{lt}/E_{l}&-\nu_{tl}/E_{t}&1/E_{t}&0&0&0\\\ 0&0&0&1/G_{tt}&0&0\\\
0&0&0&0&1/G_{lt}&0\\\
0&0&0&0&0&1/G_{lt}\end{bmatrix}\begin{bmatrix}\sigma_{11}\\\ \sigma_{22}\\\
\sigma_{33}\\\ \sigma_{23}\\\ \sigma_{13}\\\ \sigma_{12}\end{bmatrix}$ (3)

where $\gamma_{23},\gamma_{13},\gamma_{12}$ are engineering shear strains, and

$G_{tt}=\frac{E_{t}}{2(1+\nu_{tt})},\qquad\nu_{tl}=\nu_{lt}\frac{E_{t}}{E_{l}}$
(4)

are the shear modulus in the transverse plane, and Poisson’s ratio in
transverse-to-longitudinal directions, respectively.

Each ply is considered to be in 3D stress states, among which, shear stresses
associated with through thickness directions ($\sigma_{23}$ and $\sigma_{13}$)
are neglected, namely, the ply-level effective stresses are

$\underset{-}{\sigma^{(i)}_{local}}=\left[\sigma^{(i)}_{11},\sigma^{(i)}_{22},\sigma^{(i)}_{33},0,0,\sigma^{(i)}_{12}\right]^{T}$
(5)

and such an assumption on the ply-level stress state is based on the following
considerations: 1) This is a ’piece’ of filament-wound laminate whose
boundaries are constrained by the surrounding materials. FWC structures are
typically axi-symmetric in geometry and the loading (pressure) applied on the
inner surface is also axi-symmetric. And henceforth the shear deformation
associated with through thickness directions ($\sigma_{23},\sigma_{13}$) are
considered to be negligible. 2) Even for a high performance COPV, pressure
load through the thickness direction is at least one order of magnitude lower
than the intra-ply stresses. Typical burst pressure of commercial COPVs are
less than $20,000\,$psi ($\approx$ 138MPa), which is significantly smaller
than the fiber tensile stresses. Therefore the magnitudes of
$\sigma_{23},\sigma_{13}$ are expected to be small as compared to fiber
tensile stresses. 3) This study involves in-situ analysis of damages in
pressure-bearing structures, so axial stress in the thickness direction,
$\sigma_{33}$, cannot be neglected.

###  2.4 Interfacial shear stress and force equilibrium

For composite prepreg tows, $E_{l}$ is typically one order of magnitude higher
than $E_{t}$ and the modulus of resin. Therefore, similar to the
aforementioned micro-mechanical shear-lag analysis, assume that the axial
strains and stresses are uniform through the cross section of the broken tow,
so that they only depends on $x$ (distance away from the plane of tow break).
Then based on the previously established notations, for the broken tow, the
axial stress along the local fiber orientation is $\sigma^{(k,0)}_{11}(x)$.
Let $\delta u^{(k,0)}_{11}(x)$ be the relative displacement between the broken
tow and its adjacent neighbors, namely, $\delta u^{(k,0)}_{11}(x)$ is the
extra, perturbation displacement caused by the tow break. Initially when there
is no broken tow, $\delta u^{(k,0)}_{11}(x)=0$, and all the tows have a common
axial strain and uniform longitudinal stress $\sigma^{(k)}_{11,\infty}$.

The relation between $\delta u^{(k,0)}_{11}$ and $\sigma^{(k,0)}_{11}$ is
thus:

$\sigma^{(k,0)}_{11}(x)=\dfrac{d}{dx}\left[\delta
u^{(k,0)}_{11}(x)\right]\cdot E_{l}+\sigma^{(k)}_{11,\infty}$ (6)

where the contribution on $\sigma_{11}$ from the other strain components are
neglected in the stress-strain relation. To account for the essential in-situ
effects, transverse and through-thickness stresses need to be included because
they are correlated with sliding/clamping constraints that will further affect
the debonding and stress recovery lengths.

Take a small differentiable tow element of length $dx$, the equilibrium of
forces acting on such an element requires that:

$\frac{d}{dx}\left[\sigma^{(k,0)}_{11}(x)\right]\cdot(wh)=2\cdot\uptau_{\Leftrightarrow}(x)\cdot
h+2\cdot\uptau_{\Updownarrow}(x)\cdot w$ (7)

where $\uptau_{\Leftrightarrow}$ and $\uptau_{\Updownarrow}$ are the intra-ply
and inter-ply shear tractions along the interfaces between the broken tow and
its surrounding neighbors, respectively.

Since the definition of $\delta u^{(k,0)}_{11}(x)$ is the perturbation
displacement of the broken tow, a straightforward way to define interfacial
shear strains around the broken tow is:

$\gamma^{(k,0)}_{\Leftrightarrow}(x)=\frac{\delta
u^{(k,0)}_{11}(x)}{t_{\Leftrightarrow}}\,,\qquad\gamma^{(k,0)}_{\Updownarrow}(x)=\frac{\delta
u^{(k,0)}_{11}(x)}{t_{\Updownarrow}}$ (8)

Let $G_{\Leftrightarrow}$ and $G_{\Updownarrow}$ be the shear modulus of
corresponding interfaces, whose magnitude should lie within that of the shear
modulus of pure resin or the transverse stiffness of the tows, $E_{t}$.
Considering that the thickness of interfaces are relatively small, thereby
define, respectively:

$K_{s,\Leftrightarrow}=\frac{G_{\Leftrightarrow}}{t_{\Leftrightarrow}}\;,\qquad
K_{s,\Updownarrow}=\frac{G_{\Updownarrow}}{t_{\Updownarrow}}$ (9)

as the effective shear stiffness of the intra-ply and inter-ply interfaces in
pure Mode-II or Mode-III shear deformations. This gives a direct relation
between the relative displacement and the interfacial shear stresses as:

$\displaystyle\uptau^{(k,0)}_{\Leftrightarrow}(x)$
$\displaystyle=\gamma_{\Leftrightarrow}(x)\cdot
G_{\Leftrightarrow}=\frac{\delta
u^{(k,0)}_{11}(x)}{t_{\Leftrightarrow}}=K_{s,\Leftrightarrow}\cdot\delta
u^{(k,0)}_{11}(x)$ (10) $\displaystyle\uptau^{(k,0)}_{\Updownarrow}(x)$
$\displaystyle=K_{s,\Updownarrow}\cdot\delta u^{(k,0)}_{11}(x)$ (11)

which has the similar form as the commonly used traction-separation rules in
cohesive zone modeling (CZM) techniques.

Sine the interfacial shear stresses depend on status of the interface (still
bonded or has failed), whereas the elastic traction-separation response
(Eq.10) is only valid for bonded domains, it remains to define the frictional
sliding forces in the debonded interfaces. Let $d_{\Leftrightarrow}(x)$,
$d_{\Updownarrow}(x)$ be the damage variables ($d\in\interval{0}{1}$) of
intra-ply and inter-ply interfaces, respectively, with $d=0$ indicating that
an interface is intact and $d=1$ means complete interface debond. This kind of
application of damage variable is also commonly used in CZM techniques in FEA
where the interface penalty stiffness of a cohesive element is typically
degraded by a multiplier depending on the damage variables (e.g., $1-d$) to
account for interfacial damage and failure status.

As shown in Fig.6, assume that the damage evolution of the interface is rapid,
such that the damage variable reaches one instantly when the failure criterion
has been reached. Namely, once the interfacial relative displacement has
reached a critical value, it breaks instantly. After debonding, the interface
may still carry frictional forces depending on the contact pressure between
the debonded interface. For pure shear modes (Mode-II and III), assume that a
frictional force may exist after debonding (Fig.6B), where $\tau^{f}_{s}$ and
$\tau^{f}_{n}$ are the interface strengths in shear and normal modes,
respectively, while $\delta^{f}_{s}$ and $\delta^{f}_{n}$ are the
corresponding separation to cause interfacial debond, respectively.

Figure 6: Schematic plot of interfacial behavior in pure shear modes (Mode-II,
Mode-III), where $d$ is the effective damage variable of the interface, $K$ is
the penalty stiffness, $\delta^{f},\tau^{f}$ are the failure (debonding)
separation and traction, respectively. (A) Typical traction-separation
response, where the interface debonding occurs upon reaching the failure
criterion. (B) If a debonded interfacial region is still in contact and under
negative normal strain, pressure-dependent friction may exist in the
interface.

This way, the shear stresses (in both the bonded and debonded interfaces)
between the broken tow and its neighbors are expressed in piecewise forms
based on the status of the interfaces:

$\uptau^{(k,0)}_{\Leftrightarrow}(x)=\left\\{\begin{array}[]{lr}\delta
u^{(k,0)}_{11}(x)\cdot
K_{s,\Leftrightarrow}\>,&d_{\Leftrightarrow}=0\vphantom{\bigg{|}}\\\
\mu_{\Leftrightarrow}\cdot\langle-\sigma^{(k)}_{22,\infty}\rangle\>,&d_{\Leftrightarrow}=1\vphantom{\big{|}}\end{array}\right.\>,\qquad\uptau^{(k,0)}_{\Updownarrow}(x)=\left\\{\begin{array}[]{lr}\delta
u^{(k,0)}_{11}(x)\cdot
K_{s,\Updownarrow}\>,&d_{\Leftrightarrow}=0\vphantom{\bigg{|}}\\\
\mu_{\Updownarrow}\cdot\langle-\bar{\sigma}^{(k)}_{33}\rangle\>,&d_{\Leftrightarrow}=1\vphantom{\bigg{|}}\end{array}\right.$
(12)

where $\langle\,\rangle$ denotes the Macaulay bracket:

$\langle x\rangle=\max(0,x)=\begin{cases}x\;,&x\geq 0\\\ 0\;,&x<0\end{cases}$
(13)

while $\mu_{\Leftrightarrow}$ and $\mu_{\Updownarrow}$ are the corresponding
friction coefficients.

###  2.5 Governing equations and boundary conditions

Since interfacial shear tractions depends on the interface status (Eq.12),
specific expressions of the force equilibrium (Eq.7) also have to be written
in a piecewise form. Thus it is necessary to determine the inequality/equality
relation between the debonding lengths, $L_{\Leftrightarrow}$ and
$L_{\Updownarrow}$, to obtain the domain and expression of differential
equations, note that the debonding lengths are also unknown values to be
solved.

The force equilibrium of the broken tow (Eq.7) and the physical meaning of
perturbation displacement (Eq.6) indicates that $\delta u^{(k,0)}_{11}$ is a
non-negative, non-increasing function whose magnitude vanishes as
$x\to+\infty$. If the shear tractions are small, the stress recover length
would be long, and $\delta u^{(k,0)}_{11}$ decrease gradually; if the shear
tractions are large, stresses caused by the broken tow would recover to far-
field values in a shorter distance.

As shown in Fig.6, since debonding will only occur upon reaching the interface
strength, and all interfaces in the analytical analysis are assumed in pure
shear mode, the resultant boundary condition requires that the elastic shear
stresses reach the corresponding strengths at the ’debonding tips’, the
expressions are:

$\displaystyle\uptau^{(k,0)}_{\Leftrightarrow}(x)\bigg{|}_{x=L^{+}_{\Leftrightarrow}}$
$\displaystyle=\delta u^{(k,0)}_{11}(L^{+}_{\Leftrightarrow})\cdot
K_{s,\Leftrightarrow}=\tau^{f}_{s,\Leftrightarrow}$ (14)
$\displaystyle\uptau^{(k,0)}_{\Updownarrow}(x)\bigg{|}_{x=L^{+}_{\Updownarrow}}$
$\displaystyle=\delta u^{(k,0)}_{11}(L^{+}_{\Updownarrow})\cdot
K_{s,\Updownarrow}=\tau^{f}_{s,\Updownarrow}$ (15)

where $\tau^{f}_{s,\Leftrightarrow}$ and $\tau^{f}_{s,\Updownarrow}$ are the
elastic strengths of intra-ply and inter-ply interfaces in pure shear mode,
respectively. The superscript ${}^{\prime}+^{\prime}$ on $L_{\Leftrightarrow}$
and $L_{\Updownarrow}$ is used to show that function/value need to be
evaluated at the bonded domain where damage variables are zero.

Therefore, based on the non-decreasing feature of $\delta u^{(k,0)}_{11}$ and
the conditions given in Eq.14, it is observed that the relation between
$L_{\Leftrightarrow}$ and $L_{\Updownarrow}$ depends on the interfacial pure-
mode shear strengths and corresponding shear stiffness, which can be
categorized into three cases:

$\begin{cases}\text{Case 1:}\quad
L_{\Leftrightarrow}>L_{\Updownarrow}\>,&\text{if}\quad\tau^{f}_{s,\Leftrightarrow}/K_{s,\Leftrightarrow}<\tau^{f}_{s,\Updownarrow}/K_{s,\Updownarrow}\vphantom{\bigg{|}}\\\
\text{Case 2:}\quad
L_{\Leftrightarrow}<L_{\Updownarrow}\>,&\text{if}\quad\tau^{f}_{s,\Leftrightarrow}/K_{s,\Leftrightarrow}>\tau^{f}_{s,\Updownarrow}/K_{s,\Updownarrow}\vphantom{\bigg{|}}\\\
\text{Case 3:}\quad
L_{\Leftrightarrow}=L_{\Updownarrow}\>,&\text{if}\quad\tau^{f}_{s,\Leftrightarrow}/K_{s,\Leftrightarrow}=\tau^{f}_{s,\Updownarrow}/K_{s,\Updownarrow}\vphantom{\bigg{|}}\end{cases}$
(16)

Now the piecewise governing equations can be written for all cases. Take Case
1 as example, the intervals are $\rinterval{0}{L_{\Updownarrow}}$,
$\rinterval{L_{\Updownarrow}}{L_{\Leftrightarrow}}$,
$\rinterval{L_{\Leftrightarrow}}{\infty}$, and the corresponding governing
equations are:

$\displaystyle E_{l}\cdot wh\cdot\dfrac{d^{2}}{dx^{2}}\delta
u^{(k,0)}_{11}(x)=$ $\displaystyle
2h\cdot\mu_{\Leftrightarrow}\cdot\langle{-}\sigma^{(k)}_{22,\infty}\rangle+2w\cdot\mu_{\Updownarrow}\cdot\langle{-}\bar{\sigma}^{(k)}_{33}\rangle\>,$
$x\in\rinterval{0}{L_{\Updownarrow}}\vphantom{\bigg{|}}$ (17a) $\displaystyle
E_{l}\cdot wh\cdot\dfrac{d^{2}}{dx^{2}}\delta u^{(k,0)}_{11}(x)=$
$\displaystyle
2h\cdot\mu_{\Leftrightarrow}\cdot\langle{-}\sigma^{(k)}_{22,\infty}\rangle+2w\cdot
K_{s,\Updownarrow}\cdot\delta u^{(k,0)}_{11}(x)\>,$
$x\in\rinterval{L_{\Updownarrow}}{L_{\Leftrightarrow}}\vphantom{\bigg{|}}$
(17b) $\displaystyle E_{l}\cdot wh\cdot\dfrac{d^{2}}{dx^{2}}\delta
u^{(k,0)}_{11}(x)=$ $\displaystyle 2h\cdot K_{s,\Leftrightarrow}\cdot\delta
u^{(k,0)}_{11}(x)+2w\cdot K_{s,\Updownarrow}\cdot\delta u^{(k,0)}_{11}(x)\>,$
$x\in\rinterval{L_{\Leftrightarrow}}{+\infty}\vphantom{\bigg{|}}$ (17c)

which contain three second-order differential equations. Now define the
following constants:

$\displaystyle A_{1}=E_{l}(w\cdot h),\quad
A_{2}=2h\cdot\mu_{\Leftrightarrow}\cdot\langle-\sigma^{(k)}_{22,\infty}\rangle,\quad
A_{3}=2w\cdot\mu_{\Updownarrow}\cdot\langle-\bar{\sigma}^{(k)}_{33}\rangle,\quad$
(18) $\displaystyle A_{4}=2h\cdot K_{s,\Leftrightarrow},\quad A_{5}=2w\cdot
K_{s,\Updownarrow}$

which are correlated with material properties and in-situ stresses, and are
all non-negative constants given the corresponding physical meanings. Then by
setting $\delta u\equiv\delta u^{(k,0)}_{11}(x)$ for simplicity, the governing
equations for Case 1 can be expressed as:

$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{1}{A_{1}}(A_{2}+A_{3}),$
$x\in\rinterval{0}{L_{\Updownarrow}}\vphantom{\bigg{|}}$ (19a)
$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{1}{A_{1}}(A2+A5\cdot\delta u)\>,$
$x\in\rinterval{L_{\Updownarrow}}{L_{\Leftrightarrow}}\vphantom{\bigg{|}}$
(19b) $\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{A4+A5}{A_{1}}\cdot\delta u,$
$x\in\rinterval{L_{\Leftrightarrow}}{+\infty}\vphantom{\bigg{|}}$ (19c)

There are eight undetermined values in the governing equations, two of which
are the unknown debonding lengths that will also determine the intervals of
piecewise function. There are exactly eight constraints to solve for these
unknown values:

  1. (a)

Axial stress (in fiber orientation) of the broken tow is zero at plane of tow
break ($x=0$).

  2. (b)

Axial stress recovers to the far-field level after a sufficiently large
distance ($x\to+\infty$).

  3. (c)

Axial stress is continuous at domain-dividing point $x=L_{\Leftrightarrow}$ of
the piecewise function.

  4. (d)

Axial stress is continuous at domain-dividing point $x=L_{\Updownarrow}$ of
the piecewise function.

  5. (e)

Perturbation displacement is continuous at $x=L_{\Leftrightarrow}$.

  6. (f)

Perturbation displacement is continuous at $x=L_{\Updownarrow}$.

  7. (g)

Intra-ply interfacial shear stress reaches debonding strength
$\tau^{f}_{s,\Updownarrow}$ at ’debonding tip’ $L^{+}_{\Updownarrow}$.

  8. (h)

Inter-ply interfacial shear stress reaches debonding strength
$\tau^{f}_{s,\Leftrightarrow}$ at ’debonding tip’ $L^{+}_{\Leftrightarrow}$.

The expressions of these eight constraints are:

where the superscript $+$ or $-$ on domain-dividing points
($L_{\Leftrightarrow}$ and $L_{\Updownarrow}$) indicate the positive or
negative side of the interval to evaluate the piecewise function. For example,
$x=L^{+}_{\Leftrightarrow}$ means to evaluate Eq.19c at point
$L_{\Leftrightarrow}$, whereas $x=L^{-}_{\Leftrightarrow}$ means to evaluate
Eq.19b at point $L_{\Leftrightarrow}$.

###  2.6 Solutions of the shear-lag models

To solve for all eight unknowns, it is found that a polynomial equation of the
debonding lengths have to be obtained first, then the remaining six unknowns
can be expressed by either one of the debonding lengths. For Case 3, the total
number of unknowns reduces to six but the solution procedure is similar.

Stay with Case 1 as an example, the general solution takes form:

$\delta u^{(k,0)}_{11}(x)=\begin{cases}\dfrac{A_{2}+A_{3}}{2A_{1}}\cdot
x^{2}-\dfrac{\sigma^{(k)}_{11,\infty}}{E_{l}}\cdot
x+C_{2}\>,&x\in\rinterval{0}{L_{\Updownarrow}}\vphantom{\bigg{|}}\\\
-A_{2}/A_{5}+C_{3}\cdot\exp\left(\sqrt{A_{5}/A_{1}}\cdot
x\right)+C_{4}\cdot\exp\left(-\sqrt{A_{5}/A_{1}}\cdot
x\right)\;,&x\in\rinterval{L_{\Updownarrow}}{L_{\Leftrightarrow}}\vphantom{\bigg{|}}\\\
C_{6}\cdot\exp\left(-\sqrt{(A_{4}+A_{5})/A_{1}}\cdot
x\right)\;,&x\in\rinterval{L_{\Leftrightarrow}}{+\infty}\vphantom{\bigg{|}}\end{cases}$
(21)

The solutions for debonding lengths are:

$\displaystyle L_{\Updownarrow}$
$\displaystyle=-\dfrac{1}{2A_{8}}\left(A_{9}+\sqrt{A_{9}^{2}-4A_{8}A_{10}}\right)>0\,,\vphantom{\bigg{|}}$
(22) $\displaystyle L_{\Delta}$ $\displaystyle\equiv
L_{\Updownarrow}-L_{\Leftrightarrow}=\frac{1}{2}\sqrt{A_{1}/A_{5}}\cdot\ln\left(\dfrac{A_{6}}{A_{7}}\cdot\dfrac{A_{1}B_{3}+B_{1}L_{\Updownarrow}}{A_{1}B_{2}-B_{1}L_{\Updownarrow}}\right)<0\,,\vphantom{\bigg{|}}$
$\displaystyle L_{\Leftrightarrow}$
$\displaystyle=L_{\Updownarrow}-L_{\Delta}>L_{\Updownarrow}\vphantom{\bigg{|}}$

where $L_{\Delta}$ is a term showing the difference between two debonding
lengths (negative in Case 1), and the remaining four undetermined constants
can be expressed by the two debonding lengths:

$\displaystyle
C_{2}=-\dfrac{A_{2}}{A_{5}}+\dfrac{1}{2A_{5}K_{s,\Leftrightarrow}}\left[A_{6}\exp\left(-\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)+A_{7}\exp\left(\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)\right]-\dfrac{A_{2}+A_{3}}{2A_{1}}\cdot
L_{\Updownarrow}^{2}+\dfrac{\sigma^{(k)}_{11,\infty}}{E_{l}}\cdot
L_{\Updownarrow}\,,\vphantom{\bigg{|}}$ (23) $\displaystyle
C_{3}=\dfrac{A_{7}}{2A_{5}K_{s,\Leftrightarrow}}\exp\left(-\sqrt{A_{5}/A_{1}}\cdot
L_{\Leftrightarrow}\right)\,,\quad
C_{4}=\dfrac{A_{6}}{2A_{5}K_{s,\Leftrightarrow}}\exp\left(\sqrt{A_{5}/A_{1}}\cdot
L_{\Leftrightarrow}\right)\,,\vphantom{\bigg{|}}$ $\displaystyle
C_{6}=\dfrac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}\exp\left(\sqrt{(A_{4}+A_{5})/A_{1}}\cdot
L_{\Leftrightarrow}\right)$

where both Eq.22 and Eq.23 contain eight extra constants ($A_{6}$ to $A_{10}$,
$B_{1}$ to $B_{3}$), which are dependent values of $A_{1}$ to $A_{5}$ and the
far-field stress $\sigma^{(k)}_{11,\infty}$. Detailed expressions for these
constants and the essential steps to obtain solutions are discussed in
appendices.

To strictly ensure uniqueness and physical validity of the shear-lag analysis,
there are two required conditions on stress states. The first involves the
far-field tensile stress level at which the tow breaks:

$\left(\sigma^{(k)}_{11,\infty}\right)^{2}>\left(\sigma^{*}_{11,\infty}\right)^{2}\equiv\frac{A_{2}E_{l}^{2}}{A_{1}}\left(\frac{A_{2}}{A_{5}}+2\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)-\frac{E_{l}^{2}}{A_{1}}\left[\frac{A_{6}A_{7}}{A_{5}K_{s,\Leftrightarrow}^{2}}-A_{5}\left(\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)^{2}\right]$
(24a) where the defined term, $\sigma^{*}_{11,\infty}$, is the minimum
threshold value of far-field stress. The second condition involves the
frictional sliding forces in the debonded/damaged interfaces:
$2h\cdot\mu_{\Leftrightarrow}\cdot\langle-\sigma^{(k)}_{22,\infty}\rangle+2w\cdot\mu_{\Updownarrow}\cdot\langle-\bar{\sigma}^{(k)}_{33}\rangle=A_{2}+A_{3}>0$
(24b)

The physical meaning of the first condition (Eq.24a) is: if the far-field
stress level at which the tow breaks is relatively low, then only the
interfacial elastic shear stresses are adequate to counterbalance and recover
the axial tensile stress lost by the broken tow, such that the ’debond-and-
slip’ interfacial failure mode would not take place. The physical meaning of
the second condition is: there must be at least some sort of sliding
(frictional) forces in the debonded interfaces, so that the broken tow will
reach equilibrium after a sufficiently large distance. Otherwise, assume that
the far-field stress is large enough to cause interface damage/debonding
(Eq.24a is satisfied) but no frictional forces are accumulated in the debonded
interfaces ($A_{2}+A_{3}=0$), then the interfacial debonding will propagate
continuously such that the damage evolution would become unstable, and no
static equilibrium state would be reached. Mathematically, these two
conditions will ensure that the solutions of debonding lengths are real and
positive values. Derivation of Eq.24 is shown in A. It is also verified that
the solution is uniquely determined when these two conditions are met, the
proof is presented in B.

###  2.7 A discussion on the interface debonding phenomena

In the above section, we explain the physical and mathematical meanings of the
two required conditions (Eq.24) to ensure the validity of shear-lag analysis
and the uniqueness of solutions. It is then readily seen that, given the
typical dimension and modulus of prepreg tows, the required minimum threshold
value of far-field stress ($\sigma^{*}_{11,\infty}$ in Eq.24a) is
approximately two orders of magnitude lower than the typical ultimate tensile
stress of tows. This implies that interface debonding is essentially a certain
damage mode when a tensile tow break occurs. Namely, the released energy and
perturbation forces caused by the broken tow is sufficiently large to trigger
debonding along the resin-rich interfaces around the broken tow.

This phenomena can also be explained by fracture-mechanics based studies
involving a crack at (or approaching) the interface between two continuum
domains [36, 37, 38, 39, 40, 41]. In such scenarios, the crack may penetrate
into the neighboring domain (across interface) or deflect and turn into the
bi-material interface. And there is a ’competition’ between crack deflection
and penetration when determining the path and mode of crack propagation.

A conclusion particularly useful for this study is that when the elastic
properties of two continuum domains are identical, crack deflection is almost
surely to occur (instead of penetration) if the strength or fracture toughness
of the continuum material is more than around four times the mixed-mode
toughness of the interface between two domains [36, 39]. This circumstance
applies here because the longitudinal stiffness and fracture toughness of
composite tows/strands are typically at least ten times higher than the
properties of epoxy resins. Nevertheless, this does not mean that the adjacent
tows will not break: tensile failure in neighboring tows (crack penetration)
can definitely occur at certain locations if the overall stress level is high;
also, interface cracks can possibly turn and penetrate into tows at some ’weak
spots’ caused by imperfections or flaws along the tow edges. As mentioned,
studies involving multiple tow breaks are not in the focus of this paper,
instead, it should be performed once a study of a single tow break is
established.

###  2.8 Summary

In this section, a shear-lag based analysis of an in-situ tow break in a
representative uni-directional filament-wound composite is presented. By
considering different inter- and intra-ply debonding lengths, the governing
equations are categorized accordingly into different domains, whereby the
solutions of axial stress and perturbation displacement of the broken tow are
expressed in a piecewise form. Besides, essential variables such as far-field
stress states, interfacial shear strengths and frictional coefficients, are
included in the analysis to justify the consideration of in-situ effects. In
addition, to ensure that the shear-lag analysis is fundamentally meaningful,
two required conditions on the in-situ stress states are also determined. The
physical meanings of these two conditions are explained, and it is proved that
the solutions are uniquely determined when they are satisfied. The general
solution for the other two cases are given in C and D.

##  3 Stress concentration profile around an in-situ tow break

Following the shear-lag analysis, it remains an essential task to evaluate the
overload profile and stress concentration factor (SCF) around locations of
interest. Due to the complexity of such a problem, closed form analytical
solutions of the stress fields are arguably impractical, and yet it is still a
challenging task to develop reasonable approximations.

In this study, we apply a novel approach that combines classic solid mechanics
solutions with the shear-lag analysis results from the previous section.
Specifically, the idea is to: 1) firstly treat the interfacial shear stresses
as ’loads’ that generate/exert overload onto the neighboring domains, by
considering the fact that interfacial shear stresses are the main ’load
transferor’ that redistribute the overload to its surrounding intact
neighbors; and then 2) predict the stress concentration field by accounting
for the interfacial shear tractions based on their magnitude and locations.
Detailed steps and derivations are demonstrated in this section.

###  3.1 Stress field due to a tangential point load on the surface of a
transversely-isotropic half space

The resultant stress field by tangential point load applied on the surface of
a transversely-isotropic infinite half space is firstly investigated. In
particular, the expressions adopted here is based on the derivations by Lin et
al. [42]. An earlier work by Pan and Chou [43] using Green’s function
solutions was also investigated, but it seems that several constants
(’$k_{i}$’ in Eq.30 in their study) were either undefined or misprinted (also
implicitly mentioned by Kachanov et al. [44]).

Figure 7: Schematic plot of tangential load applied on the surface of an
transversely isotropic elastic infinite half space. (A) Tangential point load
applied at $(\xi_{1},\xi_{2},0)$, with resultant stress component
$\sigma_{11}$ at a point of interest $(\chi_{1},\chi_{2},\chi_{3})$ given by
Eq.25. (B) Distributed tangential point load applied over a rectangular
region. Eq.30 gives the resultant stress component $\sigma_{11}$ at a point
$(\chi_{1},\chi_{2},\chi_{3})$, where $(\eta_{1},\eta_{2},0)$ is the centroid
of the rectangular region, and $l_{1},l_{2}$ are its lengths in
$\vec{e_{1}},\vec{e_{2}}$ directions, respectively.

As illustrated in Fig.7A, consider a unit tangential point load
$\vec{F}=1\vec{e_{1}}$ applied at a point $(\xi_{1},\xi_{2},0)$ on the surface
of a transversely-isotropic, infinite 3D half space ($z>0$). The resultant
axial stress along $\vec{e_{1}}$ at a point of interest
$(\chi_{1},\chi_{2},\chi_{3})$ is:

$\Omega_{11}(\chi_{1},\chi_{2},\chi_{3}\,;\,\xi_{1},\xi_{2})=2c_{66}a_{13}\left(\dfrac{R_{2}}{D_{3}^{*}}\right)_{,12}+\sum_{\alpha=1}^{2}a_{1\alpha}\theta_{\alpha}\left[2c_{66}\left(\dfrac{R_{2}}{D_{\alpha}^{*}}\right)_{,12}+\dfrac{g_{\alpha}R_{1}}{D_{\alpha}^{3}}\right]$
(25)

in which the derivative term has form:

$\left(\dfrac{R_{\alpha}}{D_{j}^{*}}\right)_{,12}=R_{\alpha}R_{1}R_{2}\left(\dfrac{1}{D_{j}^{3}(D_{j}^{*})^{2}}+\dfrac{2}{D_{j}^{2}(D_{j}^{*})^{3}}\right)-\dfrac{\delta_{\alpha
2}R_{1}+\delta_{\alpha 1}R_{2}}{D_{j}(D_{j}^{*})^{2}}\>,\qquad\alpha=1,2,\quad
j=1,2,3$ (26)

where $\delta_{\alpha\beta}$ is the Kronecker delta operator. The expressions
for all stress components can be found in Refs.[42, 44] (notations may vary).
Here only $\sigma_{11}$ is of primary interest. The variables in $\Omega_{11}$
that involve coordinates are:

$\displaystyle R_{\alpha}=\chi_{\alpha}-\xi_{\alpha}\>,\quad
D_{j}=\sqrt{R_{1}^{2}+R_{2}^{2}+(\theta_{j}\chi_{j})^{2}}\>,\quad
D_{j}^{*}=D_{j}+\theta_{j}\chi_{j}\>,\qquad\alpha=1,2,\quad j=1,2,3$ (27)

whereas the properties only involve stiffness constants are:

$\displaystyle\theta_{1}=\left[\dfrac{\sqrt{c_{11}c_{33}}+c_{13}+2c_{44}}{4c_{33}c_{44}(\sqrt{c_{11}c_{33}}-c_{13})^{-1}}\right]^{\frac{1}{2}}+\left[\dfrac{\sqrt{c_{11}c_{33}}-c_{13}-2c_{44}}{4c_{33}c_{44}(\sqrt{c_{11}c_{33}}+c_{13})^{-1}}\right]^{\frac{1}{2}}\>,$
(28)
$\displaystyle\theta_{2}=\left[\dfrac{\sqrt{c_{11}c_{33}}+c_{13}+2c_{44}}{4c_{33}c_{44}(\sqrt{c_{11}c_{33}}-c_{13})^{-1}}\right]^{\frac{1}{2}}-\left[\dfrac{\sqrt{c_{11}c_{33}}-c_{13}-2c_{44}}{4c_{33}c_{44}(\sqrt{c_{11}c_{33}}+c_{13})^{-1}}\right]^{\frac{1}{2}}\>,$
$\displaystyle\theta_{3}=\left(\dfrac{c_{66}}{c_{44}}\right)^{\frac{1}{2}}\>,\quad
h_{1}=\dfrac{c_{11}-c_{44}\theta_{1}^{2}}{c_{13}+c_{44}}\>,\quad
h_{2}=\dfrac{c_{11}-c_{44}\theta_{2}^{2}}{c_{13}+c_{44}}\>,$ $\displaystyle
b_{1}=\theta_{1}(c_{13}-c_{33}h_{1})\>,\quad
b_{2}=\theta_{2}(c_{13}-c_{33}h_{2})\>,$ $\displaystyle
g_{1}=c_{11}-c_{13}h_{1}\>,\quad g_{2}=c_{11}-c_{13}h_{2}\>,$ $\displaystyle
a_{11}=\dfrac{b_{2}}{2\pi(b_{1}g_{2}-b_{2}g_{1})}\>,\quad
a_{12}=\dfrac{-b_{1}}{2\pi(b_{1}g_{2}-b_{2}g_{1})}\>,\quad
a_{13}=\dfrac{1}{2\pi c_{44}\theta_{3}}$

in which $c_{ij}$ are components of the compliance matrix of an transversely
isotropic material.

The point load solution adopted here requires that $\vec{e_{3}}$ is normal to
the plane of isotropy, i.e., elastic properties in $\vec{e_{1}}$ and
$\vec{e_{2}}$ directions are identical. Whereas for the purpose of this work,
$\vec{e_{1}}$ is expected to be the longitudinal fiber/tow orientation along
which the modulus is much higher than the two transverse directions. However,
a particularly important concern in this study is that the tow break is ’in-
situ’, wherein the following aspects need to be properly treated: 1) Surfaces
of the broken tow are either bonded or remain in contact with neighboring
domains, and thus, planes where interfacial shear tractions are applied may
not be considered as completely free half surfaces. 2) Sizes of tow break and
interface cracks are relatively small compared to typical dimension of the
component-scale filament structure, therefore a tow break are not expected to
cause large deformations in the surround composite along the transverse
directions. 3) Winding patterns may lead to complicated local stacking layups,
so that the fiber/tow orientations of adjacent layers may align closely along
the transverse direction of the center layer containing the tow break (for
instance, a $90^{\circ}$ winding layer above the center $0^{\circ}$ layer).
Therefore, to address the above issue and ensure compatibility with Eq.25, the
elastic properties in $\vec{e_{1}}$ and $\vec{e_{2}}$ directions are assumed
to be identical and comparable to the longitudinal modulus of the tows
$E_{l}$, while the properties in $\vec{e_{3}}$ direction are equivalent to the
transverse modulus of tows $E_{t}$.

This way, the compliance matrix takes the form:

$\begin{bmatrix}\varepsilon_{11}\\\ \varepsilon_{22}\\\ \varepsilon_{33}\\\
\gamma_{23}\\\ \gamma_{31}\\\
\gamma_{12}\end{bmatrix}=\begin{bmatrix}1/E_{l}&-\nu_{ll}/E_{l}&-\nu_{tl}/E_{t}&0&0&0\\\
-\nu_{ll}/E_{l}&1/E_{l}&\nu_{tl}/E_{t}&0&0&0\\\
-\nu_{lt}/E_{l}&-\nu_{lt}/E_{l}&1/E_{t}&0&0&0\\\ 0&0&0&1/G_{tl}&0&0\\\
0&0&0&0&1/G_{tl}&0\\\
0&0&0&0&0&2(1+\nu_{ll})/E_{l}\end{bmatrix}\begin{bmatrix}\sigma_{11}\\\
\sigma_{22}\\\ \sigma_{33}\\\ \sigma_{23}\\\ \sigma_{31}\\\
\sigma_{12}\end{bmatrix}$ (29)

where the properties are as defined in Eq.3. Note that the newly introduced
Poisson’s ratio, $\nu_{ll}$, has to be determined empirically. Briefly, the
concept of $\nu_{ll}$ is the in-plane Poisson’s ratio of a thin plate stacked
with sufficiently large amount of tows in randomly distributed orientations,
therefore a reasonable approximation is to take the typical value of balanced
woven fabric composite made of materials with comparable modulus. For
instance, an approximated value for balanced Carbon/Epoxy woven composites is
$\nu_{ll}=0.045$, and $0.13$ for E-Glass [45] composites, respectively.
Nevertheless, the effect of picked values for Poisson’s ratios is found to be
negligible comparing to other properties. Then upon definition of the
compliance matrix, the expressions for the required stiffness constants in
Eq.28 can be readily obtained by taking inverse of the compliance matrix in
Eq.29.

###  3.2 Stress field due to uniformly distributed tangential load over a
rectangular region on the surface of a transversely-isotropic half space

For the purpose of this study, a specific interest is the axial stress field
(in $\vec{e_{1}}$ direction) caused by uniformly distributed tangential load
over a rectangular region, because this is the ’shape’ of debonded domain
among which interfacial shear tractions occur. As schematically shown in
Fig.7B, by integrating the point load solution (Eq.25) over a rectangular
region [42], and then by applying necessary coordinate transformation, the
resultant axial stress caused by tangential load over a rectangular region on
the surface of half space can be expresses as:

$\Lambda_{11}(\chi_{1},\chi_{2},\chi_{3}\,;\,\eta_{1},\eta_{2}\,;\,l_{1},l_{2})=\quad\Omega_{11}(\chi_{1}{+}\eta_{1},\chi_{2}{+}\eta_{2},\chi_{3}\,;\,{-}\frac{l_{1}}{2},{-}\frac{l_{2}}{2})+\Omega_{11}(\chi_{1}{+}\eta_{1},\chi_{2}{+}\eta_{2},\chi_{3}\,;\,\frac{l_{1}}{2},\frac{l_{2}}{2})\\\
-\Omega_{11}(\chi_{1}{+}\eta_{1},\chi_{2}{+}\eta_{2},\chi_{3}\,;\,{-}\frac{l_{1}}{2},{+}\frac{l_{2}}{2})-\Omega_{11}(\chi_{1}{+}\eta_{1},\chi_{2}{+}\eta_{2},\chi_{3}\,;\,\frac{l_{1}}{2},{-}\frac{l_{2}}{2})$
(30)

where $(\xi_{1},\xi_{2},0)$ is the coordinate of centroid of the rectangular
region, and $l_{1},l_{2}$ are its lengths in the $\vec{e_{1}}$ and
$\vec{e_{2}}$ directions, respectively.

###  3.3 Critical stress concentration profile along intra-ply neighboring
tows

By utilizing the above formulas, overload caused by the broken tow can be
estimated by accounting for the contributions from interfacial shear
tractions. In particular, the following procedures are used to compute
overload profile at a point of interest:

  * •

Take each of the eight debonded interfaces around the broken tow,
respectively, as the ’surface of transversely isotropic elastic half space’.

  * •

Based on the normal directions of all debonded interfaces and the position of
the point of interest, apply coordinate transformation, and then compute the
resultant overload stress caused by interfacial frictional shear forces in
each debonded interface, respectively.

  * •

By combining together the contribution from each debonded interface, the
overload stress at any point can be obtained.

Note that the solution of stress field ($\Omega_{11}$ in Eq.25) have
singularities at the location of point load. After applying integration over a
rectangular region ($\Lambda_{11}$ in Eq.30), logarithmic singularity still
exists at the boundaries (corner and edges) of the rectangular region,
therefore, direct evaluation based on numerical integration is not valid at
those isolated locations.

Nevertheless, the main interest in this study is the ’critical’ stress
concentration profile in the plane of tow break ($x=0$), or more specifically,
those along the y- and z-axis, because these are the ’plane’ and ’paths’ where
SCF are most prominent. Since the crack opening displacement is positive
($\delta_{cod}>0$), overload profile in the entire tow break plan (including
the y- and z-axis) does not have singularities. Besides, since the shear
stress in the deformed-but-bonded interfaces are away from tow breakage and
also decay rapidly, the contribution from elastic shear stresses in the bonded
interfaces are neglected from the quantification of overload profile if the
point of interest is in the plane of tow break.

Following the above discussion, the ’critical’ stress magnitude in the intra-
ply neighboring tows along the y-axis can be evaluated by:

$\sigma^{(k)}_{11,\text{yaxis}}|_{y>w/2}=\sigma^{(k)}_{11,\infty}+\mu_{\Leftrightarrow}\cdot\langle{-}\sigma^{(k)}_{22,\infty}\rangle\cdot\Lambda^{\Sigma,\text{intra}}_{11,\text{pert},y}+\mu_{\Updownarrow}\cdot\langle{-}\bar{\sigma}^{(k)}_{33}\rangle\cdot\Lambda^{\Sigma,\text{inter}}_{11,\text{pert},y}$
(31)

where $\sigma^{(k)}_{11,\infty}$ and $\sigma^{(k)}_{22,\infty}$ are the far-
field longitudinal and transverse stress magnitudes in the local coordinate of
tows, $\bar{\sigma}^{(k)}_{33}$ is the average through-thickness stress across
the $k$th ply (as defined in Eq.1).

In Eq.31, $\Lambda^{\Sigma,\text{intra}}_{11,\text{pert}}$ and
$\Lambda^{\Sigma,\text{inter}}_{11,\text{pert}}$ are the
${}^{\prime}11^{\prime}$ component of overload (perturbation) stresses along
$y$ axis caused by unit shear tractions in the intra-ply and inter-ply
interfaces, respectively, which take forms:

$\displaystyle\begin{split}\Lambda^{\Sigma,\text{intra}}_{11,\text{pert},y}={}&\Lambda_{11}(0,0,y{-}\frac{w}{2}\,;\,\delta_{cod}{+}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},h)-\Lambda_{11}(0,0,y{-}\frac{w}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},h)\\\
&+\Lambda_{11}(0,0,y{+}\frac{w}{2}\,;\,\delta_{cod}{+}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},h)-\Lambda_{11}(0,0,y{+}\frac{w}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},h)\end{split}$
(32) $\displaystyle\Lambda^{\Sigma,\text{inter}}_{11,\text{pert},y}={}$
$\displaystyle
2\left[\Lambda_{11}(0,y,\frac{h}{2}\,;\,\delta_{cod}{+}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)-\Lambda_{11}(0,y,\frac{h}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)\right]$
(33)

where the input variables in stress functions ($\Lambda_{11}$ in Eq.30) are
replaced with the material and geometrical parameters correlated to the
modulus of tows and solved deformations from the shear-lag algorithms. The
four components in $\Lambda^{\Sigma,\text{intra}}_{11,\text{pert}}$ refer to
the four debonded interface regions in the same ply ($y=\pm w/2$ planes),
where coordinate involving $y$ are placed as the third variable due to
definition of $\Lambda_{11}$ and coordinate transformation. The two terms in
$\Lambda^{\Sigma,\text{inter}}_{11,\text{pert}}$ refer to the two debonded
interface regions between the broken tow and its neighboring ply on the top
($z=h/2$) planes, and since the debonded domains in the above and below
surfaces of the broken tow (in $z=\pm h/2$ planes) are symmetric about the
$z=0$ plane, a coefficient of 2 is applied. Difference in signs of
$\Lambda_{11}$ are due to the opposite directions of the shear tractions
relative to the positive direction of y-axis.

###  3.4 Critical stress concentration profile along inter-ply neighboring tow
regions

For the stress redistribution between intra-ply tows as shown above, it is
relatively straightforward since tows in the same ply are parallel to each
other. In contrast, tow orientations (winding angles) can vary across
different plies, and hence quantification of overload profile requires extra
consideration. This is an important aspect for typical FWCs because winding
patterns not only determine the effective properties of composites, but also
affect the crack propagation and damage mechanisms.

As shown schematically in Fig.3C, tows in other plies are regarded as ’tow
regions’ if winding angles are different. Despite that, critical stress
concentration is still expected in the plane of broken tow along z-axis. So
the idea here is to: 1) compute the ’reference’ overload profile for the case
of uni-directional tows where all plies are assumed to have the same winding
angle as the $k\text{th}$ ply, and then 2) to predict the ’real’ overload
profile based on basis transformation corresponding to a certain winding
angle.

Therefore, by firstly assuming that all winding layers are locally uni-
directional, the stress concentration field (denoted by superscript ’scf,uni’)
along the z-axis, can be expressed as:

$\sigma^{\text{scf,uni}}_{11,\text{zaxis}}|_{z>h/2}=\mu_{\Leftrightarrow}\cdot\langle{-}\sigma^{(k)}_{22,\infty}\rangle\cdot\Lambda^{\Sigma,\text{intra}}_{11,\text{pert,z}}+\mu_{\Updownarrow}\cdot\langle{-}\bar{\sigma}^{(k)}_{33}\rangle\cdot\Lambda^{\Sigma,\text{inter}}_{11,\text{pert,z}}$
(34)

where $\Lambda^{\Sigma,\text{intra}}_{11,\text{pert,z}}$ and
$\Lambda^{\Sigma,\text{inter}}_{11,\text{pert,z}}$ are the
${}^{\prime}11^{\prime}$ component of overload (perturbation) stresses along
z-axis caused by unit shear tractions in the intra-ply and inter-ply
interfaces around the broken tow, respectively, which take forms:

$\displaystyle\begin{split}\Lambda^{\Sigma,\text{intra}}_{11,\text{pert,z}}={}&\Lambda_{11}(0,0,z{-}\frac{h}{2}\,;\,+\delta_{cod}{+}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)-\Lambda_{11}(0,0,z{-}\frac{h}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)\\\
&+\Lambda_{11}(0,0,z{+}\frac{h}{2}\,;\,{+}\delta_{cod}{+}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)-\Lambda_{11}(0,0,z{+}\frac{h}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Updownarrow}}{2},0\,;\,L_{\Updownarrow},w)\end{split}$
(35) $\displaystyle\Lambda^{\Sigma,\text{inter}}_{11,\text{pert},z}={}$
$\displaystyle 2\left[\
\;\Lambda_{11}(0,z,\frac{w}{2}\,;\,\delta_{cod}{+}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},h)-\Lambda_{11}(0,z,\frac{w}{2}\,;\,{-}\delta_{cod}{-}\frac{L_{\Leftrightarrow}}{2},0\,;\,L_{\Leftrightarrow},w)\right]$
(36)

Then for the $i$th ply ($i\neq k$), the stress concentration in the local
fiber orientation, denoted by superscript ’scf,local’, is approximated by
simple change of basis:

$\sigma^{\text{scf,local}}_{11,\text{zaxis}}|_{z\in[z_{i},z_{i+1}]}=\cos\theta^{i}\cdot\sigma^{\text{scf,uni}}_{11,\text{zaxis}}|_{z\in[z_{i},z_{i+1}]}$
(37)

where $\theta^{i}\in[-90^{\circ},90^{\circ}]$ is the valid interval for ply
winding angles. It is straight forward to see that for $\pm 90^{\circ}$ plies,
the stress concentration in local fiber orientation becomes zero, and reaches
maximum when the winding angle is the same as broken tow $0^{\circ}$.

###  3.5 Summary

In this section, by combining: 1) the solution of interface debonding lengths
and tow break opening displacement from shear-lag analysis, 2) solutions of
resultant stresses caused by tangential loads applied on the surface of a
transversely-isotropic linear elastic half space, and 3) superposition concept
and coordinate transformations, we derived formulas to quantify critical
stress concentration profiles in the plane of tow break. Such an approach
requires more calculations than simple empirical formulas, but the discrete
cracks and damage patterns are treated in a more ’high-fidelity’ way. In the
following section, the proposed approach will be compared with a series of
finite element studies.

##  4 Finite element validation using AGDM

Numerical simulations to predict tow break opening, debonding lengths, and
stress concentration profiles are carried out by implementing an auto-
generated geometry-based discrete finite element model (AGDM) recently
proposed by Liu and Phoenix [46], which has been developed in terms of Abaqus
processing commands, and hence is compatible with Abaqus/CAE environment. AGDM
realizes accurate coupling between all three types of typical damage modes in
laminated composite structures (fiber/yarn break, matrix crack, and
delamination) having arbitrary stacking layup angles, and hence it is adopted
here to effectively model tow/yarn break, debonding between yarns (matrix
crack), and inter-ply debonding (delamination).

Briefly speaking, in order to represent a crack in an explicit way and ensure
accurate displacement field across the crack tip, two types of approaches may
be implemented: 1) a solid element must be partitioned into sub-domains, or 2)
an embedded cohesive element (between solid elements) must deform and fail, to
ensure real discontinuity and fracture surfaces in the numerical space. For
interaction and coupling of multiple cracks, these discontinuity (fracture
surfaces) must have matched mesh at the potential intersecting/deflecting
locations [46]. Unlike typical approaches that rely on bonded mesh to realize
matching of mesh, AGDM is created as an assembly of discrete continuum domains
and small-thickness interfaces with matching of nodes at potential crack
coupling locations, and thereby ensures accurate displacement jump across
crack tips and coupling between cracks without requiring priory-generated
(bonded) solid elements. This advantage firstly helps to create numerical
model for laminates with arbitrary stacking layup angles, and secondly
provides the flexibility to use extra fine meshes (for cohesive and/or
continuum elements) to capture high-fidelity stress gradients. Details on the
algorithms and generation procedures of AGDM can be found in Refs.[46, 47].

An earlier numerical study attempting to quantify SCF around a broken tow was
also reported by Liu and Phoenix [13], but only limited success was obtained
in single-layer models (without inter-ply debonding) due to non-matched mesh
between inserted cohesive elements at potential crack deflection locations.
Consequently, the multi-ply scenarios were not included due to converging
issues and artificial stress profiles mostly caused by ambiguous interpolation
rules between continuum and cohesive domains. Later on, AGDM was developed
after observing the limitations in existing methods, and thereafter is
implemented in this study to verify the newly derived analytical approach in a
high-fidelity way.

###  4.1 Material properties and numerical inputs

The properties of prepreg tows used for this study are based on data sheet of
a commercially available, high strength Toray® T1100g 12k tows. The
mechanical, dimensional, and interfacial properties of T1110 layers are listed
in Table 1, where the same nomenclature as in previous sections is used. The
cross-sectional area is approximated based on the nominal void volume fraction
and number of filaments, width and height of tow cross section is then
obtained by assuming a width-to-height ratio of $3:1$. It should be noted that
several parameters (such as the shear modulus of tows $G_{lt}$) are estimated
based on reported value of comparable materials. And the main purpose of FEA
here is to compare with the prediction from the derived algorithms in Section
3 using the same input values.

The interfacial behaviors (i.e., traction-separation rules in bonded
interfaces, and pressure-dependent sliding in debonded interfaces,
respectively) as explained in Section 2.4 for intra-ply and inter-ply cohesive
elements are realized by user material subroutines (VUMAT).

Table 1: Mechanical, dimensional, and interface properties of T1100g 12k resin-impregnated tows |  Definition | Variable | Value
---|---|---|---
Mechanical constants | Tensile modulus of tow in local longitudinal orientation | $E_{l}$ | 185.0 GPa
Tensile modulus of tow in transverse orientations | $E_{t}$ | 10.0 GPa
Shear modulus of tow in longitudinal-to-transverse orientation | $G_{lt}$ | 6.5 GPa
Poisson’s ratio of tow in longitudinal-to-transverse orientation | $\nu_{lt}$ | 0.25
Poisson’s ratio of tow in transverse-to-transverse orientation | $\nu_{tt}$ | 0.40
Approximated planar Poisson’s ratio of randomly stacked plies | $\nu_{ll}$ | 0.045
Cross-section dimension | Representative cross-sectional area | $A_{tow}$ |  $0.39$ $\text{mm}^{2}$
Representative cross-sectional width | $w$ |  $1.09$ mm
Representative cross-sectional height (along stacking direction) | $h$ |  $0.36$ mm
Intra-ply properties | Interface penalty stiffness | $K_{s,\Leftrightarrow}$ |  $10^{6}$ N/$\text{mm}^{3}$
Interface strength in pure normal mode (Mode-I) | $\tau^{f}_{n,\Leftrightarrow}$ | 80.0 MPa
Interface strength in pure shear modes (Mode-II or -III) | $\tau^{f}_{s,\Leftrightarrow}$ | 100.0 MPa
Interface frictional coefficient (after debonding) | $\mu_{\Leftrightarrow}$ | 0.30
Inter-ply properties | Interface penalty stiffness | $K_{s,\Updownarrow}$ |  $10^{6}$ N/$\text{mm}^{3}$
Interface strength in pure normal mode (Mode-I) | $\tau^{f}_{n,\Updownarrow}$ | 80.0 MPa
Interface strength in pure shear modes (Mode-II or -III) | $\tau^{f}_{s,\Updownarrow}$ | 120.0 MPa
Interface frictional coefficient (after debonding) | $\mu_{\Updownarrow}$ | 0.45

###  4.2 Model setup and boundary conditions

The representative FWC laminate is 100mm long along the x-axis (also the fiber
orientation), a length found to be sufficiently large such that the
longitudinal stress in broken tow recovers to far-field level before reaching
boundaries of the model. There are 4 neighboring tows on each side of the
broken tow in the center, so the width of the laminate (dimension along the
y-axis) is $w\times 9=9.81mm$. There are in total 13 plies, 6 above and 6
below the $k$th ply, respectively, so the total thickness of laminate is
$h\times 13=4.68mm$. The center (origin) of the tow break plane is placed to
overlap with the origin of global Cartesian coordinate system (Fig.4 and
Fig.5).

Fig.8A shows the continuum tows in the FE model. Fig.8B shows the embedded
intra- and inter-ply small-thickness interfaces (to be meshed to cohesive
elements) among tows to capture potential interface debonding. Fig.8C contains
two figures showing the mesh of continuum and cohesive domains, respectively.
Since the tow break occurs at the center of the 7th ply, tows in the 5th to
9th plies are all assigned with a smaller element size to ensure that there
are at least eight fully-integrated solid elements across a tow’s width, and
at least two layers of elements across a tow’s height. This helps to capture
high-fidelity stress gradients around the broken tow, as will be shown in the
next subsection.

Since the present study is focused on the equilibrium state whereby any
dynamic effect is neglected, the broken tow is represented by an initial slit
(discontinuity) with extra small gap ($10^{-6}$mm, see Fig.8D). This way,
after the model reaches assumed stress states, the resultant debonding
lengths, break opening displacement, and surrounding stress fields are
equivalent to the static equilibrium state where a tow break would result at
the same stress states. More specifically, under the framework of AGDM, one
single tow breakage interface domain at the center of the 6th ply is embedded
into the center tow. This interface is then removed from the numerical model
so that the model initially contains a ’closed’ tow break. Then during
numerical analysis, as the far-field stresses increases, the tow break
surfaces will open up and lead to debonding due to excessive interfacial shear
deformations.

Figure 8: Finite element model and mesh. (A) Continuum tows (no mesh) with
partitions along the $x=0$ and $z=0$ planes to ensure existing of element
nodes in these planes (after mesh) for evaluation of critical overload
profile. (B) Embedded small-thickness interfaces (no mesh) between the
continuum tows for potential interfacial debonding and frictions. (Ci) Mesh of
continuum tows, where the tows in the plies close to the center ply are
assigned with smaller element sizes to realize fine stress gradient. (Cii)
Mesh of embedded interface domains (into cohesive elements) between the
continuum tows. (D) Mesh around the broken tow. In this study, breakage
(fracture) is accounted for by a discontinuity (removed cohesive element) in
the center of the broken tow. The initial gap between the fracture surfaces is
$10^{-6}$mm. Matching of mesh is ensured to realize accurate crack deflection
between the tow break and both intra-ply and inter-ply interfaces.

The boundary conditions are determined based on the targeted far-field stress
magnitudes (Eq.1). Once these targeted stresses are determined, the
corresponding strains are calculated based on the effective properties of the
representative laminate (Eq.29). Then the targeted strains are applied in
terms of velocity load on all the element nodes on the outer surfaces of the
laminate model. This method of applying boundary conditions is specifically
adopted since it is found to alleviate kinetic energy as compared to directly
applying stresses. Kinetic energy is monitored throughout the numerical
analysis, and ensured to be always less than $1\%$ of the strain energy of the
model, to justify that the system can be considered as quasi-static.

Assume that when the tow breaks, the tensile stress along the fiber
orientation is $\sigma^{(k)}_{11,\infty}=1000\text{MPa}$. And the in-situ
transverse stresses ($\sigma^{(k)}_{22,\infty}$, $\sigma^{(k)}_{33,\infty}$)
are picked from the possible combinations of three values: -50, -100, and
-150MPa, resulting in 9 cases in total (Table.2). Note that
$\bar{\sigma}^{(k)}_{33}$ is in general negative for FWC structures (e.g.,
COPVs) due to applied internal pressure load. For uni-directional tows in this
study, the transverse stress $\sigma^{(k)}_{22,\infty}$ is also considered to
be negative here due to Poisson’s effect. In reality, a piece of a
representative overwrap model may be under bi-axial tension.

###  4.3 Examples of FEA results

Here we briefly show several results on deformation, stress redistribution and
concentration, and also interface debonding from FEA. For example, Fig.9A
shows slice views at the mid-plane ($z=0$) of the model, it can be observed
that the stress lost by the broken tow is redistributed to the neighboring
units, and the tensile stress gradually recovers to the far-field level.
Fig.9B shows the predicted interface debonding in terms of failed and damaged
cohesive elements, where an element shown in red indicates complete failure
(its damage variable reaches one). The debonding lengths can thus be obtained
by measuring the length of the failed elements from the plane of tow break to
the furthest x-coordinate with failed elements. Fig.9C contains two examples
showing the slice views at plane of two break, it can be seen that the in-situ
stress states affect the stress redistribution mechanisms and subsequently the
overload profile and SCFs.

Figure 9: Finite element results showing deformation, interface debonding, and
stress profile around the broken tow. (A) Slice view of the deformation around
the tow break at the $z=0$ plane, showing break opening and axial stress
profile ($\sigma_{11}$). (B) Predicted interface debonding in the form of
failed cohesive elements, wherein the failed elements (with damage variable
equals one) are shown in red, and the intact elements are shown in blue. (C)
Slice views at plane of tow break ($x=0$ plane) with two different in-situ
stress states, showing effect of the stress states on the stress concentration
profile and factors.

###  4.4 Comparison between numerical and analytical predictions

To verify the accuracy and robustness of the derived algorithms, we compare
the analytical predictions with the numerical results. A specific interest is
the overload profile and SCF in the plane of tow break, because SCFs are
expected to be most prominent in this plane, and henceforth the maximum errors
between FEA and theory would also occur in this plane since stresses will
decay from this plane. More specifically, the y-axis and z-axis are expected
to be the ’path’ along which the intra- and inter-ply stress concentrations
are most significant, respectively.

Fig.10 shows a comparison between the predicted tensile stress along fiber
orientation ($\sigma_{11}$) in the tow break plane for Case A.1
($\sigma^{(k)}_{22,\infty}=-50$, $\sigma^{(k)}_{33,\infty}=-50$). The stress
profile predicted by theory (Eq.30) agrees reasonable well with the prediction
by FEA. Moreover, the maximum stresses (and hence the maximum SCF) also
matches closely with each other.

Figure 10: Comparison between FEA and theory on the predicted tensile stress
in tow break plane, with in-situ stresses $\sigma^{(k)}_{11,\infty}=1000$,
$\sigma^{(k)}_{22,\infty}=-50$, and $\bar{\sigma}^{(k)}_{33}=-50$ (Case A.1 in
Table.2). (A) Stress profile by FEA. (B) Predicted stress profile by
analytical approach (Eq.30).

Fig.11 shows comparison between SCFs (predicted by FEA and theory) along
y-axis and x-axis for all cases. As expected, maximum stresses occur in
regions closest to the broken tow, and decrease when moving away from the tow
breakage, and then eventually drop back to far-field levels. It is also
observed that the assigned boundary conditions result in exactly the targeted
far-field stress levels (e.g., $\sigma^{(k)}_{11,\infty}=1000\text{MPa}$ in
fiber orientation), which verifies the accuracy of using velocity-based
boundary conditions in FEA to obtain the targeted stress levels for each case.

Figure 11: Comparison between FEA and theory on the predicted stress
concentration magnitude along the positive y-axis and z-axis. Markers are from
the element nodal values in FEA, continuous lines are predictions by theory.
(Ai) SCF along y-axis covering two adjacent tow widths. (Aii) Magnified view
of comparison in (Ai) showing SCF close to the broken tow (about a quarter of
tow width). (Bi) SCF along the z-axis (stacking direction). (Bii) Magnified
view of data in (Bi) showing SCF close to the broken tow.

In addition to the SCFs, Table.2 shows comparison between the predictions of
several important values, max SCF, tow break opening displacement, debonding
length, for all cases. In which, the max SCF values (and corresponding errors)
are evaluated at closest coordinates to the broken tow ($y=w/2$, $z=h/2$,
leftmost data in Fig.11), and henceforth the error computed at these two
points are also where the maximum error occurs. From Fig.11 and Table.2, it is
found that excellent agreement is realized in all cases. The errors tend to
increase moderately if one or both of the transverse in-situ stresses
($\sigma^{(k)}_{22,\infty}$, $\sigma^{(k)}_{33,\infty}$) are higher. Note that
the large transverse in-situ stress values (100MPa, 150MPa for
$\sigma^{(k)}_{33,\infty}$) are picked to validate the robustness of theory
and FEA. Whereas in practice, typical magnitudes of inner pressures for COPVs
are less than 10000 Psi ($\approx$ 72MPa).

Table 2: Numerical case studies on uni-directional plies

Case | In-situ stress states (MPa) |  | Max SCF |  Break opening |  Ratio between debonding lengths
---|---|---|---|---|---
| $\sigma^{(k)}_{11,\infty}$ | $\sigma^{(k)}_{22,\infty}$ | $\bar{\sigma}^{(k)}_{33}$ |  | Intra-ply | Inter-ply | 2$\delta_{cod}$ | $L_{\Leftrightarrow}/L_{\Updownarrow}$
A.1 | 1000 | -50 | -50 | Theory | 1.160 | 1.251 | 0.0348 | 0.952
FEA | 1.175 | 1.277 | 0.0352 | 0.949
Difference | 1.30$\%$ | 2.09$\%$ | 1.05$\%$ | -0.28$\%$
A.2 | 1000 | -50 | -100 | Theory | 1.203 | 1.548 | 0.0191 | 0.913
FEA | 1.194 | 1.472 | 0.0201 | 1.00
Difference | -0.77$\%$ | -4.92$\%$ | 5.00$\%$ | 9.56$\%$
A.3 | 1000 | -50 | -150 | Theory | 1.225 | 1.861 | 0.0132 | 0.873
FEA | 1.203 | 1.624 | 0.0148 | 1.030
Difference | -1.81$\%$ | -12.70$\%$ | 12.16$\%$ | 17.92$\%$
A.4 | 1000 | -100 | -50 | Theory | 1.289 | 1.260 | 0.029 | 0.932
FEA | 1.264 | 1.278 | 0.0283 | 0.882
Difference | -1.93$\%$ | 1.40$\%$ | -4.05$\%$ | -5.33$\%$
A.5 | 1000 | -100 | -100 | Theory | 1.346 | 1.556 | 0.0174 | 0.885
FEA | 1.283 | 1.470 | 0.0176 | 0.857
Difference | -4.66$\%$ | -5.53$\%$ | 1.13$\%$ | -3.19$\%$
A.6 | 1000 | -100 | -150 | Theory | 1.377 | 1.868 | 0.0123 | 0.838
FEA | 1.291 | 1.622 | 0.0133 | 0.896
Difference | -6.26$\%$ | -13.16$\%$ | 7.78$\%$ | 6.98$\%$
A.7 | 1000 | -150 | -50 | Theory | 1.377 | 1.868 | 0.0123 | 0.838
FEA | 1.291 | 1.622 | 0.0133 | 0.896
Difference | -4.67$\%$ | 1.09$\%$ | -5.73$\%$ | -9.81$\%$
A.8 | 1000 | -150 | -100 | Theory | 1.494 | 1.563 | 0.0160 | 0.791
FEA | 1.369 | 1.466 | 0.0158 | 0.796
Difference | -8.33$\%$ | -6.17$\%$ | -0.91$\%$ | 0.66$\%$
A.9 | 1000 | -150 | -150 | Theory | 1.531 | 1.874 | 0.0116 | 0.712
FEA | 1.370 | 1.617 | 0.0119 | 0.810
Difference | -10.56$\%$ | -13.72$\%$ | 2.75$\%$ | 13.86$\%$

  * *

Total number of plies is 13, tow breakage occurs in the center tow in the 7th
ply.

  * *

Errors of max SCF are evaluated at $(x=0,y=w/2,z=0)$ and $(x=0,y=0,z=h/2)$ for
intra- and inter-ply domains, respectively, because these points are where
maximum errors occur.

###  4.5 Summary

In this section, we performed high-fidelity finite element simulations to
evaluate the stress concentration profile around a tow break under different
stress states, and compared the numerical results with the predictions by the
derived theories. Overall, excellent agreement in essential metrics are
observed (Table 2). It should be emphasized that neither FEA results nor
analytical predictions should be considered as the ’correct reference values’
because closed-form solutions are arguably not attainable, and henceforth
under this circumstance, quantitatively close agreement between these two
approaches validate the effectiveness of the derived algorithms.

Despite the transverse in-situ stresses are negative in all cases in this
study, there are no specific restrictions on the sign and magnitude of them in
the analytical algorithms. For example, if there is zero intra-ply frictions
after debonding, then by theory, virtually all of the overload will be
distributed into the adjacent plies through inter-ply frictional forces. Such
a ’zero slipping force’ scenario may not be practical since the filament
overwrap should remain as an integrated part until destructive failure. Also,
after debonding, the frictional forces may not necessarily be linearly
proportional to the contact pressure. It may also depend on other factors such
as surface roughness, stacking layup angles, etc. For example, it is observed
in an experimental study (also Toray T1100) that the inter-tow angles directly
affect the frictional behaviors [48].

In specific, the primary reason for using negative transverse stresses in FEA
is to ensure numerical stability. Since the interfaces are represented by the
small-thickness cohesive elements embedded between continuum tows, zero or
positive in-situ transverse stresses (on the laminate boundaries) will lead to
splitting between tows, which would further result in local instabilities.
This is because the deformation mechanism of the cohesive elements is based on
the relative displacements between pairs of nodes on its top and bottom
surfaces, and subsequently, a cohesive element having both positive normal
separation and shear deformation would tend to rotate and distort, which leads
to dynamic effects and stress oscillations across the debonded interfaces.
Non-negative interfacial friction in this study also requires that the top and
bottom surfaces of a cohesive element is in contact with negative normal
forces, and henceforth loading cases with zero or positive
$\sigma^{(k)}_{22,\infty},\sigma^{(k)}_{33,\infty}$ are not investigated in
the finite element simulations. To alleviate the dynamic effect caused by
distortion of cohesive elements, specially modified cohesive zone modelling
techniques may be required, which is subjective to future studies.

##  5 Conclusion

In this study, we derived an empirical analytical approach to predict the
overload profile around a broken tow in FWC structures. This kind of problem
is challenging due to complicated geometries, interfacial debonding and
sliding between domains, and also the load-redistribution mechanisms under
different stress states. To resolve this, we firstly applied a shear-lag
analysis to solve for tow break opening, perturbation axial displacement of
the broken tow, intra- and inter-ply debonding lengths (Section 2).
Subsequently, based on the first principals of the load-sharing mechanisms
between different domains, we derived formulas that approximate the resultant
overload stresses around the broken tow under assumed in-situ stress levels
(Section 3). Thereafter, we compared the analytical predictions of critical
SCFs under different stress states to the values predicted by high-fidelity
finite element analysis, and observed excellent agreement between them
(Section 4).

The derived algorithms conveniently predict the overload stress magnitude at a
location of interest in FWCs, and henceforth, could be used as a complementary
evaluation at meso scale to improve theoretical statistical strengths
prediction of such composites. On the other side, the algorithms could be
implemented in development of failure criteria for continuum elements in
component-scale finite element models. Also, this framework may be adapted to
structures with similar formations, such as tape-based composites, layered
structures printed by direct ink writing, etc.

There are also limitations in this study. Firstly, the representative cross-
section of tows is assumed to be rectangular, whereas in practice it can be
somewhere between oval to rectangular. Secondly, the characterization of
resin-rich areas between tows, despite being straightforwardly observable from
failure pattern of FWC structures, is somehow empirical, and may need to be
further verified by experimental and tomographic studies. Thirdly, due to the
potential numerical instabilities in debonded interfaces (cohesive elements)
when the normal strain is zero or positive, the in-situ transverse stresses in
the case studies are set to negative values (Table.2). That said, the
effectiveness of the derived algorithms have been verified by varying of the
stresses states, and thus the algorithms can be used for general stress states
despite lack of corresponding numerical validation at the moment. Moreover,
the interfacial response is characterized as linear traction-separation before
debonding, followed by pressure-dependent sliding after debonding (Fig.6,
while in practice epoxy resin interfaces can exhibit certain extent of non-
linearity (before debonding). Lastly, inter-ply overload magnitudes for FWCs
with general stacking layup angles are predicted by applying multiplication of
trigonometric terms to the results on uni-directional tows (Eq.37, which may
need to be verified by additional numerical studies in the future.

Potential subsequent work can be categorized into two directions. The first is
specific applications of the proposed analytical algorithms, including the
aforementioned development of a theoretical framework for statistical strength
prediction of filament-wound composites. The second direction involves
refinement and additional validation of the present study, such as stress
concentration profiles between plies with different winding angles.

## Acknowledgments

A new empirical analytical approach is developed for predicting the stress
concentration profile around an in-situ tow break in filament-wound
composites. A shear-lag analysis is firstly performed to solve for the
perturbational axial displacement of the broken tow and resultant debonding
lengths. Solution of stress field caused by tangential load on the surface of
a elastic half space is then utilized in combination with superposition
concepts to obtain the overload magnitudes in the neighboring tows.
Subsequently, high-fidelity finite element analysis on a representative uni-
directional laminate model under different stress states is performed, and
excellent overall agreement is observed between analytical and numerical
predictions. The proposed method takes into account essential aspects such as
transversely isotropic material properties, in-situ stress states and their
effect on the interfacial frictional forces in the debonded interfaces, and
thus provides a convenient way to evaluate the stress concentration factors in
damaged filament-wound composites. In addition, this approach can be applied
to yield auxiliary failure evaluation criteria for statistical strength
prediction or finite element modeling of filament-wound composites or similar
structures. The authors acknowledge funding support from research grant
provided by National Institute of Standards and Technology (NIST) under
agreement ID: 70NANB14H323. The first author also had received support from
Cornell Graduate School.

## Data availability

Code for analytical predictions, input files (.inp) for FEA (Abaqus/Explicit),
and corresponding user subroutines (VUMAT) used for this study are available
from the authors upon request.

## References

  * [1] P. Liu, J. Zheng, Recent developments on damage modeling and finite element analysis for composite laminates: A review, Materials & Design 31 (8) (2010) 3825–3834.  doi:https://doi.org/10.1016/j.matdes.2010.03.031.
  * [2] M. Zhang, H. Lv, H. Kang, W. Zhou, C. Zhang, A literature review of failure prediction and analysis methods for composite high-pressure hydrogen storage tanks, International Journal of Hydrogen Energy 44 (47) (2019) 25777–25799.  doi:https://doi.org/10.1016/j.ijhydene.2019.08.001.
  * [3] W. Li, H. Lv, L. Zhang, P. He, C. Zhang, Experiment, simulation, optimization design, and damage detection of composite shell of hydrogen storage vessel-a review, Journal of Reinforced Plastics and Composites 42 (11-12) (2023) 507–536.  doi:https://doi.org/10.1177/07316844221132744.
  * [4] L. K. Grunenfelder, T. Centea, P. Hubert, S. R. Nutt, Effect of room-temperature out-time on tow impregnation in an out-of-autoclave prepreg, Composites Part A: Applied Science and Manufacturing 45 (2013) 119–126.  doi:10.1016/j.compositesa.2012.10.001.
  * [5] A. Cherniaev, I. Telichev, Experimental and numerical study of hypervelocity impact damage in composite materials fabricated by filament winding, International Journal of Impact Engineering 98 (2016) 19–33.  doi:10.1016/j.ijimpeng.2016.07.001.
  * [6] S. G. Schechter, L. K. Grunenfelder, S. R. Nutt, Design and application of discontinuous resin distribution patterns for semi-pregs, Advanced Manufacturing: Polymer and Composites Science 6 (2) (2020) 72–85.  doi:10.1080/20550340.2020.1736864.
  * [7] J. Zheng, P. Liu, Elasto-plastic stress analysis and burst strength evaluation of al-carbon fiber/epoxy composite cylindrical laminates, Computational Materials Science 42 (3) (2008) 453–461.  doi:https://doi.org/10.1016/j.commatsci.2007.09.011.
  * [8] L. Wang, C. Zheng, S. Wei, Z. Wei, Micromechanics-based progressive failure analysis of carbon fiber/epoxy composite vessel under combined internal pressure and thermomechanical loading, Composites Part B: Engineering 89 (2016) 77–84.  doi:10.1016/j.compositesb.2015.11.018.
  * [9] S. Alam, G. Yandek, R. C. Lee, J. Mabry, A study of residual burst strength of composite over wrapped pressure vessel due to low velocity impact, International Journal of Pressure Vessels and Piping 194 (2021) 104511.  doi:https://doi.org/10.1016/j.ijpvp.2021.104511.
  * [10] McLaughlan P. B., Grimes-Ledesma L. R., Composite Overwrapped Pressure Vessels, A Primer, Tech. Rep. March (2011).
  * [11] J. H. S. Almeida, M. L. Ribeiro, V. Tita, S. C. Amico, Damage and failure in carbon/epoxy filament wound composite tubes under external pressure: Experimental and numerical approaches, Materials and Design 96 (2016) 431–438.  doi:10.1016/j.matdes.2016.02.054.
  * [12] L. Zu, H. Xu, H. Wang, B. Zhang, B. Zi, Design and analysis of filament-wound composite pressure vessels based on non-geodesic winding, Composite Structures 207 (2019) 41–52.  doi:10.1016/j.compstruct.2018.09.007.
  * [13] J. Liu, S. L. Phoenix, Meso-scopic model for progressive failure of resin impregnated yarns in composite overwrapped pressure vessels (COPVs), Journal of Space Safety Engineering 6 (4) (2019) 225–234.  doi:10.1016/j.jsse.2019.10.007.
  * [14] J. M. Hedgepeth, Stress concentrations in filamentary structures, National Aeronautics and Space Administration, 1961.
  * [15] J. M. Hedgepeth, P. Van Dyke, Local stress concentrations in imperfect filamentary composite materials, Journal of composite materials 1 (3) (1967) 294–309.  doi:https://doi.org/10.1177/002199836700100305.
  * [16] I. J. Beyerlein, S. L. Phoenix, Stress concentrations around multiple fiber breaks in an elastic matrix with local yielding or debonding using quadratic influence superposition, Journal of the Mechanics and Physics of Solids 44 (12) (1996) 1997–2039.  doi:10.1016/S0022-5096(96)00068-3.
  * [17] I. J. Beyerlein, S. L. Phoenix, A. M. Sastry, Comparison of shear-lag theory and continuum fracture mechanics for modeling fiber and matrix stresses in an elastic cracked composite lamina, International Journal of Solids and Structures 33 (18) (1996) 2543–2574.  doi:10.1016/0020-7683(95)00172-7.
  * [18] M. Ibnabdeljalil, W. A. Curtin, Strength and reliability of fiber-reinforced composites: Localized load-sharing and associated size effects, International Journal of Solids and Structures 34 (21) (1997) 2649–2668.  doi:10.1016/S0020-7683(96)00179-5.
  * [19] I. Babuška, B. Andersson, P. J. Smith, K. Levin, Damage analysis of fiber composites Part I: Statistical analysis on fiber scale, Computer Methods in Applied Mechanics and Engineering 172 (1) (1999) 27–77.  doi:10.1016/S0045-7825(98)00225-4.
  * [20] C. M. Landis, R. M. McMeeking, A shear-lag model for a broken fiber embedded in a composite with a ductile matrix, Composites Science and Technology 59 (3) (1999) 447–457.  doi:https://doi.org/10.1016/S0266-3538(98)00091-8.
  * [21] C. M. Landis, R. M. McMeeking, Stress concentrations in composites with interface sliding, matrix stiffness and uneven fiber spacing using shear lag theory, International Journal of Solids and Structures 36 (28) (1999) 4333–4361.  doi:10.1016/S0020-7683(98)00193-0.
  * [22] S. Mahesh, S. L. Phoenix, I. J. Beyerlein, Strength distributions and size effects for 2D and 3D composites with Weibull fibers in an elastic matrix, International Journal of Fracture 115 (1) (2002) 41–85.  doi:10.1023/A:1015729607223.
  * [23] S. Mahesh, J. C. Hanan, E. Üstündag, I. J. Beyerlein, Shear-lag model for a single fiber metal matrix composite with an elasto-plastic matrix and a slipping interface, International Journal of Solids and Structures 41 (15) (2004) 4197–4218.  doi:10.1016/j.ijsolstr.2004.02.050.
  * [24] Y. Swolfs, H. Morton, A. E. Scott, L. Gorbatikh, P. A. Reed, I. Sinclair, S. M. Spearing, I. Verpoest, Synchrotron radiation computed tomography for experimental validation of a tensile strength model for unidirectional fibre-reinforced composites, Composites Part A: Applied Science and Manufacturing 77 (2015) 106–113.  doi:10.1016/j.compositesa.2015.06.018.
  * [25] A. Engelbrecht-Wiggans, S. L. Phoenix, A stochastic model based on fiber breakage and matrix creep for the stress-rupture failure of unidirectional continuous fiber composites, International Journal of Fracture 217 (1-2) (2019) 1–34.  doi:10.1007/s10704-019-00359-9.
  * [26] S. L. Phoenix, I. J. Beyerlein, 1.19 - Statistical Strength Theory for Fibrous Composite Materials, in: A. Kelly, C. Zweben (Eds.), Comprehensive Composite Materials, Pergamon, Oxford, 2000, pp. 559–639.  doi:10.1016/B0-08-042993-9/00056-5.
  * [27] J. C. Thesken, P. L. Murthy, S. Phoenix, N. Greene, J. L. Palko, J. Eldridge, J. Sutter, R. Saulsberry, H. Beeson, A theoretical investigation of composite overwrapped pressure vessel (copv) mechanics applied to nasa full scale tests, Tech. rep. (2009).
  * [28] S. Pradhan, A. Hansen, B. K. Chakrabarti, Failure processes in elastic fiber bundles, Reviews of Modern Physics 82 (1) (2010) 499–555.  arXiv:0808.1375, doi:10.1103/RevModPhys.82.499.
  * [29] P. Xu, J. Zheng, P. Liu, Finite element analysis of burst pressure of composite hydrogen storage vessels, Materials & Design 30 (7) (2009) 2295–2301.  doi:https://doi.org/10.1016/j.matdes.2009.03.006.
  * [30] T. Q. Hua, H. S. Roh, R. K. Ahluwalia, Performance assessment of 700-bar compressed hydrogen storage for light duty fuel cell vehicles, International Journal of Hydrogen Energy 42 (40) (2017) 25121–25129.  doi:10.1016/j.ijhydene.2017.08.123.
  * [31] J. P. Canal, A. Micuzzi, H. Logarzo, A. Terlisky, R. Toscano, E. Dvorkin, On the finite element modeling of COPVs, Computers and Structures 220 (2019) 1–13.  doi:10.1016/j.compstruc.2019.04.007.
  * [32] R. Pramod, C. K. Krishnadasan, N. Siva Shanmugam, Design and finite element analysis of metal-elastomer lined composite over wrapped spherical pressure vessel, Composite Structures 224 (2019) 111028.  doi:10.1016/j.compstruct.2019.111028.
  * [33] Q. Zhang, H. Xu, X. Jia, L. Zu, S. Cheng, H. Wang, Design of a 70 mpa type iv hydrogen storage vessel using accurate modeling techniques for dome thickness prediction, Composite Structures 236 (2020) 111915.  doi:https://doi.org/10.1016/j.compstruct.2020.111915.
  * [34] R. Rafiee, F. Abbasi, Numerical and experimental analyses of the hoop tensile strength of filament-wound composite tubes, Mechanics of Composite Materials 56 (2020) 423–436.  doi:https://doi.org/10.1007/s11029-020-09894-2.
  * [35] Z. Xia, W. Curtin, P. Peters, Multiscale modeling of failure in metal matrix composites, Acta materialia 49 (2) (2001) 273–287.  doi:https://doi.org/10.1016/S1359-6454(00)00317-7.
  * [36] Mi.-Y. He, J. W. Hutchinson, Crack deflection at an interface between dissimilar elastic materials, International Journal of Solids and Structures 25 (9) (1989) 1053–1067.  doi:10.1016/0020-7683(89)90021-8.
  * [37] V. Gupta, A. S. Argon, Z. Suo, Crack deflection at an interface between two orthotropic media, Journal of Applied Mechanics 59 (1992) 79–87.  doi:10.1115/1.2899511.
  * [38] D. L. Tullock, I. E. Reimanis, A. L. Graham, J. J. Petrovic, Deflection and penetration of cracks at an interface between two dissimilar materials, Acta Metallurgica Et Materialia 42 (9) (1994) 3245–3252.  doi:10.1016/0956-7151(94)90422-7.
  * [39] M. Y. He, A. G. Evans, J. W. Hutchinson, Crack deflection at an interface between dissimilar elastic materials: Role of residual stresses, International Journal of Solids and Structures 31 (24) (1994) 3443–3455.  doi:10.1016/0020-7683(94)90025-6.
  * [40] J. P. Parmigiani, M. D. Thouless, The roles of toughness and cohesive strength on crack deflection at interfaces, Journal of the Mechanics and Physics of Solids 54 (2) (2006) 266–287.  doi:10.1016/j.jmps.2005.09.002.
  * [41] E. Martin, B. Poitou, D. Leguillon, J. M. Gatt, Competition between deflection and penetration at an interface in the vicinity of a main crack, International Journal of Fracture 151 (2) (2008) 247–268.  doi:10.1007/s10704-008-9228-0.
  * [42] W. Lin, C. H. Kuo, L. M. Keer, Analysis of a transversely isotropic half space under normal and tangential loadings, Journal of Tribology 113 (2) (1991) 335–338.  doi:10.1115/1.2920625.
  * [43] Y.-C. Pan, T.-W. Chou, Point Force Solution for an Infinite Transversely Isotropic Solid, Journal of Applied Mechanics 43 (4) (1976) 608–612.  doi:10.1115/1.3423941.
  * [44] M. Kachanov, B. Shafiro, I. Tsukrov, Handbook of Elasticity Solutions, Springer Science & Business Media, 2003.  doi:10.1007/978-94-017-0169-3.
  * [45] D. Gay, Composite Materials: Design and Applications, Third Edition, 3rd Edition, CRC Press, Boca Raton, 2014.  doi:10.1201/b17106.
  * [46] J. Liu, S. L. Phoenix, An auto-generated geometry-based discrete finite element model for damage evolution in composite laminates with arbitrary stacking sequence, Composites Part C: Open Access 3 (2020) 100066.  doi:https://doi.org/10.1016/j.jcomc.2020.100066.
  * [47] J. Liu, A discrete modeling approach for progressive failure analysis of composite laminates and filament-wound pressure vessels, Ph.D. thesis, Cornell University (2020).
  * [48] M. Tourlonias, M. A. Bueno, G. Fassi, I. Aktas, Y. Wielhorski, Influence of friction angle between carbon single fibres and tows: Experimental analysis and analytical model, Composites Part A: Applied Science and Manufacturing 124 (Sep. 2019).  doi:10.1016/j.compositesa.2019.105478.

###### List of Figures

  1. 1 Basic procedures for resin-impregnation of carbon fiber tows and filament winding process (this plot shows wet-winding).
  2. 2 Schematic plot of material formation in a resin-impregnated tow and dimensions of its representative cross section.
  3. 3 Coordinate system and notations for plies and tows. (A) Intra-ply slice view ($z=0$ plane) around the tow break, showing notations of the neighboring tows in the $k$th ply. (B) Through-thickness slice view in the plane of tow break ($x=0$ plane), showing notations of the neighboring tow regions in different plies. (C) If the adjacent plies have different winding angles, the effective tow region may cover a portion of multiple tows. Here as an example, assume the first ply above the broken tow has a winding angle of $-30^{\circ}$, then superscript ($k+1,j$) corresponds to an effective region covering multiple tows.
  4. 4 Intra-ply cut-view of an in-situ tow break occurred at center of the $k$th ply with winding angle $\theta^{k}=0^{\circ}$.
  5. 5 Inter-ply cross-sectional cut-view of an in-situ tow break and corresponding interface debonding.
  6. 6 Schematic plot of interfacial behavior in pure shear modes (Mode-II, Mode-III), where $d$ is the effective damage variable of the interface, $K$ is the penalty stiffness, $\delta^{f},\tau^{f}$ are the failure (debonding) separation and traction, respectively. (A) Typical traction-separation response, where the interface debonding occurs upon reaching the failure criterion. (B) If a debonded interfacial region is still in contact and under negative normal strain, pressure-dependent friction may exist in the interface.
  7. 7 Schematic plot of tangential load applied on the surface of an transversely isotropic elastic infinite half space. (A) Tangential point load applied at $(\xi_{1},\xi_{2},0)$, with resultant stress component $\sigma_{11}$ at a point of interest $(\chi_{1},\chi_{2},\chi_{3})$ given by Eq.25. (B) Distributed tangential point load applied over a rectangular region. Eq.30 gives the resultant stress component $\sigma_{11}$ at a point $(\chi_{1},\chi_{2},\chi_{3})$, where $(\eta_{1},\eta_{2},0)$ is the centroid of the rectangular region, and $l_{1},l_{2}$ are its lengths in $\vec{e_{1}},\vec{e_{2}}$ directions, respectively.
  8. 8 Finite element model and mesh. (A) Continuum tows (no mesh) with partitions along the $x=0$ and $z=0$ planes to ensure existing of element nodes in these planes (after mesh) for evaluation of critical overload profile. (B) Embedded small-thickness interfaces (no mesh) between the continuum tows for potential interfacial debonding and frictions. (Ci) Mesh of continuum tows, where the tows in the plies close to the center ply are assigned with smaller element sizes to realize fine stress gradient. (Cii) Mesh of embedded interface domains (into cohesive elements) between the continuum tows. (D) Mesh around the broken tow. In this study, breakage (fracture) is accounted for by a discontinuity (removed cohesive element) in the center of the broken tow. The initial gap between the fracture surfaces is $10^{-6}$mm. Matching of mesh is ensured to realize accurate crack deflection between the tow break and both intra-ply and inter-ply interfaces.
  9. 9 Finite element results showing deformation, interface debonding, and stress profile around the broken tow. (A) Slice view of the deformation around the tow break at the $z=0$ plane, showing break opening and axial stress profile ($\sigma_{11}$). (B) Predicted interface debonding in the form of failed cohesive elements, wherein the failed elements (with damage variable equals one) are shown in red, and the intact elements are shown in blue. (C) Slice views at plane of tow break ($x=0$ plane) with two different in-situ stress states, showing effect of the stress states on the stress concentration profile and factors.
  10. 10 Comparison between FEA and theory on the predicted tensile stress in tow break plane, with in-situ stresses $\sigma^{(k)}_{11,\infty}=1000$, $\sigma^{(k)}_{22,\infty}=-50$, and $\bar{\sigma}^{(k)}_{33}=-50$ (Case A.1 in Table.2). (A) Stress profile by FEA. (B) Predicted stress profile by analytical approach (Eq.30).
  11. 11 Comparison between FEA and theory on the predicted stress concentration magnitude along the positive y-axis and z-axis. Markers are from the element nodal values in FEA, continuous lines are predictions by theory. (Ai) SCF along y-axis covering two adjacent tow widths. (Aii) Magnified view of comparison in (Ai) showing SCF close to the broken tow (about a quarter of tow width). (Bi) SCF along the z-axis (stacking direction). (Bii) Magnified view of data in (Bi) showing SCF close to the broken tow.

###### List of Tables

  1. 1 Mechanical, dimensional, and interface properties of T1100g 12k resin-impregnated tows
  2. 2 Numerical case studies on uni-directional plies

##  Appendix A Derivation of solutions in shear-lag analysis

The steps to obtain solutions are straightforward, and are the same for all
three cases, with Case 3 ($L_{\Leftrightarrow}=L_{\Updownarrow}$) requiring
less algebras. Therefore, only steps for Case 1 are briefly explained here.

As shown in Eq.16, Case 1 means that:

$L_{\Leftrightarrow}>L_{\Updownarrow}\>,\qquad\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}<\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}$
(38)

And from the governing equations, the general solution takes form:

$\delta
u^{(k,0)}_{11}(x)=\left\\{\begin{array}[]{ll}\dfrac{A_{2}+A_{3}}{2A_{1}}\cdot
x^{2}+C_{1}\cdot x+C_{2}\>,&x\in\left[0\>,L_{\Updownarrow}\right)\\\
-A_{2}/A_{5}+C_{3}\cdot\exp\left(\sqrt{A_{5}/A_{1}}\cdot
x\right)+C_{4}\cdot\exp\left(-\sqrt{A_{5}/A_{1}}\cdot
x\right)\>,&x\in\left[L_{\Updownarrow}\>,L_{\Leftrightarrow}\right)\vphantom{\bigg{|}}\\\
C_{5}\cdot\exp\left(\sqrt{(A_{4}+A_{5})/A_{1}}\cdot
x\right)+C_{6}\cdot\exp\left(-\sqrt{(A_{4}+A_{5})/A_{1}}\cdot
x\right)\>,&x\in\left[L_{\Leftrightarrow}\>,+\infty\right)\vphantom{\bigg{|}}\end{array}\right.$
(39)

where $A_{1}$ to $A_{5}$ are non-negative constants involving material
properties and stress states (as defined in Eq.18). The eight boundary
conditions (BCs) are listed in Eq.20 with sub-equations numbered from 1 to 8.
Here denote them as BC-1 to BC-8 for simplicity. Then BC-1 and BC-2 directly
result in $C_{1}=-\sigma^{(k)}_{11,\infty}/E_{l}$ and $C_{5}=0$, whereby
solution is simplified to the form shown in Eq.21.

Substitute BC-8 into BC-4 and BC-6, respectively, then two of the remaining
undetermined constants can be expressed in terms of $L_{\Leftrightarrow}$ as:

$C_{3}=\frac{A_{7}}{2A_{5}K_{s,\Leftrightarrow}}\cdot\exp\left(-\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Leftrightarrow}\right)\,,\quad
C_{4}=\frac{A_{6}}{2A_{5}K_{s,\Leftrightarrow}}\cdot\exp\left(\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Leftrightarrow}\right)$ (40)

where

$A_{6}=A_{2}K_{s,\Leftrightarrow}+\left(A_{5}+\sqrt{A_{4}A_{5}+A_{5}^{2}}\right)\tau^{f}_{s,\Leftrightarrow}\,,\quad
A_{7}=A_{2}K_{s,\Leftrightarrow}+\left(A_{5}-\sqrt{A_{4}A_{5}+A_{5}^{2}}\right)\tau^{f}_{s,\Leftrightarrow}$
(41)

Then substitute Eq.40 into BC-3 and BC-7, the resulting two equations, BC-3*
and BC-7* can be expressed in terms of $L_{\Updownarrow}$ and
$(L_{\Updownarrow}-L_{\Leftrightarrow})$ as

$\displaystyle\text{BC-3*:}\quad\frac{A_{2}+A_{3}}{A_{1}}E_{l}\cdot
L_{\Updownarrow}+\frac{E_{l}\sqrt{A_{5}/A_{1}}}{2A_{5}K_{s,\Leftrightarrow}}\cdot\left[A_{6}\exp\left(-\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)-A_{7}\exp\left(\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)\right]-\sigma^{(k)}_{11,\infty}=0$ (42)
$\displaystyle\text{BC-7*:}\quad\frac{K_{s,\Updownarrow}}{2A_{5}K_{s,\Leftrightarrow}}\left[A_{6}\exp\left(-\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)+A_{7}\exp\left(\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)\right]-\frac{A_{2}}{A_{5}}K_{s,\Updownarrow}-\tau^{f}_{s,\Updownarrow}=0$

where

$L_{\Delta}\equiv L_{\Updownarrow}-L_{\Leftrightarrow}<0$ (43)

as shown before, is the difference between two debonding lengths (negative in
Case 1).

It is observed that components involving $L_{\Delta}$ are exponential terms
whose powers have the same magnitude and opposite signs
($\sqrt{A_{5}/A_{1}}L_{\Delta}$), therefore, the following operations can be
applied to cancel out the terms containing $L_{\Delta}$:

$\displaystyle\left[(\text{BC-7*})\cdot\sqrt{\frac{A_{5}}{A_{1}}}\cdot\frac{E_{l}}{K_{s,\Updownarrow}}+(\text{BC-3*})\right]=0\quad\longrightarrow\quad\exp\left(-\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)=\Psi(L_{\Updownarrow})>0\vphantom{\bigg{|}}$ (44)
$\displaystyle\left[(\text{BC-7*})\cdot\sqrt{\frac{A_{5}}{A_{1}}}\cdot\frac{E_{l}}{K_{s,\Updownarrow}}-(\text{BC-3*})\right]=0\quad\longrightarrow\quad\exp\left(\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)=\Upsilon(L_{\Updownarrow})>0\vphantom{\bigg{|}}$

where $\Psi(L_{\Updownarrow})$ and $\Upsilon(L_{\Updownarrow})$ are functions
of variable $L_{\Updownarrow}$ and other constants, the detailed expressions
are:

$\displaystyle\Psi(L_{\Updownarrow})$
$\displaystyle=-\dfrac{A_{2}+A_{3}}{A_{6}}\sqrt{\dfrac{A_{5}}{A_{1}}}K_{s,\Leftrightarrow}\cdot
L_{\Updownarrow}+\dfrac{K_{s,\Leftrightarrow}}{A_{6}}\left(A_{2}+A_{1}\sqrt{\frac{A_{5}}{A_{1}}}\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}+A_{5}\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)$
(45) $\displaystyle\Upsilon(L_{\Updownarrow})$
$\displaystyle=\dfrac{A_{2}+A_{3}}{A_{7}}\sqrt{\dfrac{A_{5}}{A_{1}}}K_{s,\Leftrightarrow}\cdot
L_{\Updownarrow}+\dfrac{K_{s,\Leftrightarrow}}{A_{7}}\left(A_{2}-A_{1}\sqrt{\frac{A_{5}}{A_{1}}}\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}+A_{5}\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)$

Then the above two relations can be simplified into a quadratic equation of
$L_{\Updownarrow}$ through:

$\Psi(L_{\Updownarrow})\cdot\Upsilon(L_{\Updownarrow})=1\quad\longrightarrow\quad
A_{8}\cdot L_{\Updownarrow}^{2}+A_{9}\cdot L_{\Updownarrow}+A_{10}=0$ (46)

where

$\displaystyle
A_{8}=\left(\frac{E_{l}}{A_{1}}\right)^{2}\cdot(A_{2}+A_{3})^{2}\geq 0\,\qquad
A_{9}=-\frac{2E_{l}}{A_{1}}\sigma^{(k)}_{11,\infty}\cdot(A_{2}+A_{3})\leq
0\vphantom{\bigg{|}}$ (47) $\displaystyle
A_{10}=\left(\sigma^{(k)}_{11,\infty}\right)^{2}-\frac{A_{2}E_{l}^{2}}{A_{1}}\left(\frac{A_{2}}{A_{5}}+2\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)+\frac{E_{l}^{2}}{A_{1}}\left[\frac{A_{6}A_{7}}{A_{5}K_{s,\Leftrightarrow}^{2}}-A_{5}\left(\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right)^{2}\right]$

are constants, and therefore the two solutions for $L_{\Updownarrow}$ are:

$L_{\Updownarrow}^{a}=\dfrac{-A_{9}-\sqrt{A_{9}^{2}-4A_{8}A_{10}}}{2A_{8}}\,,\qquad
L_{\Updownarrow}^{b}=\dfrac{-A_{9}+\sqrt{A_{9}^{2}-4A_{8}A_{10}}}{2A_{8}}$
(48)

Mathematically, in order for the solution pair to be real and positive, the
required conditions are:

$(A_{2}+A_{3})>0\,,\qquad A_{10}>0$ (49)

which directly lead to the formulas given in Eq.24b and Eq.24a, respectively.
In a physical perspective, as explained in Section 2.6, these two conditions
ensure that: 1) there are interface debonding happening, so that the shear-lag
model is valid; 2) there are non-zero frictional shear forces in the debonded
interfaces, so that the broken tow will not result in unstable debonding,
i.e., the system will reach an equilibrium state.

In the next appendix section, it will be shown that the smaller one of the
pair, $L_{\Updownarrow}^{a}$, is the only valid solution, since the other one
is against the prerequisite definition for each case. Then, the remaining
unused boundary condition, BC-5, lead to solution of the last unknown
constant, $C_{2}$. In brief, by solving for the two debonding lengths
($L_{\Updownarrow},L_{\Leftrightarrow}$), all the other undetermined constants
can be solved accordingly.

##  Appendix B Proof of the uniqueness of solutions

Following the derivation of solutions for Case 1 as shown in A, we continue to
show that the solution is uniquely determined. Proofs for the other two cases
are of the same steps and thus they are skipped here.

From the constraints in Eq.44, the relation between $L_{\Updownarrow}$ and
$L_{\Delta}$ can be obtained by:

$\dfrac{\Upsilon(L_{\Updownarrow})}{\Psi(L_{\Updownarrow})}=\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}\right)=\dfrac{A_{6}}{A_{7}}\cdot\dfrac{A_{1}B_{3}+B_{1}L_{\Updownarrow}}{A_{1}B_{2}-B_{1}L_{\Updownarrow}}\equiv\Phi(L_{\Updownarrow})$
(50)

where the three new constants are:

$\displaystyle B_{1}=(A_{2}+A_{3})A_{5}E_{l}K_{s,\Updownarrow}$ (51)
$\displaystyle
B_{2}=A_{2}E_{l}K_{s,\Updownarrow}\sqrt{\frac{A_{5}}{A_{1}}}+A_{5}K_{s,\Updownarrow}\sigma^{(k)}_{11,\infty}+A_{5}\sqrt{\frac{A_{5}}{A_{1}}}\tau^{f}_{s,\Updownarrow}E_{l}\,,$
$\displaystyle
B_{3}=A_{2}E_{l}K_{s,\Updownarrow}\sqrt{\frac{A_{5}}{A_{1}}}-A_{5}K_{s,\Updownarrow}\sigma^{(k)}_{11,\infty}+A_{5}\sqrt{\frac{A_{5}}{A_{1}}}\tau^{f}_{s,\Updownarrow}E_{l}$

Then the two candidate solutions of $L_{\Updownarrow}$, denoted as
$L_{\Updownarrow}^{a},L_{\Updownarrow}^{b}$ in Eq.48, will lead to two
corresponding values of $L_{\Delta}$ as:

$\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{a}\right)=\Phi(L_{\Updownarrow}^{a})\,,\quad\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{b}\right)=\Phi(L_{\Updownarrow}^{b})$ (52)

based on their relation as shown in Eq.50.

Now we prove that the smaller one of the two solutions,
$L^{a}_{\Updownarrow}$, is the only valid one, because the other solution,
$L^{b}_{\Updownarrow}$, will lead to a positive $L_{\Delta}$ which is against
the definition (Eq.22 or 43) for this case. Firstly, it is verified that:

$\dfrac{\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{b}\right)}{\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{a}\right)}-1=\dfrac{\Phi(L^{b}_{\Updownarrow})}{\Phi(L^{a}_{\Updownarrow})}-1=\dfrac{A_{1}B_{1}\left(B_{2}+B_{3}\right)\left(L_{\Updownarrow}^{b}-L_{\Updownarrow}^{a}\right)}{\left(A_{1}B_{3}+B_{1}L_{\Updownarrow}^{a}\right)\left(A_{1}B_{2}-B_{1}L_{\Updownarrow}^{b}\right)}>0$
(53)

so that $L_{\Delta}^{b}>L_{\Delta}^{a}$, because $e^{2\sqrt{A_{5}/A_{1}}\cdot
t}$ is a positive and strictly increasing function for any $t$. To verify
Eq.53, it is readily seen that

$A_{1}>0\,,\quad B_{1}>0\,,\quad B_{2}+B_{3}>0\,,\quad
L_{\Updownarrow}^{b}-L_{\Updownarrow}^{a}=-L_{\Delta}>0$ (54)

so the numerator is positive. Meanwhile, it can be verified that the
denominator is also positive, namely,

$\displaystyle\left(A_{1}B_{3}+B_{1}L_{\Updownarrow}^{a}\right)\left(A_{1}B_{2}-B_{1}L_{\Updownarrow}^{b}\right)$
(55)
$\displaystyle\quad=A_{1}A_{5}E_{l}^{2}\left[\sqrt{\left(\dfrac{A_{2}K_{s,\Updownarrow}+A_{5}\tau^{f}_{s,\Updownarrow}}{K_{s,\Leftrightarrow}}\right)^{2}-A_{6}A_{7}\left(\dfrac{K_{s,\Updownarrow}}{K_{s,\Leftrightarrow}}\right)^{2}}-\left(A_{2}K_{s,\Updownarrow}+A_{5}\tau^{f}_{s,\Updownarrow}\right)\right]^{2}>0$

because both $A_{1}$ and $A_{5}$ are positive, and the remaining are squared
terms. Secondly, we show that

$\left[\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{a}\right)-1\right]\cdot\left[\exp\left(2\sqrt{\frac{A_{5}}{A_{1}}}\cdot
L_{\Delta}^{b}\right)-1\right]=\left[\Phi(L^{a}_{\Updownarrow})-1\right]\cdot\left[\Phi(L^{b}_{\Updownarrow})-1\right]<0$
(56)

so that $L_{\Delta}^{a}$ and $L_{\Delta}^{b}$ have opposite signs, because
$\exp(0)=1$ and $e^{2\sqrt{A_{5}/A_{1}}\cdot t}$ is a positive, strictly
increasing function for any $t$. This way, $L_{\Delta}^{b}$ must be positive
since it is already proved that $L_{\Delta}^{b}>L_{\Delta}^{a}$. However,
$L_{\Delta}$ must be a negative value by its definition for this case,
therefore, the solution of all unknowns are uniquely determined by
$L_{\Updownarrow}^{a}$ and the other candidate solution of inter-ply debonding
length, $L_{\Updownarrow}^{b}$, is invalid. To verify Eq.56, it can be shown
that

$\displaystyle\left[\Phi(L^{a}_{\Updownarrow})-1\right]\cdot\left[\Phi(L^{b}_{\Updownarrow})-1\right]$
(57)
$\displaystyle\quad=\left[\dfrac{4A_{5}K_{s,\Leftrightarrow}\left(2A_{2}K_{s,\Updownarrow}K_{s,\Leftrightarrow}+A_{5}K_{s,\Leftrightarrow}\tau^{f}_{s,\Updownarrow}+A_{5}K_{s,\Updownarrow}\tau^{f}_{s,\Leftrightarrow}\right)}{K_{s,\Updownarrow}\left(A_{2}K_{s,\Leftrightarrow}+(A_{5}-\sqrt{A_{4}A_{5}+A_{5}^{2}})\tau^{f}_{s,\Leftrightarrow}\right)^{2}}\right]\cdot\left[\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}-\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}\right]$

where it is a multiplication of two terms, it is obvious that the first term
is positive, and by the definition of Case 1 as shown in Eq.16

$\text{Case 1:}\quad
L_{\Leftrightarrow}>L_{\Updownarrow}\quad\longrightarrow\quad\tau^{f}_{s,\Leftrightarrow}/K_{s,\Leftrightarrow}<\tau^{f}_{s,\Updownarrow}/K_{s,\Updownarrow}\quad\longrightarrow\quad\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}-\frac{\tau^{f}_{s,\Updownarrow}}{K_{s,\Updownarrow}}<0$
(58)

so that Eq.56 holds true.

Therefore, it is proved that the smaller one of the two solutions for the
inter-debonding length, $L_{\Updownarrow}^{a}$ in Eq.48 is the only valid
solution. From which, $L_{\Delta}$ is uniquely given by Eq.50 so that both
debonding lengths can be determined (Eq.22). Then the remaining undetermined
constants can be expressed in terms of the two debonding lengths as given in
Eq.23.

##  Appendix C Solution for Case 2

As shown in Eq.16, the condition for $L_{\Leftrightarrow}<L_{\Updownarrow}$
is:

$\text{Case
2:}\qquad\tau^{f}_{s,\Leftrightarrow}/K_{s,\Leftrightarrow}>\tau^{f}_{s,\Updownarrow}/K_{s,\Updownarrow}$
(59)

Similar to Case 1 in the main text, by setting $\delta
u^{(k,0)}_{11}(x)\equiv\delta u$ for simplicity, and by utilizing the
previously defined constants (Eq.18), the governing equations for Case 2 can
be expressed in form:

$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{A_{2}+A_{3}}{A_{1}},$
$x\in\rinterval{0}{L_{\Leftrightarrow}}\vphantom{\bigg{|}}$ (60a)
$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{1}{A_{1}}\left(A_{3}+A_{4}\cdot\delta u\right),$
$x\in\rinterval{L_{\Leftrightarrow}}{L_{\Updownarrow}}\vphantom{\bigg{|}}$
(60b) $\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{A_{4}+A_{5}}{A_{1}}\delta u,$
$x\in\rinterval{L_{\Updownarrow}}{+\infty}\vphantom{\bigg{|}}$ (60c)

and eight boundary conditions are comparable to those shown in Eq.20, but the
expression differ due to change of intervals.

After applying the similar procedures, the debonding lengths (for Case 2) can
be solved:

$\displaystyle L_{\Leftrightarrow}$
$\displaystyle=-\dfrac{1}{2A_{8}}\left(A_{9}+\sqrt{A_{9}^{2}-4A_{8}A^{[2]}_{10}}\right)>0\,,\vphantom{\bigg{|}}$
(61) $\displaystyle L_{\Delta}$ $\displaystyle\equiv
L_{\Updownarrow}-L_{\Leftrightarrow}=-\frac{1}{2}\sqrt{A_{1}/A_{4}}\cdot\ln\left(\dfrac{A_{6}^{[2]}}{A_{7}^{[2]}}\cdot\dfrac{A_{1}B_{3}^{[2]}+B_{1}^{[2]}L_{\Leftrightarrow}}{A_{1}B_{2}^{[2]}-B_{1}^{[2]}L_{\Leftrightarrow}}\right)>0\,,\vphantom{\bigg{|}}$
$\displaystyle L_{\Leftrightarrow}$
$\displaystyle=L_{\Updownarrow}-L_{\Delta}<L_{\Updownarrow}\vphantom{\bigg{|}}$

where the constants without the superscript $[2]$ indicates that the
expressions remains identical (case independent), and the updated constants
specific to Case 2 are:

$\displaystyle A_{6}^{[2]}$
$\displaystyle=A_{3}K_{s,\Updownarrow}+\left(A_{4}+\sqrt{A_{4}A_{5}+A_{4}^{2}}\right)\tau^{f}_{s,\Updownarrow}$
(62) $\displaystyle A_{7}^{[2]}$
$\displaystyle=A_{3}K_{s,\Updownarrow}+\left(A_{4}-\sqrt{A_{4}A_{5}+A_{4}^{2}}\right)\tau^{f}_{s,\Updownarrow}$
$\displaystyle A_{10}^{[2]}$
$\displaystyle=\left(\sigma^{(k)}_{11,\infty}\right)^{2}-\frac{A_{3}E_{l}^{2}}{A_{1}}\left(\frac{A_{3}}{A_{4}}+2\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}\right)+\frac{E_{l}^{2}}{A_{1}}\left[\frac{A_{6}^{[2]}A_{7}^{[2]}}{A_{4}K_{s,\Updownarrow}^{2}}-A_{4}\left(\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}\right)^{2}\right]$

and

$\displaystyle B_{1}^{[2]}=(A_{2}+A_{3})A_{4}E_{l}K_{s,\Leftrightarrow}$ (63)
$\displaystyle
B_{2}^{[2]}=A_{3}E_{l}K_{s,\Leftrightarrow}\sqrt{\frac{A_{4}}{A_{1}}}+A_{4}K_{s,\Leftrightarrow}\sigma^{(k)}_{11,\infty}+A_{4}\sqrt{\frac{A_{4}}{A_{1}}}\tau^{f}_{s,\Leftrightarrow}E_{l}\,,$
$\displaystyle
B_{3}^{[2]}=A_{3}E_{l}K_{s,\Leftrightarrow}\sqrt{\frac{A_{4}}{A_{1}}}-A_{4}K_{s,\Leftrightarrow}\sigma^{(k)}_{11,\infty}+A_{4}\sqrt{\frac{A_{4}}{A_{1}}}\tau^{f}_{s,\Leftrightarrow}E_{l}$

##  Appendix D Solution for Case 3

In Case 3, the debonding lengths are equivalent, henceforth the intervals for
the piecewise governing equations of the perturbation axial displacement of
the broken tow reduces to two, and the number of boundary constraints and
undetermined parameters reduce to six.

Set $L_{\Leftrightarrow}=L_{\Updownarrow}=L_{d}$, the governing equations for
Case 3 are:

$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{A_{2}+A_{3}}{A_{1}},$
$x\in\rinterval{0}{L_{d}}\vphantom{\bigg{|}}$ (64a)
$\displaystyle\dfrac{d^{2}}{dx^{2}}\delta u=$
$\displaystyle\frac{A_{4}+A_{5}}{A_{1}}\delta u,$
$x\in\rinterval{L_{d}}{+\infty}\vphantom{\bigg{|}}$ (64b)

The solutions are

$\displaystyle\delta u=$
$\displaystyle\frac{A_{2}+A_{3}}{2A_{1}}x^{2}-\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}x+C_{2},$
$x\in\rinterval{0}{L_{d}}\vphantom{\bigg{|}}$ (65a) $\displaystyle\delta u=$
$\displaystyle\frac{A_{4}+A_{5}}{A_{1}}\delta u,$
$x\in\rinterval{L_{d}}{+\infty}\vphantom{\bigg{|}}$ (65b)

where the debonding lengths are:

$L_{d}=L_{\Leftrightarrow}=L_{\Updownarrow}=\frac{A_{1}}{A_{2}+A_{3}}\left(\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}-\sqrt{\frac{A_{4}+A_{5}}{A_{1}}}\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}\right)$
(66)

and the expression of the unknown constant $C_{2}$ is:

$\displaystyle C_{2}$
$\displaystyle=-\frac{A_{2}+A_{3}}{2A_{1}}L_{d}^{2}+\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}L_{d}+\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}$
(67)
$\displaystyle=\frac{A_{1}}{2(A_{2}+A_{3})}\left(\frac{\sigma^{(k)}_{11,\infty}}{E_{l}}\right)^{2}-\frac{A_{4}+A_{5}}{2(A_{2}+A_{3})}\left(\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}\right)^{2}+\frac{\tau^{f}_{s,\Leftrightarrow}}{K_{s,\Leftrightarrow}}$