Alternating strain regimes for failure propagation in flexural systems

TL;DR

This paper investigates steady-state fracture propagation in a discrete mass-beam structure under sinusoidal loading, combining analytical and numerical methods to identify fracture regimes and dynamic features.

Contribution

It introduces an analytical model reducing the fracture problem to a Wiener-Hopf equation and analyzes steady-state fracture regimes in a discrete structure.

Findings

Minimum wave energy needed for steady-state fracture identified

Dynamic features of the structure during fracture characterized

Transient analysis confirms steady-state regimes

Abstract

We consider both analytical and numerical studies of a steady-state fracture process inside a discrete mass-beam structure, composed of periodically placed masses connected by Euler-Bernoulli beams. A fault inside the structure is assumed to propagate with a constant speed and this occurs as a result of the action of a remote sinusoidal, mechanical load. The established regime of fracture corresponds to the case of an alternating generalised strain regime. The model is reduced to a Wiener-Hopf equation and its solution is presented. We determine the minimum feeding wave energy required for the steady-state fracture process to occur. In addition, we identify the dynamic features of the structure during the steady-state fracture regime. A transient analysis of this problem is also presented, where the existence of steady-state fracture regimes, revealed by the analytical model, are verified and the associated transient features of this process are discussed.

Tables

Table 1. Table 1: Load amplitude values that initiate steady-state failure regimes for ω 0 = 3.1 subscript 𝜔 0 3.1 \omega_{0}=3.1 and r = 0.4 𝑟 0.4 r=0.4 . The first two columns show the values for the average speed and the corresponding minimum load amplitude 𝒫 min subscript 𝒫 min {{\cal P}_{\text{min}}} when the regime was first identified from the numerical scheme of MATLAB. The third column presents the analytical predictions for the load amplitude based on the results of Section  4 and the theory of [ 45 , 46 ] . In the last column, we specify the regime encountered, where “1" represents the alternating generalised strain regime and “2" indicates the pure steady-state regime.

MATLAB results		Analytical results	Regime
$\overline{v}$	${𝒫}_{\text{min}}$	${𝒫}_{\text{min}}$
0.1773	8.985	8.932	1
0.2162	9.376	9.366	2
0.2768	10.045	10.055	1
0.3848	10.715	10.686	2
0.6307	11.887	11.742	1
1.7484	23.048	23.00	2

Equations

$$|w_j(t)|<w_c,j\ge p,p\in \mathbb{Z},$$

$$M\frac{{\mathrm{d}}^2{w}_m(t)}{\mathrm{d}{t}^2}={\mathcal{V}}_m^x(0,t)-{\mathcal{V}}_{m-1}^x(a,t)+H(x-Vt)({\mathcal{V}}_m^{y,+}(0,t)-{\mathcal{V}}_m^{y,-}(a,t))+Q_m(t),$$

$${\mathcal{M}}_m^y(0,t)={\mathcal{M}}_m^y(a,t),$$

