Strictly speaking, the solution expressed by Equation (18) holds provided that the interface constant *kz* is not equal to the following "critical" value:

$$\pounds\_z = \frac{4}{B} \frac{\mathbf{d}\_1 + \mathbf{d}\_2}{\left(\mathbf{c}\_1 + \mathbf{c}\_2\right)^2} \,\mathrm{}\tag{21}$$

for which Equation (20) has coincident roots and the solution requires a different expression. However, for material properties corresponding to common fiber-reinforced laminates, the numerical values of *kz* are usually greater than ˆ *kz*, and the roots of the characteristic equation are real-valued (see Appendix B.1). For very compliant interfaces, e.g., if the present model is applied to the analysis of adhesively-bonded joints with weak adhesive, it may happen that *kz* < ˆ *kz*, and the characteristic equation has complex roots. In such cases, the present solution is still valid, provided that the complex exponential functions are considered in Equation (18) and the following. Similar considerations apply also to the unbalanced case.

### 3.2.2. Internal Forces

By substituting Equation (18) into (9), the latter into (8), and integrating with respect to *x*, we obtain the expressions for the internal forces in the upper sublaminate,

$$\begin{aligned} M\_1(\mathbf{x}) &= \frac{f\_5}{\lambda\_5} \exp\left(-\lambda\_5 \mathbf{x}\right) + \frac{f\_6}{\lambda\_6} \exp\left(-\lambda\_6 \mathbf{x}\right) - f\_7 \mathbf{x} - f\_8, \\\ Q\_1(\mathbf{x}) &= \sum\_{i=1}^4 \frac{f\_i}{\lambda\_i} \exp\left(-\lambda\_i \mathbf{x}\right) - f\_{10}, \quad \text{and} \\\ M\_1(\mathbf{x}) &= -\sum\_{i=1}^4 \frac{f\_i}{\lambda\_i^2} \exp\left(-\lambda\_i \mathbf{x}\right) + \frac{h\_1 f\_5}{\lambda\_5} \exp\left(-\lambda\_5 \mathbf{x}\right) + \frac{h\_1 f\_6}{\lambda\_6} \exp\left(-\lambda\_6 \mathbf{x}\right) - \left(f\_{10} + h\_1 f\_7\right) \mathbf{x} - f\_{12}, \end{aligned} \tag{22}$$

and in the lower sublaminate,

$$\begin{aligned} M\_2\left(\mathbf{x}\right) &= -\frac{f\_5}{\lambda\_5} \exp\left(-\lambda\_5 \mathbf{x}\right) - \frac{f\_6}{\lambda\_6} \exp\left(-\lambda\_6 \mathbf{x}\right) + f\_7 \mathbf{x} + f\_9, \\ Q\_2\left(\mathbf{x}\right) &= -\sum\_{i=1}^4 \frac{f\_i}{\lambda\_i} \exp\left(-\lambda\_i \mathbf{x}\right) + f\_{11\prime} \quad \text{and} \\ M\_2\left(\mathbf{x}\right) &= \sum\_{i=1}^4 \frac{f\_i}{\lambda\_i^2} \exp\left(-\lambda\_i \mathbf{x}\right) + \frac{h\_2 f\_5}{\lambda\_5} \exp\left(-\lambda\_5 \mathbf{x}\right) + \frac{h\_2 f\_6}{\lambda\_6} \exp\left(-\lambda\_6 \mathbf{x}\right) + \left(f\_{11} - h\_2 f\_7\right) \mathbf{x} + f\_{13\prime} \end{aligned} \tag{23}$$

where *f*8, *f*9, ..., *f*13 are integration constants.

### 3.2.3. Strain Measures

By substituting Equations (22) and (23) into (4), we obtain the expressions for the strain measures in the upper sublaminate,

1 (*x*) = (a1 + b1*h*1) *f*5 *λ*5 exp (−*λ*5*x*) + *f*6 *λ*6 exp (−*λ*6*x*) − b1 4 ∑ *i*=1 *fi λ*2 *i* exp (−*λix*)+ − [(a1 + b1*h*1) *f*7 + b1 *f*10] *x* − a1 *f*8 − b1 *f*12, *γ*1 (*x*) =c1 4 ∑ *i*=1 *fi λi* exp (−*λix*) − c1 *f*10, and *κ*1 (*x*) = − d1 4 ∑ *i*=1 *fi λ*2 *i* exp (−*λix*) + (d1*h*1 + b1) *f*5 *λ*5 exp (−*λ*5*x*) + *f*6 *λ*6 exp (−*λ*6*x*) + − [(d1*h*1 + b1) *f*7 + d1 *f*10] *x* − b1 *f*8 − d1 *f*12, (24)

and in the lower sublaminate,

2 (*x*) = − (a2 − b2*h*2) *f*5 *λ*5 exp (−*λ*5*x*) + *f*6 *λ*6 exp (−*λ*6*x*) + b2 4 ∑ *i*=1 *fi λ*2 *i* exp (−*λix*)+ + [(a2 − b2*h*2) *f*7 + b2 *f*11] *x* + a2 *f*9 + b2 *f*13, *γ*2 (*x*) = − c2 4 ∑ *i*=1 *fi λi* exp (−*λix*) + c2 *f*11, and *κ*2 (*x*) =d2 4 ∑ *i*=1 *fi λ*2 *i* exp (−*λix*) + (d2*h*2 − b2) *f*5 *λ*5 exp (−*λ*5*x*) + *f*6 *λ*6 exp (−*λ*6*x*) + − [(d2*h*2 − b2) *f*7 − d2 *f*11] *x* + b2 *f*9 + d2 *f*13. (25)

### 3.2.4. Generalized Displacements

Lastly, by substituting Equations (24) and (25) into (1) and integrating the latter with respect to *x*, we obtain the expressions for the generalized displacements in the upper sublaminate,

*Appl. Sci.* **2019**, *9*, 3560

*u*1 (*x*) = − (a1 + b1*h*1) *f*5 *λ*2 5 exp (−*λ*5*x*) + *f*6 *λ*2 6 exp (−*λ*6*x*) + b1 4 ∑ *i*=1 *fi λ*3 *i* exp (−*λix*)+ − [(a1 + b1*h*1) *f*7 + b1 *f*10] *x*2 2 − (a1 *f*8 + b1 *f*12) *x* + *f*14, *w*1 (*x*) = 4 ∑ *i*=1 d1 *λ*2 *i* − c1 *fi λ*2 *i* exp (−*λix*) − (d1*h*1 + b1) *f*5 *λ*3 5 exp (−*λ*5*x*) + *f*6 *λ*3 6 exp (−*λ*6*x*) + + [(d1*h*1 + b1) *f*7 + d1 *f*10] *x*3 6 + (b1 *f*8 + d1 *f*12) *x*2 2 − (*f*16 + c1 *f*10) *x* + *f*15, and *ϕ*1 (*x*) =d1 4 ∑ *i*=1 *fi λ*3 *i* exp (−*λix*) − (d1*h*1 + b1) *f*5 *λ*2 5 exp (−*λ*5*x*) + *f*6 *λ*2 6 exp (−*λ*6*x*) + − [(d1*h*1 + b1) *f*7 + d1 *f*10] *x*2 2 − (b1 *f*8 + d1 *f*12) *x* + *f*16, (26)

and in the lower sublaminate,

*u*2 (*x*) = (a2 − b2*h*2) *f*5 *λ*2 5 exp (−*λ*5*x*) + *f*6 *λ*2 6 exp (−*λ*6*x*) − b2 4 ∑ *i*=1 *fi λ*3 *i* exp (−*λix*)+ + [(a2 − b2*h*2) *f*7 + b2 *f*11] *x*2 2 + (a2 *f*9 + b2 *f*13) *x* + *f*17, *w*2 (*x*) = − 4 ∑ *i*=1 d2 *λ*2 *i* − c2 *fi λ*2 *i* exp (−*λix*) − (d2*h*2 − b2) *f*5 *λ*3 5 exp (−*λ*5*x*) + *f*6 *λ*3 6 exp (−*λ*6*x*) + + [(d2*h*2 − b2) *f*7 − d2 *f*11] *x*3 6 − (b2 *f*9 + d2 *f*13) *x*2 2 − (*f*19 − c2 *f*11) *x* + *f*18, and *ϕ*2 (*x*) = − d2 4 ∑ *i*=1 *fi λ*3 *i* exp (−*λix*) − (d2*h*2 − b2) *f*5 *λ*2 5 exp (−*λ*5*x*) + *f*6 *λ*2 6 exp (−*λ*6*x*) + − [(d2*h*2 − b2) *f*7 − d2 *f*11] *x*2 2 + (b2 *f*9 + d2 *f*13) *x* + *f*19, (27)

where *f*14, *f*15, ..., *f*19 are further integration constants.

### 3.2.5. Supernumerary Integration Constants

In the obtained analytical solution, nineteen integration constants appear. However, seven of them are supernumerary as a consequence of the adopted solution strategy (see Section 3.1). To determine the supernumerary integration constants, we substitute the solution for the interfacial stresses, Equation (18), and generalized displacements, Equations (26) and (27), into Equation (7). In this way, two identities are obtained, which need to be fulfilled for all values of *x*. The following seven relations among the integration constants are determined:

$$\begin{aligned} f\_{10} &= -\frac{\mathfrak{a}\_3}{Bk\_x \mathfrak{a}\_4} \mathbf{d}\_2' f\_{7'} & f\_{11} &= \frac{\mathfrak{a}\_3}{Bk\_x \mathfrak{a}\_4} \mathbf{d}\_1 f\_{7'} \\ f\_{12} &= -\frac{\left(\mathfrak{a}\_1 + \mathfrak{b}\_1 h\_1\right) \mathfrak{d}\_2 - \left(\mathfrak{b}\_2 - \mathfrak{d}\_2 h\_2\right) \mathfrak{b}\_1}{\mathfrak{a}\_4} f\_8 - \frac{\mathfrak{a}\_2 \mathfrak{d}\_2 - \mathfrak{b}\_2^2}{\mathfrak{a}\_4} f\_9, \\ f\_{13} &= \frac{\mathfrak{a}\_1 \mathfrak{d}\_1 - \mathfrak{b}\_1^2}{\mathfrak{a}\_4} f\_8 + \frac{\left(\mathfrak{a}\_2 - \mathfrak{b}\_2 h\_2\right) \mathfrak{d}\_1 - \left(\mathfrak{b}\_1 + \mathfrak{d}\_1 h\_1\right) \mathfrak{b}\_2}{\mathfrak{a}\_4} f\_9, \\ f\_{14} &= f\_{17} - \left(h\_1 + h\_2\right) f\_{19} - \frac{1}{\mathfrak{k}\_x} \left[\frac{1}{B} + \frac{\mathfrak{a}\_3}{B\mathfrak{a}\_4} \left(\mathfrak{c}\_1 \mathfrak{d}\_2 - \mathfrak{c}\_2 \mathfrak{d}\_1\right) h\_1\right] f\_{7'} \\ f\_{15} &= f\_{18}, \quad \text{and} \quad f\_{16} = f\_{19} + \frac{\mathfrak{a}\_3}{B\mathfrak{k}\_x \mathfrak{a}\_4} \left(\mathfrak{c}\_1 \mathfrak{d}\_2 - \mathfrak{c}\_2 \mathfrak{d}\_1\right) f\_{7'} \end{aligned} \tag{28}$$

where:

$$
\mathfrak{a}\_4 = \mathfrak{b}\_1 \mathfrak{d}\_2 - \mathfrak{b}\_2 \mathfrak{d}\_1 + \mathfrak{d}\_1 \mathfrak{d}\_2 (h\_1 + h\_2). \tag{29}
$$

The values of the remaining twelve independent integration constants (*f*1, *f*2, ..., *f*9 and *f*17, *f*18, and *f*19) have to be determined by imposing the boundary conditions of the particular problem being solved. It can be noted that the last three independent integration constants describe a rigid motion of the specimen and can be assessed only when appropriate kinematic restraints are prescribed. The adopted solution strategy enables the determination of the internal forces, strain measures, and relative displacements, even though the values of the last three constants remain undetermined.

### *3.3. Unbalanced Case*

### 3.3.1. Interfacial Stresses

In the unbalanced case, Equations (14) are coupled. To separate the unknowns, we solve the first of such equations w.r.t. to d*τ*/ d*x* and substitute the result into the second equation. After some simplifications, we obtain:

$$\begin{aligned} \frac{\mathbf{d}^6 \sigma}{\mathbf{dx}^6} - \beta\_1 \frac{\mathbf{d}^4 \sigma}{\mathbf{dx}^4} + \beta\_2 \frac{\mathbf{d}^2 \sigma}{\mathbf{dx}^2} - \beta\_3 \sigma &= 0, \\\frac{\mathbf{d} \tau}{\mathbf{dx}} = -\frac{1}{k\_z \beta\_0} \left( \frac{\mathbf{d}^4 \sigma}{\mathbf{dx}^4} - \alpha\_1 \frac{\mathbf{d}^2 \sigma}{\mathbf{dx}^2} + \alpha\_2 \sigma \right), \end{aligned} \tag{30}$$

where:

$$
\beta\_1 = a\_1 + a\_3, \quad \beta\_2 = a\_2 + a\_1 a\_3, \quad \text{and} \quad \beta\_3 = a\_2 a\_3 - k\_x k\_z \beta\_{0\prime}^2 \tag{31}
$$

with *α*1, *α*2, *α*3, and *β*0 already defined in Equations (15) and (16), respectively.

The general solution for the interfacial stresses in the unbalanced case can be written as follows:

$$\begin{aligned} \sigma\left(\mathbf{x}\right) &= \frac{1}{B} \sum\_{i=1}^{6} g\_i \exp\left(-\mu\_i \mathbf{x}\right), \\ \tau\left(\mathbf{x}\right) &= -\frac{k\_\mathbf{x} \beta\_0}{B} \sum\_{i=1}^{6} \frac{g\_i \exp\left(-\mu\_i \mathbf{x}\right)}{\mu\_i \left(\mu\_i^2 - \alpha\_3\right)} + \frac{g\_7}{B} \end{aligned} \tag{32}$$

where *g*1, *g*2, ..., *g*7 are integration constants and *μ*1, *μ*2, ..., *μ*6 are the non-zero roots of the characteristic equation of the differential problem (30):

$$
\mu^6 - \beta\_1 \mu^4 + \beta\_2 \mu^2 - \beta\_3 = 0. \tag{33}
$$

By substituting Equation (31) into (33), the latter can be rearranged as:

$$
\left(\mu^2 - \lambda\_1^2\right)\left(\mu^2 - \lambda\_3^2\right)\left(\mu^2 - \lambda\_5^2\right) + k\_x k\_z \beta\_0^2 = 0,\tag{34}
$$

where *λ*1, *λ*3, and *λ*5 are given by Equation (19). From Equation (34), it can be easily concluded that the roots of the characteristic equation for the unbalanced case are different from those for the balanced case. Conversely, by assuming *β*0 = 0, Equation (34) turns out to be equivalent to Equation (20).

In general, Equation (33) can be solved numerically. In Appendix B.2, the conditions are given for which *μ*1, *μ*2, ..., *μ*6 are real, along with their analytical expressions in this case. It should be noted that, strictly speaking, the solution expressed by Equation (32) is valid only when the characteristic equation has no coincident roots, which is the case for common geometric and material properties.

### 3.3.2. Internal Forces

By substituting Equation (32) into (9), the latter into (8), and integrating with respect to *x*, we obtain the expressions for the internal forces in the upper sublaminate,

$$\begin{aligned} \mathcal{N}\_{1}\left(\mathbf{x}\right) &= -k\_{\mathbf{x}}\beta\_{0}\sum\_{i=1}^{6}\frac{\mathcal{G}i}{\mu\_{i}^{2}\left(\mu\_{i}^{2}-\mathfrak{a}\_{3}\right)}\exp\left(-\mu\_{i}\mathbf{x}\right) - \mathfrak{g}\gamma\mathbf{x} - \mathfrak{g}\mathbf{s}, \\\ \mathcal{Q}\_{1}\left(\mathbf{x}\right) &= \sum\_{i=1}^{6}\frac{\mathcal{G}i}{\mu\_{i}}\exp\left(-\mu\_{i}\mathbf{x}\right) - \mathfrak{g}\_{10\prime} \quad \text{and} \\\ \mathcal{M}\_{1}\left(\mathbf{x}\right) &= -\sum\_{i=1}^{6}\frac{\mathcal{G}i}{\mu\_{i}^{2}}\left(1 + \frac{h\_{1}k\_{\mathrm{x}}\beta\_{0}}{\mu\_{i}^{2}-\mathfrak{a}\_{3}}\right)\exp\left(-\mu\_{i}\mathbf{x}\right) - \left(h\_{1}\mathfrak{g}\gamma + \mathfrak{g}\_{10}\right)\mathbf{x} - \mathfrak{g}\_{12}. \end{aligned} \tag{35}$$

and in the lower sublaminate,

$$\begin{aligned} N\_2\left(\mathbf{x}\right) &= k\_x \pounds\_0 \sum\_{i=1}^6 \frac{g\_i}{\mu\_i^2 \left(\mu\_i^2 - \alpha\_3\right)} \exp\left(-\mu\_i \mathbf{x}\right) + g\_7 \mathbf{x} + g\_9, \\ Q\_2\left(\mathbf{x}\right) &= -\sum\_{i=1}^6 \frac{g\_i}{\mu\_i} \exp\left(-\mu\_i \mathbf{x}\right) + g\_{11}, \quad \text{and} \\ M\_2\left(\mathbf{x}\right) &= \sum\_{i=1}^6 \frac{g\_i}{\mu\_i^2} \left(1 - \frac{h\_2 k\_x \pounds\_0}{\mu\_i^2 - \alpha\_3}\right) \exp\left(-\mu\_i \mathbf{x}\right) - \left(h\_2 \mathbf{g}\_7 - \mathbf{g}\_{11}\right) \mathbf{x} + g\_{13}. \end{aligned} \tag{36}$$

where *g*8, *g*9, ..., *g*13 are integration constants.

### 3.3.3. Strain Measures

By substituting Equations (35) and (36) into (4), we obtain the expressions for the strain measures in the upper sublaminate,

1 (*x*) = − 6 ∑ *i*=1 *gi μ*2 *i kxβ*0 *μ*2 *i* − *α*3 (a1 + b1*h*1) + b1 exp (−*μix*) + − [(a1 + b1*h*1) *g*7 + b1*g*10] *x* − a1*g*8 − b1*g*12, *γ*1 (*x*) =c1 6 ∑ *i*=1 *gi μi* exp (−*μix*) − c1*g*10, and *κ*1 (*x*) = − 6 ∑ *i*=1 *gi μ*2 *i* d1 + *kxβ*0 *μ*2 *i* − *α*3 (d1*h*1 + b1) exp (−*μix*) + − [(d1*h*1 + b1) *g*7 + d1*g*10] *x* − b1*g*8 − d1*g*12, (37)

and in the lower sublaminate,

2 (*x*) = 6 ∑ *i*=1 *gi μ*2 *i kxβ*0 *μ*2 *i* − *α*3 (a2 − b2*h*2) + b2 exp (−*μix*) + + [(a2 − b2*h*2) *g*7 + b2*g*11] *x* + a2*g*9 + b2*g*13, *γ*2 (*x*) = − c2 6 ∑ *i*=1 *gi μi* exp (−*μix*) + c2*g*11, and *κ*2 (*x*) = 6 ∑ *i*=1 *gi μ*2 *i* d2 − *kxβ*0 *μ*2 *i* − *α*3 (d2*h*2 − b2) exp (−*μix*) + − [(d2*h*2 − b2) *g*7 − d2*g*11] *x* + b2*g*9 + d2*g*13. (38)

### 3.3.4. Generalized Displacements

Lastly, by substituting Equations (37) and (38) into (1) and integrating the latter with respect to *x*, we obtain the expressions for the generalized displacements in the upper sublaminate,

*u*1 (*x*) = 6 ∑ *i*=1 *gi μ*3 *i kxβ*0 *μ*2 *i* − *α*3 (a1 + b1*h*1) + b1 exp (−*μix*) + − [(a1 + b1*h*1) *g*7 + b1*g*10] *x*2 2 − (a1*g*8 + b1*g*12) *x* + *g*14, *w*1 (*x*) = 6 ∑ *i*=1 *gi μ*2 *i* d1 *μ*2 *i* − c1 + *kxβ*0 (d1*h*1 + b1) *μ*2 *i μ*2 *i* − *α*3 exp (−*μix*) + + [(d1*h*1 + b1) *g*7 + d1*g*10] *x*3 6 + (b1*g*8 + d1*g*12) *x*2 2 − (c1*g*10 + *g*16) *x* + *g*15, and *ϕ*1 (*x*) = 6 ∑ *i*=1 *gi μ*3 *i kxβ*0 *μ*2 *i* − *α*3 (d1*h*1 + b1) + d1 exp (−*μix*) + − [(d1*h*1 + b1) *g*7 + d1*g*10] *x*2 2 − (b1*g*8 + d1*g*12) *x* + *g*16, (39)

and in the lower sublaminate,

*u*2 (*x*) = − 6 ∑ *i*=1 *gi μ*3 *i kxβ*0 *μ*2 *i* − *α*3 (a2 − b2*h*2) + b2 exp (−*μix*) + + [(a2 − b2*h*2) *g*7 + b2*g*11] *x*2 2 + (a2*g*9 + b2*g*13) *x* + *g*17, *w*2 (*x*) = − 6 ∑ *i*=1 *gi μ*2 *i* d2 *μ*2 *i* − c2 − *kxβ*0 (d2*h*2 − b2) *μ*2 *i μ*2 *i* − *α*3 exp (−*μix*) + + [(d2*h*2 − b2) *g*7 − d2*g*11] *x*3 6 − (b2*g*9 + d2*g*13) *x*2 2 + (c2*g*11 − *g*19) *x* + *g*18, and *ϕ*2 (*x*) = 6 ∑ *i*=1 *gi μ*3 *i kxβ*0 *μ*2 *i* − *α*3 (d2*h*2 − b2) − d2 exp (−*μix*) + − [(d2*h*2 − b2) *g*7 − d2*g*11] *x*2 2 + (b2*g*9 + d2*g*13) *x* + *g*19, (40)

where *g*14, *g*15, ..., *g*19 are further integration constants.

### 3.3.5. Supernumerary Integration Constants

As for the balanced case (see Section 3.2.5), nineteen integration constants appear in the analytical solution for the unbalanced case. Again, seven of them are supernumerary and can be determined by substituting the solution for the interfacial stresses, Equation (32), and generalized displacements, Equations (39) and (40), into Equation (7). In this way, two identities are obtained, which need to be fulfilled for all values of *x*. Consequently, the following seven relations among the integration constants are determined:

*g*10 = − *α*3 *Bkxα*4 d2 + *β*0 b2−d2*h*2 *Bα*4 *g*7, *g*11 = *α*3 *Bkxα*4 d1 − *β*0 b1+d1*h*1 *Bα*4 *g*7, *g*12 = −(a1+b1*h*1)d2−(b2−d2*h*2)b1 *α*4 *g*8 − a2d2−b2 2 *α*4 *g*9, *g*13 = a1d1−b2 1 *α*4 *g*8 + (a2−b2*h*2)d1−(b1+d1*h*1)b2 *α*4 *g*9, *g*14 = *g*17 − (*h*1 + *h*2) *g*19 − 1 *Bkx* + *α*3 *Bkxα*4 (c1d2 − c2d1) *h*1 + *β*0 b1c2−c1b2+c1d2*h*2+c2d1*h*1 *Bα*4 *h*1 *g*7, *g*15 = *g*18, and *g*16 = *g*19 + *α*3 *Bkxα*4 (c1d2 − c2d1) + *β*0 b1c2−c1b2+c1d2*h*2+c2d1*h*1 *Bα*4 *g*7, (41)

where *α*4 is furnished by Equation (29).

The values of the remaining twelve independent integration constants (*g*1, *g*2, ..., *g*9 and *g*17, *g*18, and *g*19) have to be determined by imposing the boundary conditions of the particular problem being solved. Again, the last three independent integration constants describe a rigid motion of the specimen and can be assessed only when appropriate kinematic restraints are prescribed. The adopted solution strategy enables the determination of internal forces, strain measures, and relative displacements, even though the values of the last three constants remain undetermined.

### **4. Delamination Crack Analysis**

### *4.1. Energy Release Rate*

Under I/II mixed-mode fracture conditions, the *energy release rate* can be written as G = GI + GII, where GI and GII are the contributions related to Fracture Modes I and II, respectively. In LEFM, the energy release rate identifies the path-independent J -integral introduced by Rice [57]. Choosing an integration path encircling the interface between sublaminates, as detailed in Appendix C, we obtain:

$$\mathcal{J} = \frac{1}{2}k\_z \Delta w^2(0) + \frac{1}{2}k\_x \Delta u^2(0) \tag{42}$$

or, by recalling Equation (5),

$$\mathcal{J} = \frac{\sigma\_0^2}{2k\_z} + \frac{\tau\_0^2}{2k\_x} \,\mathrm{}^{\prime}\tag{43}$$

where *σ*0 = *σ*(0) and *τ*0 = *τ*(0) are the peak values of the interfacial stresses at the delamination front. The two addends in Equation (43) correspond to the Mode I and Mode II contributions to the energy release rate, respectively [41]. Since compressive normal stresses at the crack tip (*σ*0 < 0) do not promote crack opening, the modal contributions to the energy release rate can be written as

$$\mathcal{G}\_{\mathrm{I}} = \mathcal{H}(\sigma\_0) \frac{\sigma\_0^2}{2k\_z} \quad \text{and} \quad \mathcal{G}\_{\mathrm{II}} = \frac{\tau\_0^2}{2k\_x} , \tag{44}$$

where H(·) denotes the Heaviside step function. By substituting Equation (18) into (44), we obtain the following expressions for the balanced case:

$$\mathcal{G}\_{\rm I} = \frac{\mathcal{H}(\sigma\_0)}{2k\_z B^2} \left(\sum\_{i=1}^4 f\_i\right)^2 \quad \text{and} \quad \mathcal{G}\_{\rm II} = \frac{1}{2k\_x B^2} \left(f\_5 + f\_6\right)^2. \tag{45}$$

Likewise, by substituting Equation (32) into (44), we obtain the following expressions for the unbalanced case:

$$\mathcal{G}\_{\rm I} = \frac{\mathcal{H}(\sigma\_0)}{2k\_z B^2} \left( \sum\_{i=1}^6 g\_i \right)^2 \quad \text{and} \quad \mathcal{G}\_{\rm II} = \frac{1}{2k\_x B^2} \left( k\_x \beta\_0 \sum\_{i=1}^6 \frac{\mathcal{G}i}{\mu\_i \left( \mu\_i^2 - \alpha\_3 \right)} \right)^2. \tag{46}$$

Lastly, to characterize the ratio between the modal contributions to the energy release rate, we introduce the *mode-mixity angle* [11],

$$\psi = \text{sign}\left(\tau\_0\right) \arctan\sqrt{\frac{\mathcal{G}\_{\text{II}}}{\mathcal{G}\_{\text{I}}}}\,\tag{47}$$

ranging from 0◦ (pure Mode I) to 90◦ (pure Mode II).

### **5. Results**

### *5.1. Examples*

By way of illustration, we applied the developed solution to the analysis of the DCB test (Figure 3). Laminated specimens with the following geometric sizes were considered: *L* = 100 mm, *B* = 25 mm, and *H* = 3.2 mm. A delamination length *a* = 25 mm and a reference load *P* = 100 N were used in all of the following calculations. Boundary conditions (12) and (13) were specialized to the present case by taking the crack tip internal forces as *NA* 1 = *NA* 2 = 0, *QA* 1 = −*QA* 2 = *P*, and *MA* 1 = −*MA* 2 = *Pa* and assuming a CTE length = *L* − *a*.

**Figure 3.** Double cantilever beam test.

Three laminated specimens with different stacking sequences were ideally built as follows. We started from a typical carbon-epoxy fiber-reinforced lamina of thickness *tp* = 0.2 mm and elastic moduli shown in Table 1 [61]. Based on this, three stacking sequences were designed corresponding to uni-directional (UD), cross-ply (CP), and multi-directional (MD) sublaminates (Table 2). Lastly, the DCB test specimens were obtained by variously assembling the sublaminates one over the other:

where the double slash (//) represents the position of the delamination crack.

**Table 1.** Lamina elastic moduli.

The sublaminate plate stiffness matrices, calculated as explained in Appendix A, are reported in Table 3. It can be noted that no bending-extension coupling was exhibited by the UD sublaminate (all terms *Bij* = 0), while bending-extension coupling was shown by both the CP (terms *Bxx* = −*Byy* = 0) and MD sublaminates (all terms *Bij* = 0, except *Bxs* = *Bys* = 0).

**Table 3.** Sublaminate plate stiffnesses.

The sublaminate equivalent beam stiffnesses, calculated via Equations (A10) and (A19), are given in Table 4. It can be noted that the UD sublaminate had the highest extensional and bending stiffnesses, A and D, while the MD sublaminate had the highest bending-extension coupling stiffness, B. The shear stiffness, C, showed little variations among the chosen stacking sequences.

**Table 4.** Sublaminate equivalent beam stiffnesses.

A crucial point for the successful implementation of the present model is the appropriate choice of the values of the elastic interface constants. These can be assigned by conducting a compliance calibration with experimental tests or numerical simulations, if available [62]. In a previous study on asymmetric DCB test specimens [22], compliance calibration with finite element simulations gave:

$$k\_x \approx 16.2 \frac{G\_{13}}{H} \quad \text{and} \quad k\_z \approx 6.9 \frac{E\_3}{H'} \tag{48}$$

which, in the present case, yielded *kx* = 21,845 N/mm3 and *kz* = 19,016 N/mm3. Such values were adopted for the present examples.

Figure 4 shows the plots of the interfacial stresses and internal forces as functions of the *x*-coordinate in the neighborhood of the crack tip for the UD//UD specimen. This specimen was symmetric; hence, the unbalance parameter was *β*0 = 0, and the analytical solution for the balanced case applied (see Section 3.2). For the DCB test, the specimen response also corresponded to pure Mode I fracture conditions: actually, Figure 4a shows that the shear stresses were identically null, while the normal stresses attained a peak value at the crack tip cross-section. Figure 4b shows that also the sublaminate axial forces were null. Instead, Figure 4c,d shows the trends of the sublaminate shear forces and bending moments, which started from the assigned values at the crack tip cross-section and then approached zero asymptotically. Shear forces and bending moments had symmetric trends in the upper and lower sublaminates because of the specimen symmetry.

Figures 5 and 6 show the plots of the interfacial stresses and internal forces as functions of the *x*-coordinate in the neighborhood of the crack tip for the UD//CP and UD//MD specimens, respectively. Both specimens were asymmetric with unbalance parameter *β*0 = 0. Hence, the analytical solution for the unbalanced case applied (see Section 3.3). For the DCB test, the specimen response corresponded to I/II mixed-mode fracture conditions: actually, both Figures 5a and 6a show non-zero

normal and shear stresses, which attained peak values at the crack tip cross-section. All of the internal forces started from the assigned values at the crack tip cross-section and then approached zero asymptotically.

**Figure 4.** UD//UD specimen: (**a**) interfacial stresses, (**b**) axial forces, (**c**) shear forces, and (**d**) bending moments in the crack tip neighborhood.

**Figure 5.** *Cont.*

**Figure 5.** UD//CP specimen: (**a**) interfacial stresses, (**b**) axial forces, (**c**) shear forces, and (**d**) bending moments in the crack tip neighborhood.

Table 5 reports the Mode I and II contributions to the energy release rate computed for the three analyzed specimens, along with the total energy release rate and mode-mixity angle. It is apparent that the relative amount of the Mode II contribution increased with the specimen asymmetry.

**Table 5.** Energy release rate and mode mixity.

**Figure 6.** UD//MD specimen: (**a**) interfacial stresses, (**b**) axial forces, (**c**) shear forces, and (**d**) bending moments in the crack tip neighborhood.

### *5.2. Effects of Interface Stiffness*

Figure 7 shows a schematic of the double cantilever beam specimen loaded by uneven bending moments (DCB-UBM) test [63]. Here, this particular test configuration was used to illustrate the effects of interface stiffness, i.e., of the values of the interface elastic constants, *kx* and *kz*, on the solution according to the proposed model.

**Figure 7.** Double cantilever beam specimen with uneven bending moments.

The same laminated specimens considered in the examples of Section 5.1 were again analyzed, but subjected to different loading conditions. The upper sublaminate was loaded by a unit bending couple, *M*1 = 1 N m, while the lower sublaminate was loaded by a variable bending couple, *M*2. Figure 8 shows the dimensionless Mode I and Mode II contributions to the energy release rate, *γ*I = GI/G and *γ*II = GII/G, as functions of the bending moment ratio, *M*2/*M*1.

For the UD//UD specimen (Figure 8a), the modal contributions to G turned out to be independent of the interface stiffness. This result was due to the mid-plane symmetry of the specimen and the independence of the crack tip bending moments on the delamination length. In this case, also the mode mixity depended uniquely on the bending moment ratio, *M*2/*M*1. In particular, it can be noted how pure Mode I fracture conditions (*γ*II = 0) were obtained for *M*2 = −*M*1, while pure Mode II fracture conditions (*γ*I = 0) were obtained for *M*2 = *M*1. Values of *M*2/*M*1 greater than one also corresponded to pure Mode II because of the incompressibility constraint introduced for GI in Equation (44).

(**b**)

(**c**)

**Figure 8.** Dimensionless modal contributions to energy release rate vs. sublaminate bending moment ratio: (**a**) UD//UD specimen, (**b**) UD//CP specimen, and (**c**) UD//MD specimen.

Instead, for both the UD//CP (Figure 8b) and UD//MD specimens (Figure 8c), the modal contributions to G depended on the values of the interface elastic constants, *kx* and *kz*. For each modal contribution, the plots show five curves with different shades of blue: from lighter to darker, and these correspond to *kx* = *kz* = 102, 103, 104, and 105 N/mm3 ; the darkest curves represent the limit case *kx* = *kz* −→ +∞, which corresponds to a rigid interface model [51]. For the specimens considered, also the mode mixity and the conditions for pure fracture modes depended on the interface stiffness.

The range of values adopted for *kx* and *kz* in the previous plots may look quite wide. Indeed, this range corresponds to real values that may be obtained when performing a compliance calibration with experimental results or numerical simulations [62]. The results obtained suggest that the elastic interface model is potentially capable of highly-accurate prediction and interpretation of experimental results, provided that the appropriate values of the interface elastic constants are selected. To this aim, it would be useful to validate the theoretical predictions obtained with experimental results, such as those obtained by Davidson et al. [64,65]. Such a comparison, similar to the one carried out by Harvey et al. to validate different mixed-mode partition theories [66,67], will be the subject of future work.

### **6. Conclusions**

An elastic interface model was introduced for the analysis of delamination in bending-extension coupled laminated beams with through-the-width delamination cracks. The analytical solution obtained applies to a wide class of multi-directional fiber-reinforced composite laminates, including cross-ply laminates, provided that their stacking sequences allow them to be modeled as plane beams. Examples illustrated the application of the model to the analysis of the double cantilever beam test in I/II mixed-mode fracture conditions. Furthermore, the effects of the interface stiffness on the model predictions were discussed with reference to the DCB specimen loaded with uneven bending moments.

By suitably adapting the boundary conditions, the outlined solution can be easily used to analyze different delamination toughness test configurations, as well as delamination behavior in real structural elements and debonding of adhesively-bonded joints.

The proposed model is potentially capable of highly-accurate prediction and interpretation of experimental results, provided that the appropriate interface stiffness is defined. In this respect, it will be useful to compare the model predictions with the results of experimental tests available in the literature. Future developments also include the extension of the model to fully-coupled laminates, for which plate theory must be used, and I/II/III mixed-mode fracture conditions are expected.

**Author Contributions:** Conceptualization, S.B., P.F., L.T., and P.S.V.; methodology, S.B., P.F., L.T., and P.S.V.; software, P.F. and P.S.V.; validation, P.S.V.; formal analysis, P.F. and L.T.; investigation, P.F. and L.T.; resources, S.B. and P.S.V.; data curation, P.F. and P.S.V.; writing, original draft preparation, P.F., L.T. and P.S.V.; writing, review and editing, S.B.; visualization, P.F. and P.S.V.; supervision, S.B.; project administration, P.S.V.; funding acquisition, S.B. and P.S.V.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:

### **Appendix A. Laminated Beam Stiffnesses**

For a homogenous and special orthotropic beam of width *B* and thickness *H*, the extensional, bending-coupling, shear, and bending stiffnesses can be calculated respectively as:

$$\mathbf{A} = E\_{\mathbf{x}} B H, \quad \mathbf{B} = 0, \quad \mathbf{C} = \frac{5}{6} \mathbf{G}\_{\mathbf{z}\mathbf{x}} B H, \quad \text{and} \quad \mathbf{D} = \frac{1}{12} E\_{\mathbf{x}} B H^3,\tag{A1}$$

where *Ex* and *Gzx* respectively are the longitudinal Young's modulus and the transverse shear modulus of the material.

For laminated beams with generic stacking sequences, the above quantities have to be calculated as explained in the following sections.

### *Appendix A.1. Extension and Bending Stiffnesses*

For a laminated plate with a generic stacking sequence, the constitutive laws are:

$$\begin{Bmatrix} N\_{\mathbf{x}} \\ N\_{\mathbf{y}} \\ N\_{\mathbf{xy}} \end{Bmatrix} = \begin{bmatrix} A\_{\mathbf{x}\mathbf{x}} & A\_{\mathbf{x}\mathbf{y}} & A\_{\mathbf{x}\mathbf{s}} \\ A\_{\mathbf{x}\mathbf{y}} & A\_{\mathbf{y}\mathbf{s}} & A\_{\mathbf{y}\mathbf{s}} \\ A\_{\mathbf{x}\mathbf{s}} & A\_{\mathbf{y}\mathbf{s}} & A\_{\mathbf{s}\mathbf{s}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\kappa}\_{\mathbf{x}} \\ \boldsymbol{\kappa}\_{\mathbf{y}} \\ \gamma\_{\mathbf{xy}} \end{Bmatrix} + \begin{bmatrix} B\_{\mathbf{x}\mathbf{x}} & B\_{\mathbf{x}\mathbf{y}} & B\_{\mathbf{x}\mathbf{s}} \\ B\_{\mathbf{x}\mathbf{y}} & B\_{\mathbf{y}\mathbf{s}} & B\_{\mathbf{s}\mathbf{s}} \\ B\_{\mathbf{x}\mathbf{s}} & B\_{\mathbf{y}\mathbf{s}} & B\_{\mathbf{s}\mathbf{s}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\kappa}\_{\mathbf{x}} \\ \kappa\_{\mathbf{y}} \\ \kappa\_{\mathbf{xy}} \end{Bmatrix} \tag{A2}$$

and:

$$\begin{Bmatrix} M\_{\mathbf{x}} \\ M\_{\mathbf{y}} \\ M\_{\mathbf{xy}} \end{Bmatrix} = \begin{bmatrix} B\_{\mathbf{x}\mathbf{x}} & B\_{\mathbf{xy}} & B\_{\mathbf{x}\mathbf{s}} \\ B\_{\mathbf{xy}} & B\_{\mathbf{yy}} & B\_{\mathbf{ys}} \\ B\_{\mathbf{x}\mathbf{s}} & B\_{\mathbf{ys}} & B\_{\mathbf{ss}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\epsilon}\_{\mathbf{x}} \\ \boldsymbol{\epsilon}\_{\mathbf{y}} \\ \boldsymbol{\gamma}\_{\mathbf{xy}} \end{Bmatrix} + \begin{bmatrix} D\_{\mathbf{x}\mathbf{x}} & D\_{\mathbf{xy}} & D\_{\mathbf{xs}} \\ D\_{\mathbf{xy}} & D\_{\mathbf{yy}} & D\_{\mathbf{ys}} \\ D\_{\mathbf{xx}} & D\_{\mathbf{ys}} & D\_{\mathbf{ss}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\kappa}\_{\mathbf{x}} \\ \boldsymbol{\kappa}\_{\mathbf{y}} \\ \boldsymbol{\kappa}\_{\mathbf{xy}} \end{Bmatrix},\tag{A3}$$

where *Nx*, *Ny*, and *Nxy* are the in-plane forces and *Mx*, *My*, and *Mxy* are the bending moments (per unit plate width); correspondingly,  *x*,  *y*, and *γxy* are the in-plane strains and *κx*, *κy*, and *κxy* are the bending curvatures; *Aij*, *Bij*, and *Dij*, with *i*, *j* ∈ {*x*, *y*,*s*}, are the coefficients of the extensional, bending-extension coupling, and bending stiffness matrices, respectively, to be computed following classical laminated plate theory [49]:

$$\begin{Bmatrix} A\_{ij} \\ B\_{ij} \\ D\_{ij} \end{Bmatrix} = \sum\_{k=1}^{n} Q\_{ij}^{(k)} \begin{Bmatrix} z\_k - z\_{k-1} \\ \frac{1}{2} (z\_k^2 - z\_{k-1}^2) \\ \frac{1}{3} (z\_k^3 - z\_{k-1}^3) \end{Bmatrix},\tag{A4}$$

where *Q*(*k*) *ij* , with *i*, *j* ∈ {*x*, *y*,*s*}, are the in-plane elastic moduli of the *k*th lamina, included between the ordinates *z* = *zk*−1 and *z* = *zk*, in the global reference system, and *n* is the total number of laminae. For an orthotropic lamina, whose principal material axes, 1 and 2, are rotated by an angle *θk* with respect to the global axes *x* and *y*,

*Q*(*k*) *xx* = *Q*(*k*) 11 cos4 *θk* + 2(*Q*(*k*) 12 + 2*Q*(*k*) 66 ) cos2 *θk* sin2 *θk* + *Q*(*k*) 22 sin4 *θk*, *Q*(*k*) *xy* = *Q*(*k*) *yx* = (*Q*(*k*) 11 + *Q*(*k*) 22 − 4*Q*(*k*) 66 ) cos2 *θk* sin2 *θk* + *Q*(*k*) 12 (cos4 *θk* + sin4 *θk*), *Q*(*k*) *xs* = *Q*(*k*) *sx* = (*Q*(*k*) 11 − *Q*(*k*) 12 − 2*Q*(*k*) 66 ) cos3 *θk* sin *θk* + (*Q*(*k*) 12 − *Q*(*k*) 22 + 2*Q*(*k*) 66 ) cos *θk* sin3 *θk*, *Q*(*k*) *yy* = *Q*(*k*) 11 sin4 *θk* + 2(*Q*(*k*) 12 + 2*Q*(*k*) 66 ) cos2 *θk* sin2 *θk* + *Q*(*k*) 22 cos4 *θk*, *Q*(*k*) *ys* = *Q*(*k*) *sy* = (*Q*(*k*) 11 − *Q*(*k*) 12 − 2*Q*(*k*) 66 ) cos *θk* sin3 *θk* + (*Q*(*k*) 12 − *Q*(*k*) 22 + 2*Q*(*k*) 66 ) cos3 *θk* sin *θk*, and *Q*(*k*) *ss* = (*Q*(*k*) 11 + *Q*(*k*) 22 − 2*Q*(*k*) 12 − 2*Q*(*k*) 66 ) cos2 *θk* sin2 *θk* + *Q*(*k*) 66 (cos4 *θk* + sin4 *θk*), (A5)

where the in-plane elastic moduli of the lamina in its principal material reference system are:

*Appl. Sci.* **2019**, *9*, 3560

$$\begin{aligned} Q\_{11}^{(k)} &= \frac{E\_1^{(k)}}{1 - \nu\_{12}^{(k)} \nu\_{21}^{(k)}}, \quad Q\_{12}^{(k)} = Q\_{21}^{(k)} = \frac{\nu\_{12}^{(k)} E\_2^{(k)}}{1 - \nu\_{12}^{(k)} \nu\_{21}^{(k)}},\\ Q\_{22}^{(k)} &= \frac{E\_2^{(k)}}{1 - \nu\_{12}^{(k)} \nu\_{21}^{(k)}}, \quad \text{and} \quad Q\_{66}^{(k)} = G\_{12}^{(k)} \end{aligned} \tag{A6}$$

with *E*(*k*) 1 , *E*(*k*) 2 , *ν* (*k*) 12 , *ν* (*k*) 21 , and *G*(*k*) 12 being the in-plane elastic moduli of the lamina in engineering notation.

For a plate to be modeled as a plane beam in the *zx*-plane, the out-of-plane internal forces and moments must be *Ny* = *Nxy* = 0 and *My* = *Mxy* = 0. Besides, the shear-extension and bend-twist coupling terms must vanish: *Axs* = *Ays* = 0 and *Dxs* = *Dys* = 0. The bending-shear coupling terms must also be null: *Bxs* = *Bys* = 0. Such conditions are fulfilled by laminates with single or multiple special orthotropic layers and antisymmetric cross-ply laminates. Furthermore, antisymmetric angle-ply and generic multi-directional laminates may satisfy the above-mentioned conditions for particular stacking sequences. As a consequence, Equations (A2) and (A3) reduce to:

$$
\begin{Bmatrix} N\_x \\ 0 \\ 0 \end{Bmatrix} = \begin{bmatrix} A\_{xx} & A\_{xy} & 0 \\ A\_{xy} & A\_{yy} & 0 \\ 0 & 0 & A\_{ss} \end{bmatrix} \begin{Bmatrix} \epsilon\_x \\ \epsilon\_y \\ \gamma\_{xy} \end{Bmatrix} + \begin{bmatrix} B\_{xx} & B\_{xy} & 0 \\ B\_{xy} & B\_{yy} & 0 \\ 0 & 0 & B\_{ss} \end{bmatrix} \begin{Bmatrix} \kappa\_x \\ \kappa\_y \\ \kappa\_{xy} \end{Bmatrix} \tag{A7}
$$

and:

$$\begin{Bmatrix} M\_{\mathbf{x}} \\ 0 \\ 0 \end{Bmatrix} = \begin{bmatrix} B\_{\mathbf{x}\mathbf{x}} & B\_{\mathbf{x}\mathbf{y}} & 0 \\ B\_{\mathbf{x}\mathbf{y}} & B\_{\mathbf{y}\mathbf{y}} & 0 \\ 0 & 0 & B\_{\mathbf{s}\mathbf{s}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\epsilon}\_{\mathbf{x}} \\ \boldsymbol{\epsilon}\_{\mathbf{y}} \\ \gamma\_{\mathbf{xy}} \end{Bmatrix} + \begin{bmatrix} D\_{\mathbf{x}\mathbf{x}} & D\_{\mathbf{x}\mathbf{y}} & 0 \\ D\_{\mathbf{x}\mathbf{y}} & D\_{\mathbf{y}\mathbf{y}} & 0 \\ 0 & 0 & D\_{\mathbf{s}\mathbf{s}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\kappa}\_{\mathbf{x}} \\ \kappa\_{\mathbf{y}} \\ \kappa\_{\mathbf{xy}} \end{Bmatrix}.\tag{A8}$$

For a laminated beam of width *B*, the axial force and bending moment respectively are *N* = *BNx* and *M* = *BMx*; the axial strain and pseudo-curvature are  =  *x* and *κ* = *κx*. Thus, from Equations (A7) and (A8), the constitutive laws are derived as follows:

$$N = \mathsf{A}\varepsilon + \mathsf{B}\kappa \quad \text{and} \quad M = \mathsf{B}\varepsilon + \mathsf{D}\kappa,\tag{A9}$$

where the laminated beam stiffnesses are evaluated as:

$$\begin{aligned} \mathbf{A} &= B \left( A\_{xx} - \frac{A\_{xy}^2 D\_{yy} - 2A\_{xy} B\_{xy} B\_{yy} + A\_{yy} B\_{xy}^2}{A\_{yy} D\_{yy} - B\_{yy}^2} \right), \\ \mathbf{B} &= B \left[ B\_{xx} + \frac{A\_{xy} (B\_{yy} D\_{xy} - B\_{xy} D\_{yy}) - A\_{yy} B\_{xy} D\_{xy} + B\_{xy}^2 B\_{yy}}{A\_{yy} D\_{yy} - B\_{yy}^2} \right], \quad \text{and} \\ \mathbf{D} &= B \left( D\_{xx} - \frac{A\_{yy} D\_{xy}^2 - 2B\_{xy} B\_{yy} D\_{xy} + B\_{xy}^2 D\_{yy}}{A\_{yy} D\_{yy} - B\_{yy}^2} \right), \end{aligned} \tag{A10}$$

and the remaining plate strain measures turn out to be:

$$\begin{aligned} \epsilon\_{y} &= -\frac{A\_{xy}D\_{yy} - B\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^{2}}\epsilon - \frac{B\_{xy}D\_{yy} - B\_{yy}D\_{xy}}{A\_{yy}D\_{yy} - B\_{yy}^{2}}\kappa\_{\prime} \\ \gamma\_{xy} &= 0, \\ \kappa\_{y} &= -\frac{A\_{yy}B\_{xy} - A\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^{2}}\epsilon - \frac{A\_{yy}D\_{xy} - B\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^{2}}\kappa\_{\prime} \quad \text{and} \\ \kappa\_{xy} &= 0. \end{aligned} \tag{A11}$$

### *Appendix A.2. Shear Stiffness*

The shear stiffness of a laminated plate can be evaluated by imposing a strain energy equivalence, as usually done for homogeneous beams [56] and recently extended to laminated plates by Xie et al. [25]. Here, we detail such a procedure for a laminated beam [68–70]. To this aim, first we recall the expression of the normal stress in the *x*-direction within the *k*th lamina [49]:

$$
\sigma\_x^{(k)}(z) = Q\_{xx}^{(k)}\left(\varepsilon\_x + z\kappa\_x\right) + Q\_{xy}^{(k)}\left(\varepsilon\_y + z\kappa\_y\right) + Q\_{xs}^{(k)}\left(\gamma\_{xy} + z\kappa\_{xy}\right). \tag{A12}
$$

Next, we substitute the expressions of the plate strain measures obtained in the previous section into Equation (A12). Thus, we obtain:

$$
\sigma\_{\mathbf{x}}^{(k)}(\mathbf{z}) = \mathbb{Q}\_{\mathbf{c}0}^{(k)}\boldsymbol{\varepsilon} + \mathbb{Q}\_{\mathbf{x}0}^{(k)}\boldsymbol{\kappa} + z\left(\mathbb{Q}\_{\mathbf{c}1}^{(k)}\boldsymbol{\varepsilon} + \mathbb{Q}\_{\mathbf{x}1}^{(k)}\boldsymbol{\kappa}\right),\tag{A13}
$$

where:

$$\begin{split} \mathbf{Q}\_{\mathbf{c}0}^{(k)} &= \mathbf{Q}\_{xx}^{(k)} - \frac{A\_{xy}D\_{yy} - B\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^2} \mathbf{Q}\_{xy}^{(k)}, \quad \mathbf{Q}\_{\mathbf{x}0}^{(k)} = -\frac{B\_{xy}D\_{yy} - B\_{yy}D\_{xy}}{A\_{yy}D\_{yy} - B\_{yy}^2} \mathbf{Q}\_{xy}^{(k)} \\ \mathbf{Q}\_{\mathbf{c}1}^{(k)} &= -\frac{A\_{yy}B\_{xy} - A\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^2} \mathbf{Q}\_{xy}^{(k)}, \quad \text{and} \quad \mathbf{Q}\_{\mathbf{x}1}^{(k)} = \mathbf{Q}\_{xx}^{(k)} - \frac{A\_{yy}D\_{xy} - B\_{xy}B\_{yy}}{A\_{yy}D\_{yy} - B\_{yy}^2} \mathbf{Q}\_{xy}^{(k)}. \end{split} \tag{A14}$$

Following the classical derivation of Jourawski's shear stress formula [56], we impose the static equilibrium in the *x*-direction of a short portion of beam defined by a given value of the *z*-coordinate. Thus, by recalling the inverse constitutive Equation (4) and static equilibrium Equation (8), with no distributed loads, the following expression is obtained for the shear stress in the *k*th lamina:

$$\pi\_{\varepsilon x}^{(k)}(z) = \mathcal{Q}\left(\mathcal{R}\_k - S\_k z - T\_k z^2\right),\tag{A15}$$

where *Q* is the transverse shear force (acting in the *z*-direction) and:

$$\begin{aligned} R\_k &= S\_k z\_k + T\_k z\_k^2 + \sum\_{j=k+1}^n \left[ \mathcal{S}\_j \left( z\_j - z\_{j-1} \right) + T\_j \left( z\_j^2 - z\_{j-1}^2 \right) \right], \\ S\_k &= \mathbf{b} \mathcal{Q}\_{\mathbf{c}0}^{(k)} + \mathbf{d} \mathcal{Q}\_{\mathbf{c}0}^{(k)} \quad \text{and} \quad T\_k = \frac{1}{2} \left( \mathbf{b} \mathcal{Q}\_{\mathbf{c}1}^{(k)} + \mathbf{d} \mathcal{Q}\_{\mathbf{x}1}^{(k)} \right), \end{aligned} \tag{A16}$$

where b and d respectively are the bending-extension and bending compliances of the laminated beam, as given by Equation (3).

The shear compliance, c, can now be calculated by imposing the following energy equivalence for a short beam segment:

$$\frac{1}{2}\mathbf{c}Q^2 = \frac{1}{2}B\int\_{-h}^{h} \frac{1}{G\_{zx}^{(k)}} \left[\mathbf{r}\_{zx}^{(k)}(z)\right]^2 \,\mathrm{d}z,\tag{A17}$$

where:

$$G\_{\rm zx}^{(k)} = \frac{1}{\frac{\sin^2 \theta\_k}{G\_{23}^{(k)}} + \frac{\cos^2 \theta\_k}{G\_{31}^{(k)}}} \tag{A18}$$

is the transverse shear modulus of the *k*th lamina in the laminate reference system, depending on its out-of-plane shear moduli, *G*(*k*) 23 and *G*(*k*) 31 . By substituting Equation (A15) into (A17), the shear compliance is obtained:

$$\begin{split} \mathbf{c} = \mathcal{B} \sum\_{k=1}^{n} \frac{1}{G\_{\mathbf{z}\mathbf{x}}^{(k)}} \Big[ R\_k^2 \left( z\_k - z\_{k-1} \right) - R\_k S\_k \left( z\_k^2 - z\_{k-1}^2 \right) + \\ &+ \frac{1}{3} \left( S\_k^2 - 2R\_k T\_k \right) \left( z\_k^3 - z\_{k-1}^3 \right) + \\ &+ \frac{1}{2} S\_k T\_k \left( z\_k^4 - z\_{k-1}^4 \right) + \frac{1}{5} T\_k^2 \left( z\_k^5 - z\_{k-1}^5 \right) \Big]. \end{split} \tag{A19}$$

Lastly, the laminated beam shear stiffness is C = 1/c.

#### **Appendix B. Properties of Characteristic Equation Roots**

### *Appendix B.1. Balanced Case*

The roots, *λ*1, *λ*2, ..., *λ*6, of the characteristic Equation (20) for the uncoupled differential problem are furnished by Equation (19). They attain real values if:

$$
\alpha\_1 > 2\sqrt{\alpha\_2} \quad \text{and} \quad \alpha\_3 > 0,\tag{A20}
$$

that is, recalling the definitions of Equation (15),

$$k\_x > 0, \quad k\_z > \hat{k}\_z, \quad \text{and} \quad \mathbf{a}\_1 + \mathbf{d}\_1 h\_1^2 + \mathbf{a}\_2 + \mathbf{d}\_2 h\_2^2 > 2\mathbf{b}\_2 h\_2 - 2\mathbf{b}\_1 h\_1. \tag{A21}$$

with ˆ *kz* introduced by Equation (21). The above conditions are usually fulfilled for geometric and material properties corresponding to common fiber-reinforced laminates.

### *Appendix B.2. Unbalanced Case*

The roots, *μ*1, *μ*2, ..., *μ*6, of the characteristic Equation (33) for the coupled differential problem can be computed numerically. They attain distinct real values if the associated cubic equation has positive discriminant,

$$D = \beta\_1^2 \beta\_2^2 - 4\beta\_2^3 - 4\beta\_1^3 \beta\_3 + 18\beta\_1 \beta\_2 \beta\_3 - 27\beta\_3^2 > 0,\tag{A22}$$

and positive roots,

$$\begin{aligned} \zeta\_1 &= \frac{\beta\_1}{3} + \frac{2}{3}\sqrt{\beta\_1^2 - 3\beta\_2} \cos\frac{\theta}{3} > 0, \\ \zeta\_2 &= \frac{\beta\_1}{3} + \frac{2}{3}\sqrt{\beta\_1^2 - 3\beta\_2} \cos\left(\frac{\theta}{3} + \frac{2\pi}{3}\right) > 0, \quad \text{and} \\ \zeta\_3 &= \frac{\beta\_1}{3} + \frac{2}{3}\sqrt{\beta\_1^2 - 3\beta\_2} \cos\left(\frac{\theta}{3} + \frac{4\pi}{3}\right) > 0, \end{aligned} \tag{A23}$$

where:

$$\theta = \cos^{-1} \frac{2\beta\_1^3 - 9\beta\_1\beta\_2 + 27\beta\_3}{2\sqrt{\left(\beta\_1^2 - 3\beta\_2\right)^3}}.\tag{A24}$$

In this case, the characteristic equation roots are:

$$
\mu\_1 = -\mu\_2 = \sqrt{\zeta\_1}, \quad \mu\_3 = -\mu\_4 = \sqrt{\zeta\_2}, \quad \text{and} \quad \mu\_5 = -\mu\_6 = \sqrt{\zeta\_3}. \tag{A25}
$$

### **Appendix C.** J **-Integral Calculation for the Crack-Tip Element**

The path-independent J -integral was introduced by Rice [57] as:

$$\mathcal{J} = \int\_{\Gamma} \left[ \omega \,\mathrm{d}z - \left( t\_x \frac{\partial u}{\partial x} + t\_z \frac{\partial w}{\partial x} \right) \,\mathrm{d}\Gamma \right], \tag{A26}$$

where Γ is an arbitrary integration path surrounding the crack tip, *ω* is the *strain energy density* (per unit volume), *tx* and *tz* are the *stress vector components* (referred to the outer normal to Γ), and *u* and *w* are the *displacement vector components* along the coordinates *x* and *z*.

For the CTE, we choose an integration path, Γ, encircling the interface between sublaminates as depicted in Figure A1. The path is subdivided into three segments, Γ12, Γ23, and Γ34, along which quantities entering the J -integral have the expressions listed in Table A1.

**Figure A1.** Integration path Γ chosen for the evaluation of the J -integral for the CTE.

**Table A1.** Calculation of J -integral.

Accordingly, the following expression results:

$$\begin{split} \mathcal{J} &= \int\_{0}^{\ell} - \left[ -\tau \frac{\mathbf{d}}{\mathbf{dx}} \left( u\_{1} + h\_{1} \boldsymbol{\varrho}\_{1} \right) - \sigma \frac{\mathbf{d} \mathbf{w}\_{1}}{\mathbf{dx}} \right] \, \mathrm{d} \mathbf{x} + \\ &\quad + \int\_{0}^{t} \left[ \frac{\boldsymbol{\Phi}}{t} - \left( t\_{1} \frac{\mathbf{d} \mathbf{u}}{\mathbf{dx}} + t\_{2} \frac{\mathbf{d} \mathbf{w}}{\mathbf{dx}} \right) \right] \, \mathrm{d} \mathbf{z} + \\ &\quad + \int\_{\ell}^{0} - \left[ \tau \frac{\mathbf{d}}{\mathbf{dx}} \left( u\_{2} - h\_{2} \boldsymbol{\varrho}\_{2} \right) + \sigma \frac{\mathbf{d} \mathbf{w}\_{2}}{\mathbf{dx}} \right] (-\mathbf{d} \mathbf{x}) , \end{split} \tag{A27}$$

where *t* is the interface thickness and Φ = Φ(Δ*w*, Δ*u*) is the *interface potential energy* [27], such that:

$$
\sigma = \frac{\partial \Phi}{\partial \Delta w} \quad \text{and} \quad \tau = \frac{\partial \Phi}{\partial \Delta u} . \tag{A28}
$$

By recalling the definition of relative displacements, Equation (6), and assuming a vanishing interface thickness (*t* −→ 0), we obtain:

$$\mathcal{J} = \Phi\left[\Delta w(\ell), \Delta u\left(\ell\right)\right] - \int\_0^\ell \left[\sigma(\mathbf{x}) \frac{\mathbf{d}\Delta w(\mathbf{x})}{\mathbf{d}\mathbf{x}} + \tau(\mathbf{x}) \frac{\mathbf{d}\Delta u(\mathbf{x})}{\mathbf{d}\mathbf{x}}\right] \,\mathrm{d}\mathbf{x}.\tag{A29}$$

Next, by substituting Equation (A28) into (A29) and recalling that dΦ is an exact differential, we deduce:

$$\begin{split} \mathcal{J} &= \Phi \left[ \Delta w(\ell), \Delta u \left( \ell \right) \right] - \int\_{0}^{\ell} \left[ \frac{\partial \Phi}{\partial \Delta w} \frac{\mathrm{d} \Delta w(\mathrm{x})}{\mathrm{d} \mathrm{x}} + \frac{\partial \Phi}{\partial \Delta u} \frac{\mathrm{d} \Delta u(\mathrm{x})}{\mathrm{d} \mathrm{x}} \right] \mathrm{d} \mathrm{x} = \\ &= \Phi \left[ \Delta w(\ell), \Delta u \left( \ell \right) \right] - \int\_{\left[ \mathrm{Aw}(0), \mathrm{Au}(0) \right]}^{\left[ \left\{ \mathrm{Aw}(\ell), \mathrm{Au}(\ell) \right\} \right]} \frac{\partial \Phi}{\partial \Delta w} \mathrm{d} \Delta w + \frac{\partial \Phi}{\partial \Delta u} \mathrm{d} \Delta u = \\ &= \Phi \left[ \Delta w(\ell), \Delta u \left( \ell \right) \right] - \int\_{\left[ \left\{ \mathrm{Aw}(0), \mathrm{Au}(0) \right\} \right]}^{\left[ \left\{ \mathrm{Aw}(\ell), \mathrm{Au}(\ell) \right\} \right]} \mathrm{d} \Phi. \end{split} \tag{A30}$$

Hence,

$$
\mathcal{J} = \Phi \left[ \Delta w(0), \Delta u(0) \right], \tag{A31}
$$

which holds for a general cohesive interface with potential-based traction-separation laws. For a brittle-elastic interface, the potential energy function specializes to:

$$
\Phi(\Delta w, \Delta u) = \frac{1}{2}k\_z \Delta w^2 + \frac{1}{2}k\_x \Delta u^2. \tag{A32}
$$

By substituting Equation (A32) into (A31), we obtain:

$$\mathcal{J} = \frac{1}{2}k\_z \Delta w^2(0) + \frac{1}{2}k\_x \Delta u^2(0),\tag{A33}$$

which is identical to Equation (42).

For the present elastic interface model, the two addends in Equation (A33) naturally identify the Mode I and II contributions to the energy release rate (see Section 4.1). This result holds for both the balanced and unbalanced cases, for which the analytical solution has been derived in Section 3. For a general cohesive interface model, Equation (A31) yields the value of the J -integral, but its decomposition into Mode I and II contributions is not straightforward. The interested reader can find further hints about this issue in the papers by Wu et al. [37,38] and in the references recalled therein.

### **References**

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Wedge-Splitting Test on Carbon-Containing Refractories at High Temperatures**

### **Martin Stückelschweiger 1,\*, Dietmar Gruber 2, Shengli Jin 2 and Harald Harmuth 2**

Received: 4 July 2019; Accepted: 5 August 2019; Published: 8 August 2019

**Abstract:** The mode I fracture behavior of ordinary refractory materials is usually tested with the wedge-splitting test. At elevated temperatures, the optical displacement measurement is difficult because of the convection in the furnace and possible reactions of refractory components with the ambient atmosphere. The present paper introduces a newly developed testing device, which is able to perform such experiments up to 1500 ◦C. For the testing of carbon-containing refractories a gas purging, for example, with argon, is possible. Laser speckle extensometers are applied for the displacement measurement. A carbon-containing magnesia refractory (MgO–C) was selected for a case study. Based on the results obtained from tests, fracture mechanical parameters such as the specific fracture energy and the nominal notch tensile strength were calculated. An inverse simulation procedure applying the finite element method yields tensile strength, the total specific fracture energy, and the strain-softening behavior. Additionally, the creep behavior was also considered for the evaluation.

**Keywords:** high-temperature wedge splitting test; fracture parameters; reducing condition; carbon-containing refractories; strain-softening; fracture energy

### **1. Introduction**

The wedge-splitting test is performed on notched prismatic specimens that enable stable crack formation for relatively large specimens [1,2]. The reasons are the action of the wedge and the relatively high fracture surface to sample volume ratio. The specimen length is 100 mm, the height is also 100 mm, and the thickness 75 mm. A three-dimensional schematic representation of the sample including the load transmission parts is shown in Figure 1. The wedge-splitting test (WST) at room temperature is performed with a loading rate of 0.5 mm/min. During the test, a video extensometer detects the horizontal displacement at the measuring points. The displacement is measured at the height where the rollers contact the load transmission elements.

**Figure 1.** Schematic representation of the wedge-splitting test.

The sample rests on a linear support (3 × 3 × 90 mm3) to allow for the free vertical displacement of the sample except for the support. The linear support, as well as the wedge and the rollers, are made of corundum to withstand the testing temperatures.

From the results, the specific fracture energy *Gf'* and the nominal notch tensile strength σ*NT*, can be calculated directly from the load-displacement curve (Equations (2) and (3)). The horizontal force *FH*, which is needed for Equations (2) and (3), is calculated from the geometry of the wedge and the measured vertical force *FV* (Equation (1)):

$$F\_H = \frac{F\_V}{2\tan\frac{\theta}{2}}\tag{1}$$

In this equation, α is the wedge angle (10◦).

$$G\_{f'} = \frac{1}{A} \int\_0^\delta F\_H d\delta\_H \tag{2}$$

The specific fracture energy *Gf'* is the area below the load-displacement curve divided by the ligament area A. The measured horizontal displacement is defined with δ*H*, and the horizontal displacement at 15% of the maximum load is δ.

Here *Gf'* (Equation (2)) is the specific fracture energy directly evaluated from the load-displacement curve. This parameter includes only the major part of the specific fracture energy because a premature termination of the test is necessary to avoid the contact of the wedge with the ligament. This means that the WST is stopped before the force reaches zero.

The nominal notch tensile strength σ*NT* is calculated according to Equation (3):

$$
\sigma\_{NT} = \frac{F\_{H,\text{max}}}{bh} \left( 1 + \frac{6y}{h} \right) \tag{3}
$$

Here, *FH,max* is the maximum horizontal load, *b* and *h* are the width and the height of the ligament, and *y* is the vertical distance from the loading position to the center of the ligament.

The simplest possibility for the displacement measurement is to make use of the crosshead travel. Unfortunately, this is not accurate because deformations of the testing machine influence the results as

well as the adjustment of the load transmission parts during load application. Out of these reasons, it is advantageous to measure the displacement directly on the sample.

At elevated temperatures, one possibility is to measure the displacement with an optical dilatometer. Here a halogen lamp illuminates the sample on one side and on the other side the shadow cast by the sample against the light source is recorded by a high-resolution complementary metal-oxide-semiconductor (CMOS) camera [3]. With this method, the displacement is measured on one side of the sample only. Although the performance of the wedge splitting test at room temperature with optical displacement measurement is state of the art [4], at elevated temperatures, the optical measurement is far more difficult due to the radiation and convection in the furnace chamber.

High-temperature testing above 600 ◦C of carbon-containing refractories shall be carried out in non-oxidizing atmospheres. Otherwise, the oxidation of carbon *2C(g)* + *O2(g)* = *2CO(g)* takes place [5], and the microstructure and the properties will be changed significantly [6].

To counteract these problems, a new device was developed to perform the WST on carbon-containing materials up to 1500 ◦C in an argon atmosphere and to measure the displacement on the front and rear side of the sample. The horizontal displacements are measured with two laser speckle extensometers (LSE). In the present paper, the features of the new device and the fracture behavior of a resin-bonded magnesia carbon material at elevated temperatures are discussed.

### **2. Material and Sample Preparation**

The tested resin bonded magnesia carbon brick is composed of fused magnesia, large crystal magnesia clinker, and graphite. A typical application is the refractory lining of a steel ladle in the secondary steel metallurgy. The residual carbon content after cooking under reducing conditions was 10%. The composition showed a CaO/SiO2 (C/S) ratio of 3.8. Bricks were heat-treated at 1000 ◦C for 5 h in a coke bed prior to testing; this enhanced the carbon bond and reduced the emission of volatiles during testing. These emissions caused problems for the displacement measurement by laser speckle extensometers.

### **3. Device**

The wedge splitting test device (universal testing device with displacement measurement) was developed in close cooperation with the company ZwickRoell Testing Systems GmbH (Furstenfeld, Austria). An existing universal testing device with a rigid frame and two rotating spindles were equipped with several components in order to perform the wedge-splitting test at high temperatures. A schematic drawing of the device is illustrated in Figure 2. The load cell is integrated in the lower piston and the maximum load capacity is 300 kN at room temperature. The vertical load rate is set to 0.5 mm/min for the wedge-splitting test.

**Figure 2.** High-temperature testing device (courtesy of Messphysik).

In this device, a gas-tight furnace was installed (Figure 2, Figure 3). The furnace has two windows (front and rear) in order to measure the displacement during the test. The maximum testing temperature is 1500 ◦C. In order to protect the carbon in the samples from oxidation, the furnace chamber can be purged with an inert gas. The purging rate depends on the material and on the testing temperature and can be adjusted between 8 and 70 L/h. The adjustment of the purging rate was done with a valve at the digital flowmeter.

**Figure 3.** Gas-tight furnace (courtesy of Messphysik).

For the displacement measurement during the wedge-splitting test, the device was equipped with two laser-speckle extensometers. One is situated in the front and the other one in the rear of the furnace. The sample was illuminated with green laser light (wavelength: 532 nm), creating a speckled pattern on the sample surface. The movement of the speckled pattern was recorded with two full-screen digital cameras. In each evaluation area, a displacement measurement of the speckled pattern was performed. The increase in the distance between the measurement points is calculated. It is called horizontal displacement in this paper [7].

The setup of the furnace including the laser-speckle extensometers is shown in Figure 4 (all units in mm).

**Figure 4.** Furnace with two laser-speckle extensometers (courtesy of Messphysik).

Each extensometer has two lasers and two cameras. The green lines in Figure 4 show the course of the laser beam. The distance of the measuring points was set to 40 mm. This initial measurement length can be adjusted according to the sample geometry. The distance between the position of the camera and the sample surface is 380 mm.

The advantage of a laser-speckle-extensometer compared to a mechanical extensometer is that the sample and the extensometer fingers are not in contact. Contacting extensometer fingers could be damaged during the cracking of the sample, furthermore, a possible deformation of the extensometers could influence the results. An advantage compared to a standard laser extensometer is that an application of measuring marks to increase the contrast is not necessary because of the characteristic speckled pattern. The resolution of the strain measurement is 0.11 μm according to the specifications from ZwickRoell. The scattering of the measurement signal during a high-temperature test was below 3 μm. This value was determined together with the company ZwickRoell in the course of preliminary investigations. The horizontal displacement was measured in the front and rear side, the average value was used for the evaluation. Figure 5 shows the configuration after a test.

**Figure 5.** Wedge splitting test setup and sample after testing.

The alignment of the rolls and the wedge is essential to achieve symmetrical displacement on the front and rear side of the sample.

### **4. Results and Discussion**

This section is divided into three subsections. It provides a concise and precise description of the experimental results, their interpretations, and the conclusions that can be drawn.

### *4.1. Experiments and Evaluation*

Tests with the newly developed device were carried out successfully between room temperature and 1470 ◦C. During the high-temperature experiments, the furnace chamber was purged with argon and a purging rate up to 70 L/h. Figure 6 shows the load-displacement diagrams from room temperature to 1470 ◦C.

**Figure 6.** Load/displacement diagram of the magnesia carbon material measured at different temperatures.

As the oxidation of carbon in the sample increases with increasing temperature, a purging rate was defined for each individual temperature. The purging rate at 1100 ◦C was set to 10 L/h, at 1370 ◦C 40 L/h and to 70 L/h for the test at 1470 ◦C.

The Table 1 below shows the maximum horizontal force *FH,max*, the specific fracture energy *Gf'* and the nominal notch tensile strength σ*NT* according to the Equations (1)–(3).

**Table 1.** Fracture mechanical parameters determined from tests performed at room temperature (RT), 1100 ◦C, 1370 ◦C, and 1470 ◦C.

For the high-temperature tests the specific fracture energy *Gf'* and the nominal notch tensile strength σ*NT* show an increase with increasing temperature. The test at room temperature shows the smallest nominal notch tensile strength but the specific fracture energy *Gf'* is higher than that at 1100 ◦C and lower than that at 1370 ◦C. The ratio between the specific fracture energy *Gf'* and the notch tensile strength σ*NT* are applied as an indicator for brittleness [8].

The following figure illustrates a sample after the experiment at 1470 ◦C. The oxidation of the carbon in the sample is relatively small, only slight traces of oxidation are visible on the surface. The ligament area of the sample (the area inside the red frame) is without any visible oxidation (Figure 7).

**Figure 7.** Sample after a wedge-splitting test at 1470 ◦C.

The weight loss after the test at 1470 ◦C is below 1% (18 g out of 2250 g of a sample in virgin state) for a purging rate of 70 L/h. The reaction of C(s) with O2(g) decreases with an increase in the argon purging rate [9]. The weight loss after the test at 1370 ◦C is 0.8%. For the test at 1100 ◦C, the weight loss decreases to 0.5%.

### *4.2. Inverse Evaluation*

In order to determine other parameters of interest, such as the tensile strength *ft* and the total specific fracture energy *Gf* and to consider the influence of the creep behavior, an inverse evaluation method was applied on a two-dimensional finite element model of the wedge splitting test developed earlier [10]. The commercial software package ABAQUS (6.13, Dassault Sys Simulia Corp, Providence, RI, USA) was used for modeling. The mode I fracture was considered with an exponential strain-softening behavior applied in the ligament area (Figure 8) (Equation (4)).

$$\sigma = f\_l \left\{ 1 - \left( \frac{1 - e^{-\varphi \frac{x}{x\_{\rm nl}}}}{1 - e^{-\varphi}} \right) \right\} \tag{4}$$

**Figure 8.** Finite element model geometry and mesh.

In this equation *f* t is the tensile strength of the material, the parameter *xult* is the ultimate crack opening and ϕ is a non-dimensional parameter that defines the rate of damage evolution [11]. Cohesive elements were applied to the interfaces of the two half specimens to simulate the mode I fracture [11]. Through an adaptive nonlinear least-squares algorithm, the fracture parameters are inversely calculated by fitting the results of finite element simulations to the results of laboratory wedge splitting tests.

The following Figures 9–11 show the load-displacement diagrams and the determined strain-softening behavior for experiments between 1100 ◦C and 1470 ◦C. The green curves in the left diagrams (a) are the measured ones from the wedge splitting tests. The blue curves (a) are the curves simulated by finite element simulations. The red curves in the diagrams show the results with consideration of the creep behavior and are explained in the following chapter. The (b) figures show the strain-softening behavior determined via inverse evaluation from parameters *f* t, *xult*, and ϕ. Results from (b) correspond to (a).

**Figure 9.** Load-displacement diagram (**a**) and strain-softening behavior (**b**) at 1470 ◦C.

**Figure 10.** Load-displacement diagram (**a**) and strain-softening behavior (**b**) at 1370 ◦C.

**Figure 11.** Load-displacement diagram (**a**) and strain-softening behavior (**b**) at 1100 ◦C.

The inversely calculated curve (a) shows only a slight deviation from the tested one. The minimization procedure is automatically terminated when a minimum of residuals is reached [10]. The tensile strength and *xult* show a decrease with decreasing temperature when the creep is not considered.

Table 2 shows the inversely evaluated parameters and brittleness indicators.

**Table 2.** Fracture mechanical parameters and brittleness indicators at different temperatures calculated from ABAQUS simulations.

With increasing temperature, the specific fracture energy *Gf* increases. The tensile strength follows the same trend. At 1470 ◦C and 1370 ◦C, the maximum horizontal forces are similar but the specific fracture energy is lower at 1370 ◦C because of lower horizontal displacement.

Additionally, the ratio between σ*NT* and *f* t was calculated. The difference between σ*NT* and *f* t is expected. Reasons are the action of the notch on the one hand and on the other hand the dependence on the brittleness number [12]. As well as *Gf'*/σ*NT*, the ratio between the real specific fracture energy *Gf* and the tensile strength *f* t is an indicator of the brittleness. This ratio decreases with decreasing temperature.

### *4.3. Inverse Evaluation with the Consideration of Creep*

In this chapter, the tensile creep behavior was considered in the inverse evaluation. The creep behavior of this material was already investigated by the authors [13]. It is clear that creep has an influence on the service behavior and furthermore on the results of the high-temperature wedge splitting test. The aim of these simulations is to calculate the specific fracture energy, which is caused by crack propagation only. The contribution of creep depends on the loading rate. For low loading rates, the influence of creep is comparatively high and vice versa for high loading rates. The Norton–Bailey creep equation was applied (Equation (5)).

$$
\dot{\varepsilon}\_{cr} = \mathcal{K}(T) \sigma^{\text{\textquotedblleft}g} \varepsilon^{a}\_{cr} \tag{5}
$$

Here, σ is the applied load, . ε*cr* is the creep strain rate, ε*cr* is the creep strain, *K(T)* is a temperature-dependent function and the parameters *n* and *a* are the stress and strain exponents, respectively.

These parameters were implemented in the simulation of the wedge-splitting test and the creep behavior was considered in the whole specimen volume. The diagrams (Figures 9–11) in Section 4.2 shows the results in detail. The left diagrams (a) at each temperature show the experimental curve and the fits, both, with and without consideration of creep. The strain-softening behavior, which is plotted on the right side in Figures 9b, 10b and 11b, shows the change of the energy and the tensile strength of the material at each temperature.

Table 3 shows the change of the parameter with respect to temperature when creep is considered. The inverse calculated parameter *Gf* shows an increase of the specific fracture energy with increasing temperature. The ratio between the specific fracture energy and the tensile strength that is related to the brittleness shows, for both cases (with and without consideration of creep), an increase with increasing temperature. The difference between *Gf* with consideration of creep and *Gf* without consideration of creep (Table 2), Δ*Gf*, shows a maximum temperature of 1370 ◦C and is between 15% and 47% for the high-temperature tests.

**Table 3.** Fracture mechanical parameters at different temperatures calculated considering creep and influence on *Gf*.

### **5. Conclusions**

A new testing device to perform the WST was introduced in this paper. This new device allows fracture mechanical characterization of carbon-containing refractories up to 1500 ◦C under reducing conditions. Tests on resin bonded magnesia carbon material were successfully carried out. The oxidation of the carbon in the sample during the high-temperature test can be avoided by inert gas purging during the whole test. The displacement measurement directly on the sample shows a resolution below 3 μm and works stably over the whole temperature range. Based on the data obtained from laboratory experiments, fracture parameters can be determined and implemented in finite element simulations with an inverse estimation procedure. Besides, simulations including the tensile creep behavior were performed and the fracture parameters under pure mode I failure were determined. The influence of creep on the results obtained from the inverse evaluation for the specific fracture energy is between 15% (1100 ◦C) and 47% (1370 ◦C).

**Author Contributions:** Conceptualization, M.S. and D.G.; Methodology, M.S.; Software, S.J.; Validation, D.G., S.J. and H.H.; Formal Analysis, M.S. and. S.J.; Investigation, M.St.; Resources, M.S. and D.G.; Data Curation, D.G.; Writing—Original Draft Preparation, M.S.; Writing—Review & Editing, D.G.; Visualization, M.S.; Supervision, H.H.; Project Administration, D.G.; Funding Acquisition, H.H.

**Acknowledgments:** This research was funded by K1-MET GmbH, metallurgical competence center. The research program of the competence center K1-MET is supported by COMET (Competence Centre for Excellent Technologies), the Austrian program for competence centers. COMET is funded by the Federal Ministry for Transport, Innovation and Technology, the Federal Ministry for Science, Research and Economy, the province of Upper Austria, Tyrol, and Styria, the Styrian Business Promotion Agency. The authors gratefully thank the company ZwickRoell Testing Systems GmbH (former Messphysik Materials Testing GmbH) for the good cooperation and manufacturing of the newly developed device.

**Conflicts of Interest:** The authors declare no conflict of interest

### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **How Soft Polymers Cope with Cracks and Notches**

### **Andrea Spagnoli \*, Michele Terzano, Roberto Brighenti, Federico Artoni and Andrea Carpinteri**

Department of Engineering and Architecture, University of Parma, Parco Area delle Scienze 181/A, 43124 Parma, Italy<EMAIL_ADDRESS>(M.T<EMAIL_ADDRESS>(R.B<EMAIL_ADDRESS>(F.A<EMAIL_ADDRESS>(A.C.) **\*** Correspondence<EMAIL_ADDRESS>Tel.: +39-0521-905-927

Received: 5 February 2019; Accepted: 9 March 2019; Published: 14 March 2019

**Abstract:** Soft matter denotes a large category of materials showing unique properties, resulting from a low elastic modulus, a very high deformation capability, time-dependent mechanical behavior, and a peculiar mechanics of damage and fracture. The flaw tolerance, commonly understood as the ability of a given material to withstand external loading in the presence of a defect, is certainly one of the most noticeable attributes. This feature results from a complex and highly entangled microstructure, where the mechanical response to external loading is mainly governed by entropic-related effects. In the present paper, the flaw tolerance of soft elastomeric polymers, subjected to large deformation, is investigated experimentally. In particular, we consider the tensile response of thin plates made of different silicone rubbers, containing defects of various severity at different scales. Full-field strain maps are acquired by means of the Digital Image Correlation (DIC) technique. The experimental results are interpreted by accounting for the blunting of the defects due to large deformation in the material. The effect of blunting is interpreted in terms of reduction of the stress concentration factor generated by the defect, and failure is compared to that of traditional crystalline brittle materials.

**Keywords:** soft materials; polymers; strain rate; defect tolerance; digital image correlation; stress concentrators; notch blunting

### **1. Introduction**

A class of materials, which are relevant from the point of view of advanced applications, is represented by the so-called soft or highly deformable materials, such as elastomeric polymers, colloids, liquid-crystals polymers, gels, foams, as well as biological materials—such as soft tissues—that are roughly governed by the same mechanical principles [1]. Typically, these materials have mechanical properties falling within the following range; elastic modulus of 0.1–1.5 MPa, tensile strengths of 1–10 MPa, ultimate tensile strains up to 2000%, and fracture energy of 100–1000 J/m2. Soft materials are endowed with unique features and mechanical properties, explaining the great attention that they have been receiving from the scientific community in the last decades. In natural systems, soft tissues are a fertile source of inspiration for advanced applications, with mechanics and biology going hand-in-hand to formulate the underlying mechanical principles and develop new optimized structural materials [2–6].

The peculiar properties exhibited by soft materials are a direct consequence of their complex and entangled molecular structure, which involves millions or billions of atoms forming linear chains. We owe the fundamentals of the chemical- and physics-based mechanical behavior of this class of materials to the research work of Paul J. Flory and P.J. De Gennes [7–9]. In fully amorphous materials, the structure at the nanoscale level consists of a three-dimensional network of polymer chains, linked together at several discrete points identified as cross-links. The mechanical response of these materials at the meso- or macroscale is heavily affected by the amount of entanglement and the number of existing cross-links per unit volume, rather than by the bonding strength existing between the atoms,

as happens in fully crystalline materials (e.g., metals, ceramics, etc.). Upon the application of load, the deformation induces an alignment and unentanglement of the polymeric chains, with the initial amorphous conformation turned into a semicrystalline-like one. Such a phenomenon usually occurs at quite large deformation, thanks to the highly-oriented arrangement of the chains along the tensile direction [10–12]. This so-called strain-induced crystallization is responsible for the typical stiffening behavior which is noticed in the stress–strain curve of polymers at high strain levels [13,14].

Soft materials normally can deform several times their original length, without being damaged or ruptured. If macroscopic flaws are introduced in the form of cuts, cracks or notches, the stretchability can be affected in different measure, depending on the flaw sensitivity of the material. The concept of a material length scale, separating flaw-sensitive from flaw-insensitive rupture, has been initially proposed for hard materials [15] and later extended to soft polymers [16]. Differently from traditional crystalline solids, rupture of soft polymers occurs at large deformation, when the rearrangement of the polymeric chains leads to flaw reshaping and strengthening around the highest strained region. In particular, an existing initial sharp crack blunts significantly before propagation, with the effect of relieving the strain concentration around the tip [17–19]. The resulting feature is an enhanced defect tolerance, with some tough soft materials showing an insensitivity to flaws up to a few millimeters long, in contrast to the typical nanometer scale for hard brittle solids [20]. For instance, biological tissues such as skin are known for their extreme resistance to flaws, which makes virtually impossible to propagate tearing in a stretched sample. Experimental observations on rabbit skin showed that a small notch does not propagate but progressively blunts, due to straightening and stretching of the collagen fibers [21]. However, other soft materials, such as the so-called hydrogels, are often prone to premature fracture and have low fatigue resistance. In such materials, the rupture process is strongly influenced by the fluid interaction [22], and it has been discovered that increased toughness and fatigue resistance are obtained through the development of polymer networks containing chains of different lengths (such as the double-network hydrogels, containing both a short- and a long-chain highly stretchable network). In this fashion, the shortest chains act as sacrificial elements, while the longer ones provide the material with a further elastic behavior [23,24]. From this perspective, it appears that the polymer network characteristics play a crucial role in defining the macroscopic behavior of the material. Furthermore, soft materials are sensible to time-dependent effects, with the fracture energy depending on the rate of application of the external loads, because of the viscous energy dissipation occurring in the crack tip region [25–27]. Sometimes, under a constant applied load, a so-called delayed fracture has been observed, depending on the network structure of the soft matter [28].

In this work we present a comprehensive investigation into the defect tolerance of flawed specimens of rubber-like polymers, with a detailed summary of experimental findings recently published by the authors [29–31]. Various configurations with cracks and notches are examined in order to evaluate the macroscopic mechanical response in relation to the flaw severity. Moreover, the effect of the applied strain rate is taken into account phenomenologically. All the experimental tests are conducted under strain control, and the kinematically-related quantities are measured through a contactless Digital Image Correlation (DIC) technique. Due to the local high deformations in proximity of the defects, a severe defect remodeling with evident blunting is noticed. The main purpose is to show how crack blunting affects the rupture behavior of soft materials. Through a simplified analytical model, the increase in the curvature radius at the notch root with the remote applied load is described. Such a model is applied to the experimental results, putting into evidence that it is the blunting effect which controls the rupture process; in particular, the grade of blunting appears to determine the transition from the typical small-scale yielding failure to rupture at a constant theoretical strength.

The paper is organized as follows. Section 2 presents a collection of the results obtained from the experiments, with details on the materials and methods employed. Accurate images taken from the DIC elaboration are included. In Section 3, a detailed discussion is developed in order to give a comprehensive interpretation of the experiments, with specific reference to the notch blunting effect, which is computed through a simplified analytical model. Finally, Section 4 sets out the conclusions.

### **2. Experimental Tests on Macroscopically Flawed Thin Plates**

In this section, we present the results of tensile tests on thin plates made of various elastomeric polymers, carried out under displacement control up to complete failure. The plates contain flaws with different geometrical features, ranging from cracks to notches, denoted here as geometric discontinuities with a finite radius of curvature. Rubber-like polymers have been chosen because their behavior is highly representative of a vast range of soft materials. For instance, silicone rubber is often employed as a substitute of human skin, since it does not show the strain-induced stiffening of natural rubbers, maintaining comparable values of tensile and tearing strength [32]. The stress–strain curves obtained from tensile testing of silicone rubbers are typical of a hyperelastic, almost incompressible behavior. Several models have been proposed [33], which usually show good agreement at low to moderate stretches (generally below 1.5). During our experiments we have found that the one-parameter neo-Hookean model offers a satisfactory approximation in the stretch range of interest (see, for instance, Figure 1).

**Figure 1.** Typical curve obtained from tensile testing of a silicone rubber (experimental and fitted data). The enlarged inlet shows the region λ = 1.2–1.4, with excellent agreement between the different models.

The purpose of the experiments herein described is to evaluate the effect of initial flaws on the tensile strength of the material, considering both the effect of the flaw size and of the applied strain rate. The response of the specimens during the experimental tests is monitored by measuring the applied force and corresponding displacements, and by using the Digital Image Correlation (DIC) technique. The DIC is a contactless technique widely used to get full-field displacement and strain maps in experiments, through numerical reconstruction of the kinematic field shown by the surface points of the samples. For an optimal use of such a technique, the surface of the specimens needs to be covered with speckle patterns before testing. Several parameters affect measurement accuracy and spatial resolution, included optical measures connected to the camera and lens resolution, image magnification, mean size of the speckle pattern, and factors depending on the correlation algorithm, such as the image subset size and the gray-level interpolation. In particular, the resolution of the displacement measurement is governed by the subset size, whose lower bound is limited by the speckle size and, consequently, by the available pixels [34,35]. In the present work, the images are acquired with a high-resolution digital camera (maximum resolution of 5184 × 3456 pixels) mounted on a tripod, and lights are used to ensure a uniform illumination of the samples. The sequence of

images is treated by means of the freeware software Ncorr [36], developed in MATLAB environment, for monitoring the displacement and the strain fields within the specimen.

### *2.1. Plates Containing Elliptical Flaws*

The first set of experiments is carried out on elastomeric sheets under tension, containing elliptical flaws of length 2*a*, characterized by a finite notch root radius (Figure 2a). The plates are made of Sylgard®, a common silicone polymer having the following elastic parameters: initial Young modulus *E* = 0.84 MPa and Poisson ratio *ν =* 0.37. Three types of specimens have been prepared, containing a centered elliptical flaw with different values of the root radius *ρ* and same length 2*a*. The plate aspect ratio is kept constant at *L/W* ≈ 2. The geometric characteristics of the specimens are reported in Table 1. The flawed plates are subjected to tensile loading along the *y*-axis, applied at a constant strain rate of . *ε*0 = . *δ*/(2*L*) = 4.8·10−3 *s*−1. The tests have been interrupted before failure, since an evident notch remodeling was noticed. The response during the experimental tests is monitored by means of the DIC technique.

**Figure 2.** (**a**) Sketch of a thin plate containing a centered internal flaw. The enlarged view shows the case of an elliptical notch, with semi-axes *a* and *b*. *ρ* is the tip radius of the notch (*ρ* = 0 in the case of a crack-like flaw). (**b**) Sketch of an edge-cracked plate. The plates are subjected to a remote displacement *δ* applied along the *y*-axis.

**Table 1.** Geometric characteristics of the thin plates with elliptical notches. The rightmost column contains the values of the initial notch root radius *ρ* = *b*2/*a*.

Figure 3 shows the initial (undeformed) and the generic stretched shapes of the specimen El1 at two increasing levels of the remote applied stretch. The corresponding strain patterns obtained from the DIC analysis are illustrated in Figure 4, specifically, the Green–Lagrange strain *Eyy* parallel to the loading direction and the strain *Exx* transversal to such a direction. Three different levels of the applied remote stretch are considered, defined as

$$
\lambda\_0 = 1 + \frac{\delta}{2L} \tag{1}
$$

where *δ* is the applied displacement.

**Figure 3.** Qualitative view of the (**a**) undeformed and (**b**,**c**) increasingly deformed configurations, in the flawed specimen El1.

The DIC plots show that the material is able to comply with very high deformations, leading to a severe defect remodeling characterized by an evident blunting of the notch. A compressed region just in front and behind the elliptical hole is observed, due to the contraction effect arising in the direction normal to the applied load (see also Figure 3c). This phenomenon has a beneficial effect in terms of the strain concentration, with a sort of augmented notch blunting due to the local flexural instability of the thin plates in the compressed zones [29].

**Figure 4.** Strain maps obtained from the Digital Image Correlation (DIC) analyses on the elliptically notched specimen El1, for three increasing values of the applied remote stretch λ0. (**a**–**c**) Green–Lagrange strain *Eyy*; (**d**–**f**) corresponding maps for the strains *Exx*. In (a), the length scale is expressed in cm.

### *2.2. Plates Containing Internal Crack-Like Flaws*

A second set of experiments deals with plates containing a centered crack of length 2*a* (Figure 2a with *b* = *ρ* = 0). The samples are made of a commercial silicone rubber (TSE3478T by Momentive). The desired crack is obtained by manually cutting the samples with a sharp blade. From tensile tests on sound specimens, the initial Young modulus of the material is found approximately equal to *E* = 1.12 MPa and the Poisson ratio to *ν* = 0.42. Four series of specimens containing a centered crack are

here considered, characterized by different values of the relative crack length *a/W*, whereas the plate aspect ratio is kept constant at *L/W* ≈ 2. The geometric characteristics of the specimens are reported in Table 2. The plates are subjected to tensile loading along the *y*-axis up to failure, applied at a constant strain rate of . *ε*0 = 5.8·10−3 s−1. The response during the experimental tests is monitored by means of the DIC technique.

**Table 2.** Geometric characteristics of the specimens containing a centered crack. The rightmost column contains the ultimate remote stretch before failure.

Figure 5 shows the deformed patterns of the specimens, for three different stages of the applied loading: at the beginning of the test (undeformed, left column), at an intermediate stage (remote stretch *λ*0 = 1.29, central column), and at the final stage before failure (ultimate remote stretch, right column). The strain patterns obtained from the DIC analyses are also shown, specifically, the Green–Lagrange strain *Eyy* (parallel to the loading direction) at the intermediate and the ultimate stages. It can be noticed that the initial crack-like shape tends to blunt under loading and the applied ultimate stretch generally decreases for an increasing relative crack size *a/W*. At the intermediate stage, the DIC maps clearly show a strain concentration typical of an elliptical notch, where the maximum strain values occur in the locations corresponding to the original crack tips. At incipient failure, the strain maps exhibit a complex distribution due to the failure mechanisms developing in the vicinity of the notch roots. It is worth noticing that out-of-plane displacements occur in two limited regions close to the crack edges (the blue regions in Figure 5), because of the contraction effect arising in the direction normal to the applied load, i.e., in the *x*-direction.

**Figure 5.** Images of the precracked specimens at the initial stage (left column), and with maps of the Green–Lagrange strain *Eyy* at an intermediate stage (central column) and final stage at incipient failure (right column). The remote stretch is shown on the top of each plot. (**a**–**c**) Specimens type CC1 (*a/W* = 0.179); (**d**–**f**) specimens type CC4 (*a/W* = 0.446). In (b,e), the length scale is expressed in cm.

The distribution of the Green–Lagrange strain *Eyy* along the *x*-axis is shown in Figure 6.

**Figure 6.** Green–Lagrange strain *Eyy* distribution along the cracked section of the plate, at an intermediate (continuous lines) and ultimate values (dashed lines) of the remote applied stretch. (**a**) Specimens type CC1 (*a/W* = 0.179) and (**b**) specimens type CC4 (*a/W* = 0.446).

### *2.3. Plates Containing Edge Crack-Like Flaws*

The last set of experiments that we review is related to plates containing edge crack-like defects (Figure 2b), where we consider the effects of the relative crack size *a*/*W* and of the applied strain rate. A first group of edge-cracked plates (specimens EC1–EC4 in Table 3) has been prepared with the commercial silicone RTV 615 (Elantas Europe Srl). From tensile tests on sound specimens, the initial Young modulus of the material is equal to *E* = 1.50 MPa and the Poisson ratio is equal to *ν =* 0.42. The plates are subjected to tensile loading along the *y*-axis up to failure, applied with three different strain rates: . *ε*0 (1) = 1.9·10−3 s−1, . *ε*0 (2) = 0.48·10−3 s−1, and . *ε*0 (3) = 0.16·10−3 s−1. The response during the experimental tests is monitored by means of the DIC technique.

**Table 3.** Geometric characteristics of the edge-cracked plates. The rightmost column contains the ultimate remote stretch before failure. In specimens EC1–EC4, the fastest rate . *ε*0 (1) = 1.9·10−3 s−1 is considered.

A second group of edge-cracked plates (EC5-EC8 in Table 3) has been prepared using a different silicone rubber (Elite Double 32 by Zhermack Dental), with an estimated Young modulus equal to *E* = 1.36 MPa and the Poisson ratio is equal to *ν =* 0.42. The plates are subjected to tensile loading along the *y*-axis up to failure; such a loading is applied with a strain rate . *ε*0 = 3.8 − 4.2·10−3 s−1. Four series of specimens have been prepared for each group, with different values of the relative crack length *a*/*W*, while the plate aspect ratio is kept constant at *L/W* ≈ 1.5.

The results obtained from the DIC analyses in the specimens EC2–EC4 are illustrated in Figure 7, specifically, the Green–Lagrange strain *Eyy* (parallel to the loading direction) at incipient failure for the highest strain rate among the three considered. From a qualitative observation, it can be noticed that crack tip blunting is limited (compared, for instance, to the case of an internal crack in Figure 5).

**Figure 7.** Map of the Green–Lagrange strain *Eyy* at incipient failure, at the applied rate . *ε*0 (1). The remote stretch is shown on the top of each plot. (**a**) Specimens type EC2 (*a/W* = 0.077), (**b**) type EC3 (*a/W* = 0.192) and (**c**) type EC4 (*a/W* = 0.308). In (a), the length scale is expressed in cm.

The results show that specimens with only one ligament zone are more sensible to defects with respect to the ones containing a central crack (CC1–CC4). The eccentric load effect, which increases with the deformation, plays a crucial role in intensifying the stress close to the crack tip, leading to a premature failure of the polymer network chains within the crack process region.

### *2.4. Effects of Intrinsic Material Defects*

The group of experiments on edge-cracked plates has also involved samples with intrinsic defects, in the form of microvoids (D1–D4 in Table 4). For this purpose, we have prepared samples with the same material and sizes of the group EC1–EC4, but following a different treatment during the preparation of the silicone mixture. To prepare the samples, 50 g of component A (matrix) are thoroughly mixed with 5 g of component B (curing agent). At this point, the mixture for the material without defects is carefully degassed in vacuo and subsequently mechanically spread into the custom-made aluminum mold. This stage is followed by a second degassing, then the mixture is cured in oven at 60 ◦C overnight and finally mechanically removed from the mold. In order to obtain microbubbles embedded in the material, the silicone mixture is directly cured in oven without the degassing stages.

**Table 4.** Geometric characteristics of the edge-cracked plates with intrinsic defects. The rightmost column contains the ultimate remote stretch before failure.

A comparison of the materials, with and without the microvoids, is shown in Figure 8a,b. The average void radius is equal to *rv* = 0.35 mm, and its distribution is described by adopting a normal probability function (variance equal to 0.021 mm2). The plot of the probability function (Figure 8c) shows a good agreement with the measured distribution, and ~80% of the voids have a radius within half of the smallest edge crack length considered.

**Figure 8.** Images of the elastomeric material employed in specimens EC1–EC4 and D1–D4. Detail of the material (**a**) without and (**b**) with intrinsic defects. (**c**) Distribution of the void size cumulated probability.

Figure 9 illustrates the failure behavior of edge-cracked specimens EC1–EC4 and D1–D4, in terms of the ultimate stretch vs the relative crack size. As was expected, smaller ultimate stretches are attained for increasing values of the relative crack size *a/W*. The decreasing trend, with a roughly quadratic pattern, confirms that the material is sensitive to the presence of the initial flaws.

**Figure 9.** Ultimate stretch vs the initial relative crack size for the edge-cracked samples. (**a**) Effect of the applied strain rate in the specimens EC1–EC4 and (**b**) influence of the material defects, for a strain rate . *ε*0 = 1.9·10−3 s−1 (results from samples EC1–EC4 are shown as solid circles, and D1–D4 as hollow circles).

The effect of the strain rate is investigated in Figure 9a. The highest strain rate provides the highest ultimate stretch, although the effect tends to reduce at lower ultimate stretches, suggesting that the influence of flaw size on the ultimate stretch decreases if the deformation is applied sufficiently slow. Figure 9b compares the failure in specimens with and without intrinsic material defects. It can be inferred that microdefects tend to anticipate failure, although the results are not so clear, probably due to the uneven distribution of the microbubbles in the samples.

### **3. Discussion**

### *3.1. Model for Notch Blunting*

Cracks and notches in soft materials become distorted under loading due to large strain effects, with an evident tip blunting, as is clearly shown in the figures obtained from the DIC on the experimental tests. In order to quantify the flaw severity and account for the geometrical effect of blunting, we resort to the concept of the stress concentration factor. Considering a notched sample under remote uniform stress acting parallel to the *y*-axis, the stress concentration factor can be defined as

$$K\_l = \sigma\_{\text{max}} / \sigma\_0 \tag{2}$$

where *σ*max = *σyy*,max is the maximum notch root stress and *σ*0 is the corresponding remote stress, measured with respect to a uniformly stressed (gross) section of the sample. In the following, we propose a simplified analytical model, capable of relating the stress concentration factor with the tip radius *ρ* of the blunted notch and explore the variation of *Kt*(*ρ)* under increasing applied loading. Such an approach, conceived for the analysis of blunting of elliptical notches [29], is extended to the tip blunting of cracks (Section 3.2).

The starting point is the solution for an elliptical notch in an infinite elastic plate, having semi-axes *a* and *b*. Its equation can be written as follows

$$y(\mathbf{x}) = \frac{b}{a}\sqrt{a^2 - \mathbf{x}^2} \tag{3}$$

and the radius of curvature at the tips (*x* = ±*a*) is equal to

$$
\rho = b^2 / a
\tag{4}
$$

The stress concentration factor *Kt* for an elliptical notch is obtained from the renowned Inglis' solution [37]:

$$K\_t(\rho) = \left(1 + 2\sqrt{\frac{a}{\rho}}\right) \mathbb{Y}\left(\frac{a}{W'}, \frac{a}{b}\right) \tag{5}$$

where *Y* is a corrective factor introduced to account for the finite width of the specimens [38] (considering the initial geometry of the tested notched specimens, the values of *Y* range from 1.09 to 1.10). For the sake of simplicity, this factor is taken equal to unity in the following equations.

In order to describe the change of the stress concentration factor with the increasing deformation of the blunted tip under loading, we need to consider the variation of the radius of curvature. Such a variation depends on the deformation state in proximity of the notch root, which has to be evaluated in the deformed configuration (see Figure 10a). Let us consider a small square element of edge length equal to *h*, located in the proximity of the notch root, with edges inclined by an angle α with respect to the horizontal axis (Figure 10b). The small strain tensor *ε* in the local reference system 1–2 (with components *ε*11,*ε*22,*ε*12) is related to the corresponding tensor *ε* (with components *εtt*,*εnn*,*εtn*) in the reference system *t-n*, with its origin located at the notch root, through the well-known relationships [33]:

$$
\varepsilon\_{11} = \varepsilon^2 \varepsilon\_{tt} + s^2 \varepsilon\_{nn}, \\
\varepsilon\_{22} = s^2 \varepsilon\_{tt} + \varepsilon^2 \varepsilon\_{nn}, \\
\varepsilon\_{12} = \csc \varepsilon\_{tt} - \csc \varepsilon\_{nn} \tag{6}
$$

with *c* = cos *α*,*s* = sin *α*. Assuming plane stress conditions and the material governed by the generalized Hooke model, the strain tensor components are

$$
\varepsilon\_{tt} = -\nu \sigma\_{nn} / E , \varepsilon\_{nn} = \sigma\_{nn} / E , \varepsilon\_{tn} = 0 \tag{7}
$$

If the angle α is sufficiently small, the only non-zero component of the stress tensor is

$$
\sigma\_{nn} = \mathbb{K}\_l(\rho)\sigma\_0 \tag{8}
$$

The radius of curvature at the notch root, in the undeformed state (point A in Figure 10a,b), is approximated by the radius of the local osculating circle:

$$
\rho \simeq \frac{h}{2\cos\beta} \tag{9}
$$

with *β* = *π*/2 − *α*, and the increased radius of curvature in the deformed state (point A' in Figure 10a,b) is approximated as

$$\rho' = \frac{h(1+\varepsilon\_{22})}{2\cos(\beta+\gamma)}\tag{10}$$

where *γ* = 2*ε*12.

**Figure 10.** (**a**) Model of notch blunting. Schematic of the notch profile modification due to the large strains occurring in the material during an increment of the remote stress *dσ*0. The radius of curvature at point *A* (undeformed) and *A'* (deformed) is shown. (**b**) Deformation of a small material element in the vicinity of the notch root. (**c**) Elliptical notch, with the equivalent semi-axes obtained from the crack flank displacement under tensile stress *σ*0.

The model here described is nonlinear, because the local strains depend on *Kt* through Equations (6)–(8), where, in turn, *ρ* is also a function of the local strains through Equation (10). Following an incremental procedure, firstly the increments of strain are evaluated for an applied remote stress variation *dσ*0 through Equations (6)–(8); then, the updated radius of curvature is obtained from Equation (10) and the stress concentration factor is computed from Equation (5). In other words, following a stepwise updated Lagrangian approach, at the first increment the increased radius of curvature *ρ* of Equation (10) is calculated with respect to the reference configuration on the undeformed notch, whilst at successive increments the reference configuration, i.e., *ρ* from Equation (9), corresponds to the radius of the notch root under the current stress level. The results applied to one of the elliptically-notched plates tested experimentally are illustrated in Figure 11.

**Figure 11.** Stress concentration factor as a function of the remote applied stretch. Prediction of the 'pure' LEFM model (continuous line) and results obtained from the 'modified' model, in the center-cracked specimens CC1 and CC4. The dash-dot line refers to the elliptically-notched plate El2. The enlarged view shows the region 1.1 < *λ* < 1.3.

### *3.2. Crack Tip Blunting*

Linear Elastic Fracture Mechanics (LEFM) is grounded on the hypothesis of infinitesimal displacements, with the near-tip region being described by the K-dominated fields, which imply, for instance, a square root singularity in the strain and stress components at the crack tip. In soft materials, the large strain region near the tip can be relevant, and the local fields need to be defined within the framework of finite strain elastostatic, showing for instance stress singularities different from the well-known square-root [39]. However, tip blunting can still be described according to the simplified model introduced in Section 3.1, provided that the local strain is not too large.

According to LEFM, the deformed shape of a generic crack of semilength *a*, subjected to a uniform stress *σ*0, is described by the following expression of the crack flank displacement [40].

$$y(\mathbf{x}) = 2\frac{\sigma\_0}{E}\sqrt{a^2 - \mathbf{x}^2} \tag{11}$$

which is the equation of an ellipse, having the major semi-axis equal to *a* and the minor semi-axis equal to *b*0 = 2*a σ*0 *E* (Figure 10c). Retrieving the expression of the tip radius of an ellipse from Equation (4), the equivalent radius of curvature is

$$
\rho\_0 = 4a \left(\frac{\sigma\_0}{E}\right)^2\tag{12}
$$

A direct relationship between the applied remote stretch and the stress concentration factor predicted by the Inglis solution, Equation (5), is derived as

$$K\_l(\lambda\_0) = 1 + 2\sqrt{\frac{a}{\rho\_0}} = \frac{\lambda\_0}{\lambda\_0 - 1} \tag{13}$$

where we have used the standard relation for linear elastic materials: *λ*0 = 1 + *σ*0 *E* . The noticeable feature of such an expression is that it does not depend on the initial length of the crack. According to LEFM, the ultimate condition at failure occurs when the remote stress equals

$$
\sigma\_{\mathbb{C}} = \frac{1}{Z\left(\frac{a}{\mathbb{W}}\right)} \sqrt{\frac{G\_{\mathbb{C}}E}{\pi a}}\tag{14}
$$

where *Z a W* is a corrective factor for the finite width of the cracked specimens [41] (considering the initial geometry of the tested notched specimens, the values of *Z* range from 1.02 to 1.69) and *GC* is the fracture resistance of the material. The critical stretch, omitting the corrective factor *Z*, is then obtained as

$$
\lambda\_{\mathbb{C}} = 1 + \sqrt{\frac{G\_{\mathbb{C}}}{\pi E}} a^{-\frac{1}{2}} \tag{15}
$$

Up to this point, we have not considered the deformation of the blunted crack under loading. The predicted stress concentration factor, as is described from this 'pure' LEFM model in Equation (13), is illustrated by the continuous line in Figure 11.

The effects of the blunted tip deformation under increasing loading and the resulting stress concentration factor are described adopting the 'modified' model introduced in Section 3.1, where the current tip radius is computed on the deformed configuration from Equation (10). We assume that the initial configuration of the crack is the critical condition of LEFM, Equation (14), when *σ*0 = *σC* and the corresponding value of the radius of curvature, from Equation (12), is

$$
\rho\_C = \frac{4}{\pi} \frac{G\_C}{E} \tag{16}
$$

The results of the 'modified' model are also shown in Figure 11. We have considered the center-cracked specimens CC1 and CC4, which present two different lengths of the crack: we notice that, contrary to the LEFM model results, there is a dependence of *Kt* on *a*.

### *3.3. Application to Experimental Data*

The model is applied to account for the effect of blunting in the failure of the specimens tested during the experiments. The data collected in Figure 12 illustrate the ultimate condition, when the remote stretch is equal to *λ*0 = *λU*, with respect to the initial length *a* of the defect. In this plot we have also added the predicted trend of LEFM, as obtained from Equation (15), where an average value of the fracture resistance for silicone rubbers is taken equal to *Gc* = 1 kJ/m2 .

**Figure 12.** Ultimate remote stretches at failure *λU* as a function of the initial crack length *a*. In the case of the elliptical notches, the upward arrows indicate that failure has not occurred within the considered range of applied stretch. The continuous line corresponds to the prediction of the 'pure' LEFM solution, Equation (15). The dashed lines correspond to best-fit curves with a *λU* vs. *a* power law dependence with exponent −0.5. Results for EC1–EC4 for a strain rate . *ε*0 = 1.9·10−3 s−1.

It should be noticed that the plot in Figure 12 is by no means a representation of the flaw sensitivity in the canonical form, primarily because the smallest size of the flaw in the tested specimens is larger than the critical flaw size (see [20]). The discriminant here seems to be the grade of crack blunting: the edge-cracked specimens, which fail at smaller stretches, are well approximated by the 'pure' LEFM model, where the ultimate stretch decreases with the square root of *a*. On the contrary, the plates containing central cracks (specimens CC1–CC4) resist to larger stretches and do not agree well with the power law trends, suggesting that crack blunting plays a fundamental role.

To the latter group of specimens, we have applied the 'modified' model previously described. The results are summarized in Figure 13, where the quantity *Kt*(*λU* − 1) is plotted as a function of the crack length *a*. Note that the quantity 1 *E* (*Kt*(*λU* − 1)*E*) represents the normalized true stress at the notch root at incipient failure, and can be therefore considered a material property, independent of the presence of the flaw. From the observation of the results reported in Figure 13, it can be noticed that such a local stress seems to be an appropriate parameter for quantifying the material failure in the case of the centered-cracked plates CC1–CC4, but not in the others. In other words, where crack blunting allows rupture at larger stretches, failure can be predicted by a local stress quantity rather than by a fracture mechanics-related one.

**Figure 13.** Normalized notch root true stress at incipient failure. The effect of crack blunting is enhanced for the specimens CC1–CC4, with the failure stress being approximately constant (dashed line).

### **4. Conclusions**

The present paper investigates the defect tolerance capability of flawed silicone specimens subjected to tensile loading, adopting both a theoretical and an experimental approach. Silicone rubber has been chosen as illustrative of the typical response of soft matter in general, including other polymers, gels, and biological tissues. Different flaw geometries, sizes, and rates of the applied strain are examined in the experimental tests, with the aim of understanding their influence on the macroscopic mechanical response of the samples. The presence of internal defects, in the form of microbubbles, is also considered. The experimental response is monitored by measuring the applied force and using the Digital Image Correlation (DIC) technique to obtain precise two-dimensional displacement and full-field strain maps.

Flawed samples of soft materials subjected to tensile loading undergo a remarkable notch blunting prior to failure, which tends to reduce the stress concentration due to the presence of the flaw. Moreover, when notches are contained in nonconfined thin elements, a further notch blunting occurs as consequence of the local buckling of the material in the compressed zones, which arises normally to the loading direction. Experimental tests on thin silicone plates with elliptical notches showed that, irrespective of the initial size of the notch, very high remote stresses are supported thanks to the favorable notch profile evolution under load. With respect to the strain rate effects, we can observe that the flaw sensitivity is augmented at higher rates while, if a sufficiently slow strain rate is applied, the ultimate strain before failure is less affected by the size of the initial flaw. Indeed, slow deformation rates allow the internal microstructure of the material in high strained regions to be rearranged and to fail locally, with the consequence to smooth out the peak strain arising close to the geometric discontinuities.

A simple analytical model is proposed to account for the effect of crack blunting, in terms of the increase of the tip radius with the remote applied stretch. Although based on the assumption of linear elastic behavior, the application of the model to the experimental results has provided an interesting insight into the defect tolerance of cracked thin plates under tensile loading. Specifically, we have observed a transition from a typical flaw-size-dependent failure, as predicted by linear elastic fracture mechanics, to rupture occurring at a constant theoretical strength of the material. Our results show that, even at flaw lengths larger than the critical size found in other studies [20], silicone rubber specimens can withstand large stretches thanks to the flaw reshaping allowed by their peculiar microstructure. We are confident that the obtained results might be applied profitably to the evaluation of the safety levels of notched soft structural components, commonly found in numerous advanced applications.

**Author Contributions:** Conceptualization, A.S.; Data curation, R.B. and F.A.; Methodology, R.B.; Supervision, A.S.; Visualization, M.T.; Writing–original draft, M.T.; Writing–review & editing, A.C.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Analytical Study of Reinforced Concrete Beams Tested under Quasi-Static and Impact Loadings**

### **Sayed Mohamad Soleimani 1,\* and Sajjad Sayyar Roudsari 2**

Received: 28 May 2019; Accepted: 13 July 2019; Published: 16 July 2019

**Abstract:** During dynamic events (such as impact forces), structures fail to absorb the incoming energy and catastrophic collapse may occur. Impact and quasi-static tests were carried out on reinforced concrete beams with and without externally bounded sprayed and fabric glass fiber-reinforced polymers.  The drop height and loading rate were varied. The load-carrying capacity of reinforced concrete beams under impact loading was obtained using instrumented anvil supports (by summing the support reactions). In quasi-static loading conditions, the beams were tested in three-point loading using a Baldwin Universal Testing Machine. ABAQUS FEA software was used to model some of the tested reinforced concrete beams. It was shown that the stiffness of reinforced concrete beams decreases with increasing drop height. It was also shown that applying sprayed glass fiber-reinforced polymers (with and without mechanical stiffeners) and fabric glass fiber-reinforced polymers on the surface of reinforced concrete beams increased the stiffness. Results obtained from the software analyses were in good agreement with the laboratory test results.

**Keywords:** reinforced concrete beam; impact and quasi-static loading; retrofitting

### **1. Introduction**

Concrete structures are subjected to a variety of different loading conditions in day-to-day situations, with the two most ubiquitous of these conditions being quasi-static and impact loadings. With regard to the former, there exists an abundance of research studying the behavior of reinforced concrete (RC) elements under quasi-static load. Conversely, structures such as high-rise buildings are often exposed to damaging impact forces, including the incidental impact of objects, explosions, sudden collapse of cranes, and unregulated motion of heavy machinery. As such, it is also imperative to study the behavior of RC beams under these dynamic loading conditions. To this end, there have been numerous studies conducted on impact loading and its effects on RC structures, with the first known dynamic test being conducted on concrete in 1917 [1]. After years of inactivity, there has been a resurgence in various impact loading experiments performed on concrete over the past 50 years. Researchers such as Atchley and Furr [2], Scott et al. [3], Dilger et al. [4], Mlakar et al. [5], and Soroushian et al. [6] studied the behavior of RC elements under dynamic conditions and discovered that with incremental increases in loading rate, the ultimate stress and strain of concrete increases by about 25%. Moreover, with respect to the compressive strength of RC elements, experiments carried out by Wastein [7], and Malvar and Ross [8] suggest that the compressive strength of concrete under dynamic loads increased by 85–100%. However, the existing literature on the bending capacity of RC beams in response to impact loading remains comparatively scarce. Bertero et al. conducted experiments on simply supported RC beams [9] and noted that the flexural capacity and rigidity of RC beams show an incremental increase at high rates of strain. Similarly, Wakabayashi et al. [10]

ran dynamic experiments on beams under strain rates of 0.001 s−1 and discovered that the loading capacity of beams increased by 30% in high strain rates. Comparable experiments have also been carried out by Fujikake et al. [11] to experimentally and analytically evaluate the effect of hammer height and longitudinal reinforcement on the behavior of RC beams, using a mass-spring damper system to simulate the impact response. Furthermore, Pham and Hao [12] performed an accurate numerical simulation (less than 12% average absolute error) using an artificial neural network to predict the response of RC beams to impact loading, with concrete compressive strength, hammer velocity, and mass of hammer being considered as variables. Another important variable—namely, contact stiffness—was studied experimentally and numerically by Pham et al. [13], with results indicating that this variable can significantly influence the peak of load–displacement diagram. Moreover, contact stiffness significantly influenced the bending moment and shear force of the RC beams.

Given the advancing age of RC structures, especially those in developed countries, the importance of using advanced materials—such as fiber-reinforced polymers (FRP)—for rehabilitation purposes encompassing both quasi-static and impact loadings should be emphasized. The use of externally bonded FRP as a strengthening technique for RC structures has gained popularity over the past three decades, with a notable surge in its use being observed in recent years. As such, it is necessary to study the behavior of FRP-strengthened RC elements under the aforementioned loading conditions. With respect to quasi-static loading, there has been an increasing trend in the number of studies conducted in recent years [14,15], establishing the role of externally bonded FRP in augmenting the load-carrying capacity of RC elements. On the other hand, researchers have found that the behavior of RC elements strengthened by externally bonded FRP under impact loading conditions is slightly different compared to that under quasi-static loading, with the significance of the correlation between effective FRP-concrete bonding and performance under load being highlighted [16–21].

Considering the existing body of literature, there exists the possibility of further development pertaining to the analytical study of RC beams—specifically, those tested under impact loading, as well as those strengthened by means of sprayed glass fiber-reinforced polymer (GFRP) under quasi-static loading. This paper will attempt to develop upon these topics in order to illuminate the behavior of RC beams treated with sprayed GFRP under quasi-static loading conditions. Greater insight into the use of sprayed GFRP in strengthening RC elements affords several noteworthy advantages over externally bonded fabric FRP, including ease of implementation, cost reduction, and time efficiency. As such, the novelty of this paper lies in two primary domains—the finite element modeling of RC beams strengthened for shear by sprayed GFRP, which remains an understudied area of research, in addition to the presence of mechanical fasteners (i.e., through bolts) in associated modeling schemes.

### **2. Materials and Methods**

A total of 17 RC beams were tested in this study and the experimental results were evaluated using ABAQUS software. Nine RC beams were tested under quasi-static loading including two control and seven retrofitted beams. The other eight beams were used to study the behavior of RC beams under impact loading. Different retrofitting schemes using sprayed GFRP (with and without mechanical fasteners) were employed for shear strengthening of the RC beams tested under quasi-static loading. Fabric GFRP was employed for flexural strengthening. Samples previously experimented in Soleimani's doctoral dissertation [22] were employed for the modeling and analysis of the behavior of RC beams under quasi-static and impact loads. Tables 1–3 show the characteristics of RC beams as well as the properties of sprayed GFRP containing randomly distributed chopped fibers and unidirectional fabric GFRP.

### *2.1. Quasi-Static Loading*