$$H(x)=\left\{\begin{array}[]{ll}1\;,&\quad\mbox{ if }x\geq 0\;,\\ 0\;,&\quad\mbox{otherwise}\;.\end{array}\right.$$

$${\mathcal{M}}_m^y(\tilde{x},t)=-\frac{2E_1{I}_1}{a^3}[a^2{\theta }_m^y(t)+a(a-3\tilde{x})({\theta }_m^y(t)+{\theta }_{m+1}^y(t))-3(2\tilde{x}-a)(w_{m+1}(t)-w_m(t))],$$

$${\mathcal{V}}_m^x(\tilde{x},t)=\frac{6E_1{I}_1}{a^3}[a({\theta }_m^y(t)+{\theta }_{m+1}^y(t))-2(w_m(t)-w_{m+1}(t))],$$

$${\mathcal{V}}_m^{y,\pm }(y_{\pm },t)=\mp \frac{12E_2{I}_2}{a^3}{w}_m(t).$$

$$\begin{gathered} 6\left\{2[2w_m(t)-w_{m-1}(t)-w_{m+1}(t)]-a[{\theta }_{m+1}^y(t)-{\theta }_{m-1}^y(t)]\right\} \\ +24r{w}_m(t)H(m-Vt/a)+\frac{M{a}^3}{E_1{I}_1}\frac{{\mathrm{d}}^2{w}_m(t)}{\mathrm{d}{t}^2}=\frac{a^3}{E_1{I}_1}{Q}_m(t) \end{gathered}$$

$$3[w_{m+1}(t)-w_{m-1}(t)]+a[{\theta }_{m+1}^y(t)+{\theta }_{m-1}^y(t)+4{\theta }_m^y(t)]=0,$$

$$\tilde{t}=t\sqrt{\frac{E_1{I}_1}{M{a}^3}},\tilde{v}=V\sqrt{\frac{Ma}{E_1{I}_1}},{\tilde{w}}_m(\tilde{t})=\frac{w_m(t)}{a},{\tilde{Q}}_m(\tilde{t})=\frac{a^2}{E_1{I}_1}{Q}_m(t),$$

$$\begin{gathered} 6\left\{2[2w_m(t)-w_{m-1}(t)-w_{m+1}(t)]-[{\theta }_{m+1}^y(t)-{\theta }_{m-1}^y(t)]\right\} \\ +24r{w}_m(t)H(m-vt)+\frac{{\mathrm{d}}^2{w}_m(t)}{\mathrm{d}{t}^2}=Q_m(t), \end{gathered}$$

$$3[w_{m+1}(t)-w_{m-1}(t)]+[{\theta }_{m+1}^y(t)+{\theta }_{m-1}^y(t)+4{\theta }_m^y(t)]=0.$$

$$w_m(t)=(-1{)}^{m}w(\eta ),{\theta }^y(t)=(-1{)}^m{\theta }^y(\eta ),$$

$$\{w^{\mathrm{F}},{\theta }^{y\mathrm{F}}\}={\int }_{-\infty }^{\infty }\{w(\eta ),{\theta }^y(\eta )\}{e}^{\text{i}k\eta }d\eta ,$$

$$w_{\pm }(k)={\int }_{-\infty }^{\infty }w(\eta )e^{\text{i}k\eta }H(\pm \eta )d\eta ,\pm \text{Im }k>0,w^{\mathrm{F}}=w_{+}+w_{-}.$$

$$\begin{gathered} 6\left\{2[2w(\eta )+w(\eta -1)+w(\eta +1)]+[{\theta }^y(\eta +1)-{\theta }^y(\eta -1)]\right\} \\ +24rw(\eta )H(\eta )+v^2\frac{{\mathrm{d}}^{2}w(\eta )}{\mathrm{d}{\eta }^2}=Q(\eta ) \end{gathered}$$

$$-3[w(\eta +1)-w(\eta -1)]+[4{\theta }^y(\eta )-{\theta }^y(\eta +1)-{\theta }^y(\eta -1)]=0.$$

$${\theta }^{y\mathrm{F}}=-\frac{3\mathrm{i}\sin (k)}{2-\cos (k)}{w}^{\mathrm{F}}$$

$$h_1(k,0+\text{i}kv)w_{+}+h_2(k,0+\text{i}kv)w_{-}=Q^{\mathrm{F}},$$

$$h_j(k,Y)={\Omega }_j^2(k+\pi )+Y^2,j=1,2,$$

$${\Omega }_1(k)=\sqrt{\frac{48{\sin }^4(k/2)}{2+\cos (k)}+24r},{\Omega }_2(k)=\sqrt{\frac{48{\sin }^4(k/2)}{2+\cos (k)}}.$$

$$0+\text{i}kv=\underset{\epsilon \to +0}{\lim }(\epsilon +\mathrm{i}kv)$$

$${\omega }_1(k)={\Omega }_1(k+\pi )=\sqrt{\frac{48{\cos }^4(k/2)}{2-\cos (k)}+24r},$$

$${\omega }_2(k)={\Omega }_2(k+\pi )=\sqrt{\frac{48{\cos }^4(k/2)}{2-\cos (k)}},$$

$$w_{+}+\frac{{\Psi }_{+}(k)}{{\Psi }_{-}(k)}L(k)w_{-}=\frac{Q^{\mathrm{F}}}{h_1(k,0+\mathrm{i}kv)},$$

$${\Psi }_{+}(k)=\frac{\prod _{i=1}^l(0-\mathrm{i}(k-z_{2i}))(0-\mathrm{i}(k+z_{2i}))}{(1-\mathrm{i}k{)}^{2(l-b)}\prod _{j=1}^b(0-\mathrm{i}(k-y_{2j}))(0-\mathrm{i}(k+y_{2j}))},$$

$${\Psi }_{-}(k)=\frac{(1+\mathrm{i}k{)}^{2(l-b)}\prod _{j=0}^b(0+\mathrm{i}(k-y_{2j+1}))(0+\mathrm{i}(k+y_{2j+1}))}{\prod _{i=0}^l(0+\mathrm{i}(k-z_{2i+1}))(0+\mathrm{i}(k+z_{2i+1}))},$$

$$L(k)=\frac{{\Psi }_{-}(k)}{{\Psi }_{+}(k)}\frac{h_2(k,0+\mathrm{i}kv)}{h_1(k,0+\mathrm{i}kv)}.$$

$$\text{Re}(L(k))=\text{Re}(L(-k))\text{ and }\text{Im}(L(k))=\text{Im}(L(-k)).$$

$$L(k)=L_{+}(k)L_{-}(k),L_{\pm }(k)=\text{exp}\left(\pm \frac{1}{2\pi \mathrm{i}}{\int }_{-\infty }^{\infty }\frac{\ln L(\xi )}{\xi -k}d\xi \right),\pm \text{Im}(k)>0.$$

$$\frac{1}{L_{+}(k){\Psi }_{+}(k)}{w}_{+}+\frac{1}{{\Psi }_{-}(k)}{L}_{-}(k)w_{-}=\frac{Q^{\mathrm{F}}}{L_{+}(k){\Psi }_{+}(k)h_1(k,0+\mathrm{i}kv)}.$$

Full text

Keywords: Discrete periodic media, mass-beam structures, fracture, Wiener-Hopf technique, numerical simulations.

1 Introduction

The modelling of periodic flexural materials is a useful tool in understanding the behaviour of structures commonly found in civil engineering, such as buildings, bridges, rooftops, pipeline systems and many more. The need to understand the response of these structures, as shown in Figure 1, is greater when failure initiates and propagates through the system.

We present a simplified analytical and numerical model to represent the failure of a long bridge or a rooftop. The structure is modelled by Euler-Bernoulli beams connecting periodically placed masses, where the failure is assumed to propagate steadily within the structure.

The modelling of failure in structured materials has been widely developed for the analysis of mass-spring systems. These models help in describing several microstructural processes which occur in the fracture of materials [1]. Dynamic Mode III fracture of a square cell lattice, composed of massless springs connecting periodically placed masses, has been considered in [2]. For discrete mass-spring systems which undergo phase transition processes at a uniform rate, see [3]. Dynamic fracture modes in elastic triangular lattices were treated in [4] for a homogeneous lattice and in [5] for an anisotropic lattice. In these models, the micro-level processes involved in the fracture phenomenon can be identified, including effects such as wave radiation occurring when failure propagates. The introduction of inertial links in a lattice brings new features in terms of wave radiation processes that accompany the crack growth, as demonstrated for Mode III fracture inside an inertial square cell lattice in [6].

In addition to the analytical models considered in [2, 3, 4, 6], numerical modelling of Mode I and II crack growth inside strips of triangular lattice has been carried out in [7], where some surprising periodic patterns of fault propagation were observed.

The analytical approach developed in [8] for treating phase transition in a structured material utilises the Fourier transform of the governing equations with respect to a moving coordinate that follows the position of the phase transition front. The problem can then be reduced to a Wiener-Hopf equation along the line containing the defect, where the conditions ahead of and behind the phase transition front differ. This functional equation contains information about the dynamics of the medium on both macro- and micro-scales. The application of the Wiener-Hopf technique to solving dynamic problems for lattices with defects is of wide utility and has also been employed in [9, 10, 11, 12, 13] in analysing scattering and diffraction of waves by defects in several types of periodic media.

Lattices, composed of fundamental mechanical elements, can easily be designed to produce a variety of media capable of controlling the flow of waves for different applications. In particular, the introduction of structural heterogeneities may influence the admissible regimes where steady crack propagation is possible. The effects of inhomogeneities on Mode III crack propagation in two-dimensional lattices have been analysed in [14]. The dynamic failure of dissimilar chains has been analysed in [15, 16].

Periodically distributed inhomogeneities in mass-spring lattices have been shown to enhance fracture propagation in [17]. There, a structured interface was shown to induce a failure mechanism caused by a “knife wave” or localised deformation surrounding the crack, capable of sustaining crack propagation within the interface. Numerical simulations of a crack growing through a high-contrast interface were later carried out in [18]. Moreover, these numerical simulations revealed non-steady failure regimes such as clustering, where a crack may propagate in an unusual way.

Other failure processes propagating in a lattice have been treated in [19], where the extraction of a mass-spring chain due to a point force was considered. The study of a propagating bridge crack in an inhomogeneous lattice has been carried out in [20]. Different fracture criteria, independent of or dependent on time, can also produce a variety of fracture patterns within lattice systems and can affect the existence of admissible steady-state fracture regimes. The impact of time dependent fracture criteria on the failure of one-dimensional mass-spring chains has been investigated in [21].

The important question of admissibility of unstable and stable crack regimes for cracks propagating in discrete periodic media has been discussed in [1, 8, 22]. The existence of admissible low velocity steady-state failure regimes in a mass-spring square cell lattice [8] has recently been re-addressed in [23]. The admissibility of crack propagation regimes can also be investigated for different loading conditions. The influence of a moving load on Mode III fracture regimes in a discrete mass-spring structure has been analytically and numerically studied in [24]. A transient analysis of thermal shock induced crack propagation in triangular lattices has also been carried out in [25, 26] and in [27], where cracks with chiral coatings were considered.

Failure regimes in discrete structures with non-local interactions and their effects on the wave radiation processes are presented in [28] and for media with additional non-linearities in the connections see [29]. Macro and micro-level models for the reversible decohesion of finite elastic layers, associated with the denaturation of DNA, are considered in [30].

In contrast to the study of dynamic fracture in mass-spring systems, very few articles concentrate on the failure of systems incorporating flexural elements such as beams. The study of bending modes in 2D discrete flexural structures, composed of beams connecting masses and containing static faults, can be found in [31]. The analysis of a static Mode III crack in a three-dimensional beam lattice, representing an open cell foam, is presented in [32]. The investigation of the fracture toughness of structured materials, with applications to understanding the response of foams with cracks and broken elements, is considered in [33]. The dynamics of periodic materials with small suddenly appearing flaws have been analysed in [34]. Quasi-static damage propagation in two-dimensional beam structures under tensile loading has been considered in [35], motivating the design of fault-tolerant beam lattices in [36].

One-dimensional mass-beam chains supported by an elastic foundation, have been used to create a simplified model of a collapsing bridge [37]. This model has been further developed in [38] and used to predict the collapse rate of the San-Saba rail road bridge, Texas, in 2013. Transition waves in continuous flexural systems have been shown to exhibit surprising non-steady failure behaviour [39].

Experimentally, one can observe counter-intuitive fracture behaviour in elastic materials such as rubber sheets [40], where analytical descriptions and numerical simulations of this process have been given in [41]. The influence of lattice vibrations and wave radiation processes on the stability of crack trajectories in Silicon crystals has been investigated experimentally and numerically using molecular dynamics simulations in [42, 43, 44].

In [45], an analytical model for the failure of a discrete flexural structure, contained within an interface, was studied. The fracture was assumed to be caused by a remote sinusoidal load and the dynamic features of the structure were studied. Further, in [46] the transient failure process was numerically modelled and the results of [45] were verified.

Here, we consider failure propagation within the mass-beam structure shown in Figure 2. The failure is represented by the sequential removal of transverse links supporting the masses along the central chain, aligned with the horizontal axis. In [45, 46], the failure process was investigated under the assumption that the supporting columns of the central chain are removed if the masses achieve a positive critical displacement. The steady-state regimes achieved in this case are referred to here as “pure steady-state regimes” and assumed to be driven by a remote mechanical sinusoidal load. In terms of applications, it is realistic to consider the case when the supporting links break if a positive or negative critical displacement is reached. As we show here, along with the pure steady-state regimes observed in [45, 46], the imposed fracture criterion yields additional regimes, which can be predicted and verified numerically. These regimes correspond to the case when, at each consecutive failure in the system, a change in the sign of the bending moments and shear forces in the supporting beams can be observed. These regimes are investigated in detail here and are referred to as “alternating generalised strain regimes”. We note that for mass-spring systems, the analogous phenomenon is studied in [17].

The structure of the present article is as follows. In Section 2, we present an analytical description of the problem concerning failure propagation in the mass-beam structure subjected to a sinusoidal load. This section also includes the reduction of the problem to a functional equation of the Wiener-Hopf type, from which the alternating generalised strain regimes can be identified. The characterisation of the dispersive nature of this particular discrete system is given in Section 3. The results of Section 3 are then used in Section 4 to solve the Wiener-Hopf equation. We also identify the dynamic properties of the system during the steady propagation of the alternating generalised strain regime in Section 5. In Section 6, we present numerical simulations which support the analytical results of Sections 2–5. We give some conclusions in Section 7. Finally, Appendices A and B contain details of some derivations.

2 Model of failure within a discrete periodic flexural structure

2.1 Description of the problem

We consider a structure composed of a mass-beam chain, as shown in Figure 2(a). This structure is formed from periodically placed point masses, connected longitudinally (along $x$-axis) by massless Euler-Bernoulli beams all with Young’s modulus $E_{1}$ and second moment of area $I_{1}$. Each junction node has mass $M$ and corresponds to an index $m\in\mathbb{Z}$. This chain is assumed to be partially supported by transverse Euler-Bernoulli beams (parallel to the $y$-axis), having Young’s modulus $E_{2}$ and second moment of area $I_{2}$. These beams connect the masses to an interface where the beams are clamped. All beams have length $a$.

Inside the structure, failure is assumed to propagate with a uniform speed $V$ as result of the breakage of the transverse connections to the interfaces within the structure. This breakage occurs when the absolute value of the displacement $w_{p}$ of the $p^{\textmd{th}}$ mass (where $m=p$ corresponds to the position of the transition front in the structure, see Figure 2(a), that depends on time) reaches a critical value $w_{c}$. Thus we assume for bonds to remain intact

[TABLE]

and when the condition $w_{p}=\pm w_{c}$ is fulfilled, the transverse links at the $p^{th}$-mass break and the moving interface advances a distance $a$ in the structure to the $(p+1)^{\textmd{th}}$ node.

The condition (1) enables us to exclude the case of non-steady propagation of the interface inside the structure. Under these assumptions, we consider the case when the generalised strains (moments and shear forces) inside the transverse links alternate in sign during the failure process.

2.2 Governing equations

At a given time $t$, the position of the interface in the structure is given by $m=\lfloor Vt/a\rfloor$, where $\lfloor x\rfloor$ denotes the integer part of $x$. Here, the inequality $m\geq\lfloor Vt/a\rfloor$ corresponds to masses located in the supported part of the structure and $m<\lfloor Vt/a\rfloor$ are those in the unsupported region of the structure (see Figure 2).

The equations for the balance of linear and angular momentum for the $m^{th}$ mass, $m\in\mathbb{Z}$, are

[TABLE]

[TABLE]

where $H$ is the Heaviside function

[TABLE]

In (2) and (3), ${\cal M}^{y}_{m}(\tilde{x},t)$ and ${\cal V}^{x}_{m}(\tilde{x},t)$, denote the $y$-component of the moment vector and the shear forces in the $m^{th}$ horizontal beam, respectively, at time $t$ and position $\tilde{x}$ along the beam, where $0<\tilde{x}<a$, (see Figure 2(b)). They can be written in terms of the generalised coordinates describing the motion of the $m^{th}$ mass as

[TABLE]

[TABLE]

where $w_{m}$ is the displacement and $\theta_{m}^{y}$ is the rotation about the $y$-axis of the $m^{th}$ mass (see Appendix A for the derivation of (4) and (5)). In a similar way, the shear forces ${\cal V}^{y,+}_{j}(y_{+},t)$ and ${\cal V}^{y,-}_{j}(y_{-},t)$ correspond to those in the transverse beam located at $0<y<a$ and $-a<y<0$, respectively, for $x=ja$ where $j\geq\lfloor Vt/a\rfloor$. Here, the notation $\tilde{y}_{+}=y$ ($\tilde{y}_{-}=y+a$) represents the local coordinate for the transverse beam located at $y>0$ ($y<0$). The shear forces ${\cal V}^{y,-}_{m}$ adopt the same convention as described for ${\cal V}^{y,+}_{m}$ in Fig. 2. As discussed in Appendix A, the internal shear forces in these transverse elements take the form

[TABLE]

In (2), the term $Q_{m}(t)$ is the applied load at a node $m$. In addition, in Appendix A, we show the angular momentum balance about the $x$-axis for a mass in the supported region shows that the masses do not rotate about about the $x$-axis, i.e. $\theta^{x}_{m}=0$ for $m\in\mathbb{Z}$.

Equations (2) and (3) then become

[TABLE]

and

[TABLE]

where the contrast parameter $r=E_{2}I_{2}/E_{1}I_{1}$. In (7) and (8), we introduce the normalisations

[TABLE]

where the symbol “tilde" will be omitted in the following for ease of notation. In this case, the dimensionless governing equations for the system are then

[TABLE]

and

[TABLE]

2.3 Derivation of the Wiener-Hopf equation

Since we consider the case of alternating generalised strains, we look for the displacements and the rotations as functions of the moving coordinate $\eta=m-vt$ as follows

[TABLE]

where $w(\eta)$ and $\theta(\eta)$ play the role of envelope functions (see [17]). In addition, we assume the load takes the form $Q_{m}(t)=(-1)^{m}Q(\eta)$.

We introduce the Fourier transforms with respect to the variable $\eta$ for the quantities $w^{{\rm F}}$, $\theta^{y{\rm F}}$ as

[TABLE]

where $k$ is the dimensionless wavenumber. The following “half" transforms are also used:

[TABLE]

Here $w_{+}$ ($w_{-}$) corresponds to a function analytic in the upper (lower) half of the complex plane defined by $k$.

After substitution of (11) into (9) and (10) we obtain

[TABLE]

and

[TABLE]

The Fourier transform of (13) with respect to the moving coordinate $\eta$ leads to

[TABLE]

and together with the Fourier transform of (12) we then obtain

[TABLE]

where $Q^{\rm F}$ is the Fourier transform of the load $Q$ and

[TABLE]

Here

[TABLE]

in (15) appears as a result of the causality principle discussed in [8], and represents the passing from the transient regime to the steady-state regime in the Laplace transform. Note that the expressions (17) correspond to dispersion relations for the problem in [45, 46], where the fracture was considered for the pure steady-state failure regimes.

3 Dispersion relations and characterisation of waves in the structure

The dispersive properties of the structure can be determined from the zeros of $h_{j}$, $j=1,2,$ in (15). In replacing $Y$ by ${\rm i}\omega$ in (16) and rearranging for $\omega$ we obtain the dispersion relations

[TABLE]

for the waves ahead of the transition front and

[TABLE]

for waves behind the transition front.

Waves inside the structure are assumed to propagate as a result of either the external action from a remote load or through dissipation from the transition front. These waves can be identified from the functions in (16). In Figure 3, we give a representative example showing the dispersion curves $\omega_{1}$ and $\omega_{2}$ as functions of the wave number $k$, based on (19)–(20), along with the generic line $\omega=kv$.

Characterisation of the waves. Intersections of the line defined by $\omega=kv$ with the curve based on the relation $\omega_{1}$ in (19) represent waves propagating ahead of the transition front ($\eta\geq 0$). Those intersection points of this line with $\omega_{2}$, given by (20), correspond to waves behind the transition front ($\eta<0$).

By computing the group velocity of the wave defined by $v_{g}=d\omega/dk$ at each intersection point, it is possible to determine the direction of propagation of the wave. Any intersection point of the curve given by $\omega_{1}$ with the line $\omega=kv$ where $v_{g}<v$ or $v_{g}>v$, indicates a wave propagating in the supported region towards or away from the transition front, respectively. In a similar way, any intersection of the curve provided by $\omega_{2}$ with the ray $\omega=kv$ where $v_{g}>v$ or $v_{g}<v$ indicates a wave propagating in the unsupported region towards or away from the transition point, respectively.

Roots of the functions $h_{j}$, $j=1,2$. Here we describe, in general, possible roots of the functions $h_{j}$, $j=1,2$, which will determine the nature of the waves inside the structure. These roots are associated with the intersections of the ray $\omega=kv$ and the curves determined by the functions $\omega_{j}$, $j=1,2$.

For the analytical model presented here, we consider the speed range $0<v\leq{v_{1}=1.117}$. Here the upper bound determines the maximum speed for which the failure process can occur under the alternating generalised strain regime. This upper bound is the slope of the dotted inclined line in Figure 3. For the failure process with speed $v$ to occur, it is required that there exists one intersection point of the line $\omega=kv$ with the curve $\omega_{2}$ where $v_{g}>v$. Such a point defines a wave that can propagate behind and deliver energy to the transition front.

Zeros of $h_{1}$. The function $h_{1}(k,{\rm i}kv)$ has one, three or more pairs of simple zeros at the points $k=\pm y_{1},\pm y_{2},\dots,\pm y_{2b+1}$, with $b\in\mathbb{Z}$, $b\geq 0$. Here $b$ depends on the speed $v$. Prior to taking the limit in (18), according to [8], these zeros are located in the complex plane and possess a small imaginary part. The sign of the imaginary part is then determined by comparing the group velocity $v_{g}$ with the failure speed $v$. For the points

with $k=\pm y_{1},\dots,\pm y_{2b+1}$, the group velocity $v_{g}<v$ and 2. 2.

with $k=\pm y_{2},\dots,\pm y_{2b}$, the group velocity $v_{g}>v$.

If $v_{g}<v$ ($v_{g}>v$), these points are located in the upper (lower) part of the complex plane defined by $k$. Thus in the limit (18), we have $k=\pm y_{1}+{\rm i}0,\dots,\pm y_{2b+1}+{\rm i}0$ and $k=\pm y_{2}-{\rm i}0,\dots,\pm y_{2b}-{\rm i}0$.

Here, the points with $v_{g}<v$ can represent waves produced by a load situated far ahead of the transition front. In the loading problem considered in the next section, the points with $v_{g}>v$ will be associated with transmitted waves that propagate ahead of the transition front in the steady-state failure process. These waves can appear for particular failure speeds and an example of an interval of such speeds is shown in the shaded region in Figure 3. In this case, the boundaries of this region correspond to the failure speed $v={0.5214}$ and $v={0.703}$.

Zeros of $h_{2}$. The function $h_{2}(k,{\rm i}kv)$ can have three, five or more pairs of simple zeros at $k=\pm z_{1},\pm z_{2},\pm z_{3},\dots,\pm z_{2l+1}$, $l\in\mathbb{Z}$, $l\geq 1$. Again $l$ is a parameter which depends on the speed $v$. For

$k=\pm z_{1},\pm z_{3},\dots,\pm z_{2l+1}$, we have $v_{g}<v$, 2. 2.

$k=\pm z_{2},\pm z_{4},\dots,\pm z_{2l}$, it holds that $v_{g}>v$.

Again, in the limit in (18), we receive $k=z_{1}+{\rm i}0,\dots,z_{2l+1}+{\rm i}0$ and $k=\pm z_{2}-{\rm i}0,\dots,\pm z_{2l}-{\rm i}0$.

In the next section, in the problem concerning the loading of the structure, the wave numbers here with a negative or positive imaginary part will correspond to feeding or reflected waves, respectively, that appear behind the transition front during the steady failure regime.

4 The Wiener-Hopf equation and fracture criterion

Here, we develop the solution of (15) for the alternating generalised strain regime.

We rewrite this equation in the form

[TABLE]

with

[TABLE]

where $\prod_{j=1}^{0}=1$ and

[TABLE]

Here, $L(k)>0$ for $k\in\mathbb{R}$ and satisfies the conditions required for its factorisation. That is, for $k\in\mathbb{R}$

[TABLE]

In addition, $L(k)\to 1$ as $k\to\pm\infty$ and the index of the function $L(k)$ is zero. Then, $L(k)$ can be written using the Cauchy-type factorisation in the form:

[TABLE]

In this representation, the function $L_{+}$ ($L_{-}$) is analytic in the upper (lower) half of the complex plane defined by $k$. Then, (21) can be written in the form of the Wiener-Hopf equation

[TABLE]

4.1 The solution of the Wiener-Hopf equation

As in [8], non-trivial solutions of (26) correspond to singular points of the right-hand side. Such singular points occur when

[TABLE]

As $L(k)$ has no real zeros and, referring to Section 3 and (22), we note the above left-hand side is only zero when $k=\pm y_{2j+1}-\mathrm{i}0$, $0\leq j\leq b$, and $k=\pm z_{2j}-\mathrm{i}0$, $1\leq j\leq l$. The points $k=\pm y_{2j+1}-\mathrm{i}0$, $0\leq j\leq b$, correspond to the action of a remote force ahead of the transition front as they are associated with the inequality $v_{g}<v$. On the other hand, the points $k=\pm z_{2j}-\mathrm{i}0$, $1\leq j\leq l$, represent the loading from behind the transition front (with $v_{g}>v$).

Here, we assume that the loading of the structure takes the form of an oscillatory mechanical load. This load is situated at some point far from the transition front on the negative $x$-axis with frequency $\omega_{0}=vk_{f}$, where we note in addition $\omega_{0}=\omega_{1}(k_{f})$. Here $k_{f}$, corresponds to a feeding wave that causes the failure to propagate with constant speed $v$ (see section 6.2 for the description of how $k_{f}$ is chosen). Following [8], this allows one to rewrite the right-hand side of (26) as

[TABLE]

where $\phi$ is the phase of the feeding wave and $C$ is a constant to be linked to the amplitude of the feeding wave. The solution to the above equation is then found in the form

[TABLE]

4.2 The fracture criterion and uniqueness of the solution

In terms of the moving coordinate system, the condition (1) can be re-interpreted for the amplitude function $w(\eta)$ as

[TABLE]

The first of these conditions can be satisfied by computing the limits

[TABLE]

using (28) and (29). In doing this we obtain

[TABLE]

which allows for the determination of $\phi$ if the feeding wave amplitude, and consequently $C$, are known. With regard to the second condition in (30), we construct the asymptote of $w_{+}$ in (28) as $k\to\mathrm{i}\infty$. which in terms of the physical variables represents the expansion

[TABLE]

We have

[TABLE]

for $k\to{\rm i}\infty$. Owing to the fact that

[TABLE]

we obtain that the second condition in (30) is satisfied if

[TABLE]

This allows for $\phi$ to be uniquely determined as a solution of (32).

5 Dynamic loading and waves in the steady-state fracture regimes

We discuss the dynamic behaviour of the structure, subjected to a sinusoidal load, during the alternating generalised strain failure process. In particular, concerning the wave radiation processes associated with the moving transition front, we discuss the reflected and transmitted waves emitted from this point.

5.1 Sinusoidal remote load and the feeding wave

We assume a sinusoidal load acts at some point far from the origin on the negative $x$-axis (see Figure 2) inside the structure. The load has amplitude $P_{0}$ and frequency $\omega_{0}$. Following the normalisations adopted in Section 2.3, we introduce the normalised forms of these quantities as

[TABLE]

From here on, we omit the symbol $``\,\tilde{\,\,}\,"$. This load generates a wave incident on the moving interface with the form

[TABLE]

where $\phi$ is the feeding wave phase, $\eta=m-\omega_{0}t/k_{f}$, with $\omega_{0}/k_{f}=v$ being the phase speed of the outgoing wave and

[TABLE]

(see the Appendix B for the derivation of the amplitude in (34)).

We now determine the expression for the constant $C$ in the solution (28), (29) by considering the fracture criterion at the point $\eta=0$. We note that the form of the feeding wave based on (34) can also be found if one considers the complex residues of the simple poles at $k=\pm k_{f}$ for the inverse Fourier transform of $w_{-}$. This is derived under the assumption $\eta\to-\infty$, corresponding to a remote distance far behind the transition front and using

[TABLE]

In doing so, we have

[TABLE]

Noting that

[TABLE]

we obtain

[TABLE]

Comparing with (34), we see that

[TABLE]

5.2 Criterion for steady-state failure propagation in the alternating generalised strain regime

For the crack to propagate steadily, it is required that amplitude of the feeding wave reaches or exceeds the critical displacement, i.e.

[TABLE]

Using this condition, (32) and (39) yields the criterion

[TABLE]

for the failure to propagate steadily in the alternating generalised strain regime, where

[TABLE]

Figure 4, shows the plot of ${\cal P}_{\text{min}}$ in (41) for a representative range of load frequencies. The computations are performed for $r=0.4$, corresponding to a structure with soft supports, which will be investigated further in Section 6. For a given choice of the load frequency, the grey solid and dashed lines, based on (41) in Figure 4, predict the ${\cal P}_{\text{min}}$ required for a fracture process with speed $v=\omega/k_{f}$ to propagate through the flexural system. Here $k_{f}$ is determined from the dispersion diagram. The grey solid curve corresponds to regimes with a higher fracture speed than those connected with the grey dashed curve. These curves are associated with the ${\cal P}_{\text{min}}$ for the alternating generalised strain regime.

For illustrative purposes, we also plot the results of [45, 46] with black curves and these predict the appearance of pure steady-state modes, where the sign of the bending moments and shear forces are uniform at each instant of the failure process.

Accompanying these analytical computations is a verification of the predictions based on numerical simulations for a sufficiently large finite structure implemented in MATLAB (and described in detail in Section 6). There is a excellent agreement between the analytical results, which are based on an infinitely long medium, and the results of the MATLAB computations based on a finite medium.

5.3 Dynamic features for the alternating generalised strain regime

Here, we give the representation of the envelope functions for the reflected and transmitted waves encountered during the failure process based on (28) and (29). The expressions for the associated waveforms in the structure may then be determined by combining the results presented here with (11).

Reflected waves behind the transition front. The reflected waves can be derived from the residues of the poles $k=\pm z_{2j+1}+\mathrm{i}0$, $1\leq j\leq l$, of the function $w_{-}$ in (29). The functions $w_{r}^{(s)}$, $0\leq s\leq l$, representing the reflected waves, have the following form

[TABLE]

with amplitude $A_{r}^{(s)}$ given by

[TABLE]

where

[TABLE]

and the phase shift $\psi_{r}^{(s)}$ given as

[TABLE]

Waves transmitted ahead of the transition front. The waves that are transmitted ahead of the propagating interface are associated with the poles of the function $w_{+}$ in (28), at the points $k=y_{2j}-\mathrm{i}0$, $1\leq j\leq b$.

The transmitted waves $w_{tr}^{(s)}$, $1\leq s\leq b$, are given by:

[TABLE]

Here the amplitude $A_{tr}^{(s)}$ has the form

[TABLE]

with

[TABLE]

The phase shift $\psi_{tr}^{(s)}$ is

[TABLE]

Note, if $b=0$, then there are no propagating waves transmitted to the intact structure.

6 Numerical simulations modelling the failure process

In this section, we implement a numerical scheme in MATLAB that models the behaviour of a finite flexural structure subjected to a sinusoidal load. In addition, using the numerical scheme we also trace the propagation of failure inside this structure as a result of the action of the load. The scheme is based on the normalised governing equations (9) and (10) for the system, which are solved using the ode45 routine of MATLAB.

We use the results from the analysis of the alternating generalised strain fracture regime presented in Section 2, to show that we can predict: (i) when such regimes occur, (ii) the speed with which the fracture propagates in such regimes and (iii) the behaviour of the structure in these particular regimes.

6.1 The numerical model

A finite structure composed of 3800 nodes, connected by massless beams, is considered. As the numerical scheme is based on the dimensionless equations (9) and (10), the dimensionless length and flexural stiffness of the connecting beams can be taken equal to unity. Therefore, the position of each mass is then given by $x=n$, $1\leq n\leq 3800$. The system is initially at rest. We consider this structure composed of a region supported by transverse links and a region where the masses along the central axis of the structure are unsupported by such links. The supported region is characterised by the heterogeneity parameter $r$. The nodes corresponding to $0\leq x\leq 1999$ constitute the unsupported part of the structure, behind the interface, and the remaining 1800 nodes (corresponding to $2000\leq x\leq 3800$) represent the system in the supported region. Therefore, the interface between these two media is initially at $n=2000$.

The external applied force is taken as a sinusoidal force having the form $P_{0}\sin({\omega}_{0}t)$, where $P_{0}$ and ${\omega}_{0}$ are the dimensionless load amplitude and frequency. This load is situated at node $x=1800$. We can choose the value of ${\omega}_{0}$ based on the frequency $\omega$ of the unforced problem for the unsupported structure of infinite extent in the $x$ direction. To generate feeding waves in the unsupported structure capable of reaching the transition front, the load frequency $\omega_{0}$ should be chosen in the interval $0<\omega_{0}<\sqrt{48}$, representing the passband for the analogous infinite medium. Below, we consider $\omega_{0}=3.1$ and 5.9 representing frequencies located in the middle and close to the upper limit of the passband, producing waves capable of exciting the microstructure. In the numerical scheme, the right end of the structure is clamped at $x=3800$, whereas the node at $x=1$ is free. The structure’s size has been chosen to minimise effects due to reflections produced by the ends.

As discussed in Section 5.1, the feeding waves generated by the load help to initiate and propagate the failure process. Fracture regimes in two different types of structure are examined here:

Case 1 : A structure whose supporting transverse links are “softer” than those along the central axis of the structure ($r=0.4$)

Case 2 : A structure with transverse supports that are much stiffer than the beams along its central axis ($r=3.4$).

6.2 Case 1 - soft supports

Here we show that it is possible to encounter two fracture regimes. The first involves the regime where the failure will propagate if (1) is violated at the transition front. The second concerns the regime defined by

[TABLE]

where $p$ represents the node number where the interface is located at a given time $t$. This was studied in detail in [45, 46]. For these regimes, the sign of the generalised strains at the instant of each fracture in the structure remains constant. We refer to the regimes associated with (42) as the pure steady-state regimes.

We first show that speeds associated with the failure regimes observed can be predicted from the dispersion curves (19), (20) and those identified in [45, 46].

Failure speed predictions. A list of steady-state speeds for the considered fracture regimes can be determined using the dispersion curves discussed in Section 3. Figure 5(a) shows the dispersion curves (19), (20) plotted as functions of the wave number $k$ for the case $r=0.4$. In addition, the horizontal dashed lines represent the loading frequencies $\omega_{0}$ chosen for the MATLAB computations.

For the case $\omega_{0}=3.1$, the horizontal dashed line in Figures 5(a) and (b) associated with this frequency intersect the curves given by $\omega_{2}(k)$ and $\omega_{4}(k)$ infinitely many times. Since we require a feeding wave to reach the transition front inside the structure, we need the group velocity $v_{g}$ of this wave to be greater than the fracture speed $v$ i.e. the slope of the line $\omega=kv$. For $k>0$, we have indicated using squares the intersection points, representing possible feeding waves, where this criterion is satisfied in Figure 5(a) and (b).

Without loss of generality, we take the first of these points $(k^{*},\omega)$, shown in Figure 5(a), along the line $\omega=3.1$ with the smallest wavenumber. We connect this point to the origin by a ray whose slope is defined by $v=\omega/k^{*}$ (see Figure 5(a)). This slope predicts a possible steady-state speed of failure inside the structure. In the case of $\omega=3.1$ in Figure 5(a), such an intersection point is given by $(k^{*},\omega)=(4.915,3.1)$ and the slope of the corresponding ray is $0.6307$. That particular intersection point corresponds to an alternating generalised strain failure regime with speed $v=0.6307$ initiated by the sinusoidal load of frequency $\omega=3.1$.

Owing to the $2\pi$-periodicity of the function $\omega_{2}(k)$, a decreasing sequence of possible speeds can be deduced in the form

[TABLE]

As an example, for the case $\omega=3.1$ in Figure 5(a), this list takes the form

[TABLE]

The intersection points connected with these speeds are shown along the line defined by $\omega=3.1$ in Figure 5(a) by squares.

In the MATLAB simulations we can also identify pure steady-state regimes corresponding to the case when the failure propagation occurs as a result of the violation of (42) at the transition front. In this scenario, following [45, 46], the dispersion relations are given by  $\omega_{j+2}(k)=\Omega_{j}(k)$, $j=1,2,$ (see (16)) and are shown in Figure 5(b). Repeating the procedure outlined above, we can obtain the speeds in the pure steady-state regime. For $\omega=3.1$, the list of predicted speeds for this regime is:

[TABLE]

Predictions for the potential steady-state speeds at different load frequencies can be calculated in a similar way.

Results of the transient analysis and identification of failure regimes. Here, for $r=0.4$ and $\omega=3.1$, we analyse the failure processes in the finite structure for various load amplitudes. We show that the failure processes discussed here and in [45, 46] both can appear as a result of the action of the load.

Figure 6 shows a result of the transient computation presented here when one of the alternating generalised strain failure regimes is initiated in the finite medium. As mentioned earlier, for a given load frequency, the appearance of these regimes is dependent on the applied load amplitude. In Figure 6, we show the dynamic response of the structure when the speed of this failure regime is maximal. There, as predicted by the analytical results in Sections 2–5, the transition front propagates with a constant speed as a result of feeding waves from the load providing energy to the front. This allows the front to reach the critical displacement for the failure of the supports. During the steady propagation, waves are radiated outward from the transition front when the failure of the supporting beams is achieved. The video of the Supplementary data shows the evolution of this failure process and the behaviour of the structure resulting from the transient analysis undertaken here. Later we show that in the transient regime the failure process can occur non-steadily before reaching such states. In particular, we will also demonstrate that dynamic response of the structured medium is different when undergoing an alternating generalised strain failure mode when compared with a pure steady failure mode.

Transient behaviour of the transition front** **

The failure of the finite structure resulting from the harmonic load was performed for a large interval of load amplitudes. We denote by $m^{*}$ the index of the mass where the transition front is initially located. For the simulations $m^{*}=2000$. The index $m_{i}=m-m^{*}+1$ refers to the index of a mass in the supported region where the failure of the transverse beams at time $t_{i}$, $i\geq 1$, $i\in\mathbb{Z}$, occurs. Simultaneously, as the length of the beams along the central axis is equal to unity, $m_{i}$ also represents the position of the failure in the supported region measured with respect to the initial position of the transition front.

Figure 7(a) shows several examples of the failure position as a function of time $t$ for various load amplitudes where steady failure regimes can be observed in the simulations. Since the results are based on a transient simulation, we can expect some effects in the initial period of the failure process. For instance, for ${\cal P}=11.496$ in Figure 7(b) we see a visible oscillatory behaviour in the transition front speed initially. Similar behaviour can also be seen for the computations based on the other load amplitudes.

After these initial periods, the transition process settles and the position of failure appears to behave as a linear function of time, representing that the steady-state failure phenomenon has been achieved. After the transient failure process, we can calculate the average failure speed from the data shown. The average speed $\bar{v}$ is then computed as the slope of the line of best fit for this data. In this case, the average speeds obtained are $\bar{v}=0.1773,0.2162,0.2768,0.3848,0.6307,1.7484$, that agree with the analytical predictions for the possible failure speeds of Section 6.2. As expected from physical considerations, Figure 7 shows that as the load amplitude is increased, a higher failure speed is achieved.

Actually, one can observe further oscillations in the failure speed when it appears to have settled to a uniform state. We define the instantaneous speed $v_{i}$ for the fracture process as

[TABLE]

Using the notion of the instantaneous speed, we also can understand how quickly the failure process settles to the steady-state limit.

In Figure 8(a), we show the instantaneous speed plotted as function of the failure position $m_{i}$. Here, the instantaneous speed data presented corresponds to cases considered in Figure 7. To demonstrate the speed of convergence to the steady-state failure process, each set of data for the instantaneous front speed has been normalised by the analytical prediction for the speed ${v_{\text{pred}}}$ that was attained in Figure 7 (see Section 6.2). In this case, convergence to the steady-state speed is represented by the convergence to unity in Figure 8. It is observed that the instantaneous speeds, for nearly all the load amplitudes, oscillate about the predicted speed in each case (see the magnification in Figure 8(b)). In fact, Figure 8(b) shows that there exist micro-oscillations in speed of the failure front not necessarily seen in Figures 7(a) and (b). In addition, these oscillations in the transition front speed are not part of the analytical model studied in Section 2–5, where the speed of the front is assumed to be uniform.

Larger oscillations in the transition front speed are attained with larger values of the amplitude $\cal P$. The speed of convergence of these regimes appears to be a non-monotonic function of the load amplitude.

In Figure 8(c), we show a special case (${\cal P}=10.101$) where the instantaneous speed of the front does not converge to unity in the steady failure process, but instead settles to two distinct speeds (given by $v_{i}/v_{\rm pred}$ approximately equal to 0.79 and 1.39 after the 118th failure event). Although there is an apparent jump in the instantaneous speed of the failure, the average speed $\bar{v}$ is in fact equal to the predicted value 0.2768. Moreover, the corresponding data representing the failure position as a function of time in Figure 7(a) exhibits a straight line that is commonly associated with the steady-state failure propagation.

Dependency of the average failure speed on the load amplitude** **

In the Figure 9, we show the average fracture speed $\bar{v}$ as a function of the normalised oscillating force amplitude ${\cal P}$. There we see six plateaus corresponding to the steady-state speeds reached by the transition front observed in Figures 7 and 8. These include the speeds $\bar{v}=0.1773,0.2768,0.6307$, which correspond to alternating generalised strain regimes, and $\bar{v}=0.2162,0.3848,1.7484$ that are the speeds associated with the pure steady-state regimes.

From Figure 9, we can conclude:

1. $\bullet$

The average speed is a monotonically increasing function of the load amplitude. 2. $\bullet$

There exists several plateaus, highlighted in red and blue indicating alternating generalised strain regimes and pure steady-state regimes, respectively, where speeds predicted by the analytical model are realised. 3. $\bullet$

Only a finite collection of speeds predicted by the analytical model are realised and lower speeds exist for narrow intervals of the load amplitude or are never realised. 4. $\bullet$

Between any two pure steady-state regimes, corresponding to the fracture process associated with criteria (1), there exists a plateau of finite width at the failure speed associated with the alternating generalised strain failure regime. 5. $\bullet$

The size of plateaus for each steady-state regime corresponding to the alternating generalised strain fracture process increases with the load amplitude. The same behaviour is observed for the plateaus corresponding to the pure steady-state regimes. 6. $\bullet$

Outside the plateaus, the fracture propagates non-uniformly. Between two steady-state speeds, we encounter the clustering phenomenon, first observed in [17], where fracture propagates in regular periodic bursts.

Additionally, with reference to Section 5, we show that the analytical model can be used to efficiently predict the initial load amplitudes when steady-state failure regimes can appear. In Table 1, we present data from the MATLAB simulations when the steady-state regimes are initiated and we supply the classification of each regime. The analytical predictions for these load amplitudes, based on the right-hand side of (40) and the results of [45, 46], are also given. We note that there is an excellent match between the analytical results and those obtained from the numerical scheme implemented in MATLAB.

Note that here we have shown only a finite collection of the predicted speeds is obtained in the simulations. The intervals of the load amplitude for when these steady-state speeds can occur have the form $[{\cal P}_{\min},{\cal P}_{\rm max}]$. For the alternating generalised strain regime, ${\cal P}_{\rm min}$ is determined using (40) and (41), whereas for the pure steady-state regime it is determined using the results of [45, 46]. On the other hand, the maximum value ${\cal P}_{\rm max}$ of the load amplitude for which a steady-state regime can exist may be obtained by analysing the profile of the structure based on the analytical solution, as shown in [46]. Here, for particular failure modes, as the failure speed decreases the distance between ${\cal P}_{\rm min}$ and ${\cal P}_{\rm max}$ becomes smaller (for instance, see Figure 9). In the case when the analytical model predicts ${\cal P}_{\rm min}>{\cal P}_{\rm max}$, which is a physically unacceptable scenario, the corresponding failure regimes cannot be realised.

**Behaviour of the system during the dynamic failure processes **

For ${\cal P}=13.115$, Figure 10(a) shows the system undergoing failure at a particular time $t={1012.4}$ after the fracture of 600 pairs of transverse links. The failure here propagates steadily under the alternating generalised strain regime. The average speed of this process observed from the MATLAB simulations is $\bar{v}=0.6307$. In this case, one can observe the combination of the feeding and reflected waves behind the transition front. Ahead of this point, there exists an evanescent wave. We note that at approximately $t=1034.0$, when the next failure occurs, the profile of the structure can be obtained from that shown in Figure 10(a) mirrored about a horizontal line that corresponds to zero displacement.

In Figure 10(b), we consider a pure steady-state regime resulting from when external load amplitude ${\cal P}=33.484$. The time when the snapshot is taken is $t={403}$. The average steady-state speed achieved is $\bar{v}=1.7490$ and as the failure propagates steadily, the profile of the structure local to the transition front is preserved. Figure 10(b) is an example of the pure steady-state regimes studied in [45, 46] that also occur in the simulations presented here.

In Figures 10(a) and 10(b), we compare the MATLAB simulation data for the arrangement of the masses along the profiles with the analytical results of Section 4 and [45, 46]. In Figure 10(a), for the computations based on the analytical results, the inverse Fourier transform of (28) and (29) is computed using (39), ensuring that (33) holds for the feeding wave phase determined from (32). For the computations in Figure 10(b), the same procedure is followed based on the results of [45, 46]. We emphasise there is again an excellent agreement between the MATLAB results and the analytical predictions.

Failure induced by a high frequency.

Here, we give an example of when the alternating generalised strain regime is never realised. In this illustration, $\omega_{0}=5.9$ and $r=0.4$. As mentioned earlier in Section 6.2, the analytical model predicts a semi-infinite list of possible speeds for both the alternating generalised strain regimes and the pure steady-state regimes (also see the dispersion curves of Figure 5 with the intersection points for $\omega=5.9\,(=\omega_{0})$). However, we observe only one of the predicted steady-state speeds in the MATLAB simulations.

Figure 11(a), shows the average speed $\bar{v}$ as a function of the point load amplitude ${\cal P}$. This particular regime identified is a pure steady-state regime and the average failure speed is $\bar{v}=2.3020$. The load amplitude interval for where this regime is realised is $46.2172\leq{\cal P}\leq 52.612$. There is once again a good agreement with the predictions of analytical model in [45, 46].

In Figure 11(b), we show the position of failure within the supported region as a function of the time when ${\cal P}=46.534$, ${\cal P}=49.7$ and ${\cal P}=51.282$. These values are located within interval defining the plateau for the pure steady-state failure process. For ${\cal P}=46.534$, we see that the transition front speed is constant and equal to the predicted steady-state speed $v=2.3020$, after approximately $44$ breakages. Prior to this, the front propagates non-steadily in the transient regime. For ${\cal P}=49.7$, the pure steady-state regime is reached earlier at approximately the $13^{th}$ breakage. For ${\cal P}=51.282$, we see the failure inside the system settles to the steady-state regime after $100$ breakages.

For the same amplitudes, in Figure 11(c), we present the instantaneous speed $v_{i}$, defined by (44), as a function of the position of fracture. The instantaneous speeds again oscillate about the speed predicted by the analytical model. It is apparent from Figure 11(c) that for those load amplitudes, situated at the extremes of the steady-state plateau in Figure 11(a), the failure process will behave irregularly during a large initial period of the failure process. For ${\cal P}=49.700$, which is situated inside the plateau of Figure 11(a), the process is seen to converge to a more regular oscillatory behaviour earlier in the failure process.

6.3 Case 2 - stiff supports

Here we analyse the failure process in a structure possessing transverse links that are “stiffer” than the links along the central axis of the structure.

Stiffer supports and intermediate load frequency. Here the supports in the system are characterised by the parameter $r=3.4$. We first study the case when this structure is subjected to an oscillating force with a frequency $\omega_{0}=~{}3.1$. Although the flexural stiffness of the transverse links inside the supported region is now larger compared to that of the beams along the central axis, the list of admissible steady-state speeds for failure, outlined in Section 6.2, remains the same. The MATLAB simulations described below indicate which of these regimes are realised.

In Figure 12(a), we show the average speed $\bar{v}$ as a function of $\cal P$. We recall the results of Section 6.2, where the same load frequency was considered and the structure was supported by transverse links with a lower flexural stiffness. There, it was possible to find six failure regimes (both alternating generalised strain regimes and pure steady-state regimes) with different speeds.

When $r=3.4$, we only observe one of the predicted steady-state speeds, $v=1.7484$, which corresponds to a pure steady-state regime with the highest possible speed. The regimes analysed in Sections 2-4 are not observed. In comparing with the computations of Figure 9 for the case of a structure with softer supports, we note that the stiffer supports require larger load amplitudes to initiate and maintain the failure propagation, as expected. Moreover, the stiffer supports allow for a larger range of load amplitudes where the highest steady-state speed can be observed (compare Figure 12(a) with Figure 9). The plateau in Figure 12(a) is defined by the interval $99.452\leq{\cal P}\leq 211.972$. This agrees with the predictions based on the results in [45, 46].

Figure 12(b), shows the position of failure against time for ${\cal P}=99.687,195.62$ and $211.815$, located in the interval defining the plateau for the steady-state regime. For the load amplitude ${\cal P}=195.62$ the time interval for the initial transient failure process appears to be small and the steady-state regime is realised with the average speed $\bar{v}=1.7484$. For ${\cal P}=99.687$, the transition front converges to the steady-state regime after approximately 27 breakages. For the largest amplitude ${\cal P}=211.815$, the steady-state failure process is achieved after approximately 48 breakages. It is interesting to note that for this larger load amplitude, the non-steady failure process attributed to the transient regime is visible and the profile propagates in steps that overlap (see inset of Figure 12(b)). This failure process is known as a forerunning fracture and was first identified in the study concerning the separation of a beam from an elastic foundation [39].

In addition, the normalised instantaneous speeds $v_{i}$, calculated using (44), are shown in Figure 12(c) as a function of the failure position. These computations reveal large fluctuations in the transition front speed about the analytically predicted speed $1.7484$, (represented by unity on the vertical axis in this figure). It is evident that the failure process corresponding to the load amplitudes ${\cal P}=99.687$ and $211.815$ (located at the extremes of the plateau in Figure 12(a)) give large fluctuations in the instantaneous speed of the transition front at the beginning of the failure process. After some time, for all cases shown, the instantaneous speeds oscillate about the predicted speed value in a regular manner and the amplitude of the oscillations decrease as the fracture process develops (see inset of Figure 12(c)).

The example presented here shows that the number of failure regimes observed for different structures subjected to a remote sinusoidal load with a fixed frequency can vary. In this case, the stiffer supports have reduced the number of possible failure regimes that are realisable.

High frequency loading of structure with stiff supports. Finally, we consider a structure with stiff supports and this structure is subjected to a sinusoidal load with a high frequency, which is located near the upper boundary of the passband for unsupported part of the structure.

This particular example shows that the alternating generalised strain regime can be realised in high-frequency loading configurations. In fact, the results of the MATLAB simulations show there are two regimes encountered, as shown in Figure 13(a). In this figure, it can be seen that plateaus exist for $90.386\leq{\cal P}\leq 105.341$ and ${218.57}\leq{\cal P}\leq 376.664$. These intervals again agree well with the analytical predictions. The plateau corresponding to the alternating generalised strain regime (with predicted speed $1.0344$) is very narrow in comparison to plateau belonging to the pure steady-state regime (having the predicted speed $2.3020$). If we compare with the computations of Figure 7 for the same loading frequency, but for the structure with softer supports, we see the present configuration admits additional failure regimes at lower speeds.

Figure 13(b) shows how the position of failure varies as a function of time for load amplitudes chosen inside the plateaus of Figure 13(a). The non-steady behaviour of the transition front in the transient regime is noticeable and there exists a visible variation in the failure speed. After this initial period, the system begins to settle to the steady-state failure process. The instantaneous speed distribution obtained using (44) for the cases analysed in Figure 13(b) is presented in Figure 13(c). It is again evident that there are large oscillations in failure speed about the analytically predicted values in the initial period of the failure process. When the failure process settles, the instantaneous speeds follow a regular oscillatory pattern about the predicted speed. In this example, the computations for the pure steady-state regime (${\cal P}=279.376$) appear converge the fastest to the predicted speed (see inset of Figure 13(c)) than those for the alternating generalised strain regime. Finally, we comment on the instantaneous speed distributions, shown in Figure 13(d), for some load amplitudes within the lowest plateau corresponding to the alternating generalised strain regime. The instantaneous speeds associated with ${\cal P}=92.851$ and ${\cal P}=105.177$, at the extremes of the load amplitude interval defining this plateau, clearly require a large duration of time to settle to the steady-state failure speed. In fact in the case of ${\cal P}=105.177$, several higher speeds are observed up to the 140th instant of failure. Following this, the instantaneous speed eventually settles into steady oscillations about the predicted speed from the analytical model for all cases shown in Figure 13(d).

7 Conclusions

We have developed an analytical model characterising a particular propagation regime in a flexural system, composed of beams connecting periodically placed masses. The regime identified corresponds to the scenario where the generalised strains, represented by bending moments and shear forces in the beam connections, alternate in sign as the failure advances through the system.

Analytical results characterising the dynamics of the system during its failure and when possible failure regimes occur have been determined. In particular, the results have been shown to provide excellent predictions for the behaviour of sufficiently long finite systems subjected to sinusoidal loads that were modelled with a numerical scheme developed in MATLAB. These simulations revealed that the existence of regimes modelled here is dependent on the parameters describing the structure and sinusoidal loading. In addition, the failure regimes identified in [45, 46] can also be realised together with those studied here in the MATLAB simulations.

The model can be generalised to consider different types of loading and more complicated failure criteria, where other types of special failure phenomena may exist. Applications of the tools developed here are envisaged in civil engineering, where the dynamics and failure of large multi-structures are considered.

**Acknowledgements: **M.J.N. gratefully acknowledges the support of the EU H2020 grant MSCA-IF-2016-747334-CAT-FFLAP. I.S.J. would like to thank the EPSRC (UK) for its support through Programme Grant no. EP/L024926/1.

Appendix

A Derivation of the governing equations in terms of generalised coordinates of the masses

The displacements $W_{m}$ and $W^{\pm}_{j}$, along the horizontal and transverse beams, respectively, can be found by solving the equations

[TABLE]

and

[TABLE]

for $t>0$, $\tilde{x}=x-am$, $\tilde{y}_{+}=y$ and $\tilde{y}_{-}=y+a$. The function $W_{m}$ satisfies the boundary conditions

[TABLE]

The function $W_{j}^{+}$ is subject to the conditions

[TABLE]

whereas $W_{j}^{-}$ satisfies

[TABLE]

Thus, from the above we have

[TABLE]

and

[TABLE]

Note, in the considered phase transition problem, the transverse supports fail at the connections to rigid interface. In this case, behind the transition front, the transverse beams are attached to the masses at one end and are free at the other. At the free end, the bending moments and shear forces are zero. Hence, the displacements along these beams can be shown to be linear functions of their local $y$-co-ordinate. Consequently, they do not contribute to the linear and angular momentum balance of the masses behind the transition front (for instance, see (2)).

The bending moments in the transverse beams applied about the $x$-axis can be calculated via

[TABLE]

The angular momentum balance about the $y$-axis for a mass in the supported region then gives

[TABLE]

where the right-hand side is zero as it is assumed all masses have negligible moments of inertia. This together with (LABEL:Wp) and (A.3) implies $\theta^{x}_{m}=0$ in the supported region. The shear forces in these beams then have the form

[TABLE]

and using (LABEL:Wp) and (A.3) leads to (6). The internal bending moments and shear forces in the $m^{th}$ beam aligned with the horizontal axis, (see Figure 2) are computed via

[TABLE]

respectively. Insertion of (A.1) into these relations yields (4) and (5).

B An oscillating point force in an infinite mass-beam chain

Here, we derive the amplitude associated with the outgoing waves generated by point force in a mass-beam chain. The force is assumed to have frequency $\omega_{0}$ and amplitude $P_{0}$. According to the derivation of the governing equations in Section 2.2, the problem may be written as

[TABLE]

and

[TABLE]

Here the Kronecker delta $\delta_{i,j}$ has been used to represent the position of the force in the structure. We consider the complex solution of this problem, i.e. we look for the solutions as

[TABLE]

which allows us to consider the problem in terms of complex functions $\mathfrak{w}_{m}$ and $\vartheta^{y}_{m}$:

[TABLE]

and

[TABLE]

Here we non-dimensionalise the equations by introducing

[TABLE]

and we will also assume the solutions take the form

[TABLE]

We drop the “tilde" in going forward and derive the system

[TABLE]

and

[TABLE]

where all variables in the above equations are dimensionless. In the following, the discrete Fourier transforms are used

[TABLE]

We take the discrete Fourier transform to obtain

[TABLE]

and

[TABLE]

Combining (B.9) and (B.10) we have

[TABLE]

and applying the inverse discrete Fourier transform yields

[TABLE]

For outgoing waves from the source to exist we require that $\omega_{0}<\sqrt{48}$, representing the pass band of the structure. In this case, in the interval $-\pi<k<\pi$ the chosen load frequency can be linked with two wave numbers $\pm k_{f}\mp\text{i}0$, that define the waves propagating in the medium. These wavenumbers are simple poles of the kernel in the above integral. The above is then rewritten as

[TABLE]

where the positive function $R_{0}$ is defined in (35) and has no zeros on the real axis. Applying the residue theorem, under the assumption $m\to\infty$, the integral on the right-hand side has the asymptote

[TABLE]

Finally, combining this with (B.6) and (B.3) leads to

[TABLE]

Thus, from this we see that the amplitude of the feeding wave in the problem considered in Section 5 is

[TABLE]

where the variables in this expression are dimensionless.