# Identification of a cohesive zone model from digital images at the micron-scale

Julien Réthoré, Rafael Estevez

<figure><img src="image_1.png" /></figure>

## To cite this version:

Julien Réthoré, Rafael Estevez. Identification of a cohesive zone model from digital images at the micron-scale. Journal of the Mechanics and Physics of Solids, 2013, 61 (6), pp.1407-1420. 10.1016/j.jmps.2013.01.011. hal-00850226

HAL Id: hal-00850226

https://hal.science/hal-00850226

Submitted on 10 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Identification of a cohesive zone model from digital images
at the micron-scale

Julien Réthoré<sup>a,*</sup>, Rafael Estevez<sup>b</sup>

<sup>a</sup> Université de Lyon/INSA Lyon/CNRS/LaMCoS, 20 Avenue des Sciences, F-69621 Villeurbanne Cedex, France

<sup>b</sup> SIMaP UMR 5266, Université de Grenoble/CNRS/Grenoble-INP/UJF, 1130 rue de la Piscine, B.P. 75, F-38402 St Martin d'Hères cedex, France

We present a new methodology for the identification of a zone cohesive model that describes material failure. The material under consideration fails by crazing. The study is conducted at the micron scale in order to capture and analyze the fracture mechanism. The crack tip displacement fields are measured optically by Digital Image Correlation. The local stress intensity factors (mode I and II) and the location of the equivalent elastic crack tip are calculated during the loading. The variation of the location of the equivalent crack tip is used to track the initiation and growth of the process zone, up to the onset of crack propagation. These experimental measurements are used to define the appropriate parameters in a cohesive zone model. The methodology addresses the onset of crazing, the traction-separation profile and the maximum opening corresponding to the local nucleation of a crack. The cohesive parameters that are derived from the experimental data are consistent with results available in the literature. In addition, the model enables the characterization of the normal and tangential mode of the cohesive model.

# 1. Introduction

The description of fracture mechanisms by cohesive zone models (Barenblatt, 1962) has been a subject of increased interest in the micromechanics community since the pioneering work of Needleman (1987) and Xu and Needleman (1994). This approach enables a local description of the fracture mechanics that can be physically motivated. The fracture mechanisms is represented by a traction-opening relationship that mimics the failure process. Cohesive zone models are typically characterized by a maximum traction $T_{max}$ at the onset of decohesion and a critical opening $\delta^{cr}$ that corresponds to the nucleation of a crack locally. In metals and ceramics, the opening $\delta^{cr}$ is on the order of nanometers (Gall et al., 2000; Abraham, 2001; Yamakov et al., 2006; Kubair et al., 2009) resulting in an extension of the cohesive process zone by tens of $\delta^{cr}$ (Rice, 1980). These dimensions are smaller than one micrometer and are not accessible by direct optical observation techniques. The use of electronic microscopes could be an alternative, but the sample size required to observe a failure process make them difficult to handle.

In polymer fracture, crazing is the mechanism responsible for failure, a process that is well documented and characterized for amorphous, glassy polymers (Kramer, 1983; Kramer and Berger, 1990). The typical dimensions involved in crazing allow an analysis by optical interferometry and extensive work by Döll and coworkers (Döll, 1983; Döll and Könczöl, 1990) have established that the critical opening is 2–10 µm, resulting in a length for the cohesive process zone of some tens of micrometers. Such dimensions allow for an optical analysis of the crack tip displacement fields: by interferometry as reported by Döll, but also by means of recently developed Digital Image Correlation (DIC) techniques.

* Corresponding author. Tel.: +33 472 438787.

E-mail addresses: julien.rethore@insa-lyon.fr (J. Réthoré), Rafael.Estevez@simap.grenoble-inp.fr (R. Estevez).
Experimental techniques based on DIC for the analysis of the displacement fields are currently expanding due to the improvements in the optical devices and monitoring setups. The methods initiated with the extraction of continuous regular fields (e.g., Sutton et al., 1986; Sun et al., 2005; Besnard et al., 2006) but recent developments using discontinuous enrichment for the displacement decomposition, addressed problems with a singularities like shear bands (Réthoré et al., 2007b) and the analysis of a crack tip field (Réthoré et al., 2007a; Grégoire et al., 2009; Nguyen et al., 2011; Poissant and Barthelat, 2010). Since the earliest development of the DIC technique, fracture mechanics has been a topic of investigation (e.g., McNeill et al., 1987; Abanto-Bueno and Lambros, 2002; Réthoré et al., 2005; Roux and Hild, 2006). However, the aim has usually been to extract stress intensity factors, whereas research for extracting adhesive or cohesive properties of an interface or a failure in a material has only been addressed more recently (Abanto-Bueno and Lambros, 2005; Fedele et al., 2009; Shen and Paulino, 2011; Fuchs and Major, 2011). Theoretical developments had already been proposed in Hong and Kim (2003) for evaluating a cohesive law from elastic far-fields. In the papers mentioned above, the observational scale is far too large compared to the resolution of DIC in order to capture the local state of the material during the crazing phase. Herein, the optical setup is adapted to the scale of the cohesive zone and sub-micron resolution images of the sample surface are acquired. From these images, a direct identification strategy of a cohesive zone law for Poly(Methyl Methacrylate) (PMMA) is established.

Fracture in glassy polymers is governed by the competition between shear yielding and crazing. Shear yielding corresponds to the development of localized plasticity due to softening upon yielding. The magnitude of the shear bands is limited by a progressive re-hardening that takes place during continued deformation. Crazing is also a mechanism of localized plasticity, but at a smaller scale. It involves three characteristic stages: (i) initiation at a local critical stress state, (ii) craze thickening with the formation of craze fibrils up to (iii) craze fibril breakdown corresponding to the local nucleation of a crack. In the literature, these three stages are addressed separately with studies devoted to craze initiation (Sternstein et al., 1968; Sternstein and Myers, 1973; Oxborough and Bowden, 1973; Argon and Hannoosh, 1977), the formation of craze fibrils (Kramer, 1983; Kramer and Berger, 1990) and the craze opening profile measured by optical interferometry (Döll, 1983; Döll and Könczöl, 1990). In particular, Döll and coworkers reported the maximum craze opening for various glassy polymers. The prediction of the cohesive process zone obtained with the Dugdale model (1960) is in good agreement with those observations. However, a reduced domain in terms of loading conditions is explored in Döll and coworkers' experiments and they are restricted to steady state crack growth under a constant load.

Kramer (1983) and Kramer and Berger (1990) pointed out that crazing, and in particular, the fibrillation stage involves some plasticity. Based on Kramer's work, a rate-dependent cohesive zone model was proposed by Tijssens et al. (2000) and Estevez et al. (2000) that accounts for some plasticity and the three characteristic stages of crazing. Rate-dependence of the failure process has been experimentally demonstrated for PMMA by Saad-Gouider et al. (2006). In this study, the loading rates corresponded to a time to rupture that varied from one second and to several hours and a rate of the stress intensity factor ranging between $1 \times 10^{-5}$ MPa$\sqrt{m/s}$ and $1$ MPa$\sqrt{m/s}$. Saad-Gouider et al. (2006) have shown experimentally that under mode I, the energy release rate in PMMA is rate-dependent, and its bulk plasticity is negligible. This observation confirms the need for a rate-dependent cohesive zone model as formulated by Tijssens et al. (2000), that results in a rate-dependent Dugdale model (1960) for monotonic loadings. The calibration presented by Saad-Gouider et al. (2006) focused on the identification of the parameters involved in the description for craze initiation and craze thickening. The value of the critical opening, $\delta^{cr}$, was taken from Döll (1983) and Döll and Könczöl (1990). Although tractable, the extraction of cohesive zone model parameters requires (i) specific experiments to characterize craze initiation, (ii) numerous fracture tests at different loading rates for craze widening and related fibrillation, and (iii) optical interferometry of wedge cracks to measure the craze critical opening, $\delta^{cr}$.

In the present study, an identification of a cohesive zone model for crazing from the analysis of the displacement fields measured by DIC in the vicinity of the crack tip is presented. The material under consideration is a PMMA identical to that investigated by Saad-Gouider et al. (2006). We demonstrate that the proposed analysis provides an estimation of (i) the craze stress at the onset of crazing, (ii) the traction-opening profile during craze widening and (iii) the critical opening corresponding to craze fibril breakdown and onset of crack propagation. These parameters and profile are derived from a single test. The experiment is also simple and provides complementary information to that reported in Kausch (1983, 1990) which focused on specific features of crazing.

The paper is organized as follows. First, the experimental setup and the fracture test from which the crack tip fields are analyzed by DIC is described. In a second section, the background for the analysis of the crack tip displacement fields is presented. Next, stress intensity factor calculations, the monitoring of the crack tip position and related process zone during loading are described. These terms are used to derive an identification protocol for the traction-displacement relationship of the cohesive zone model that is based on the DIC data. The results from the proposed identification protocol are then discussed from a materials science perspective and the novelty of the present strategy is pointed out with respect to available literature.

The following notation is adopted. Tensors are denoted by bold-face symbols. Implicit summation is considered over repeated Latin indices. The summation convention is *not* used for repeated Greek indices.

## 2. Experimental setup

A commercial PMMA (Perspex) in the form of 10 mm thick plates, stored for more than three months under standard room conditions in order to avoid any ageing effects has been used in this study. This is the same material used
by Saad-Gouider et al. (2006). The molecular mass measured by gel permeation chromatography results in $M_n$ = 864 kg/mol and $M_w$ = 1843 kg/mol, in agreement with typical values for PMMA. According to Döll (1983) and Döll and Könczöl (1990), crazing and fibrillation is observed in PMMA as long as $M_n$ > 200 kg/mol; this requirement was satisfied for all samples in the present study. The primary and secondary relaxation temperatures at 1 Hz were 400 K and 284 K respectively.

The preparation of a sharp crack is identical to that reported in Saad-Gouider et al. (2006). Sharp notches are introduced following a two step procedure. A pre-notch of 0.25 mm radius is first created with a circular cutting wheel. In order to prevent heating while machining, this operation is done under a flow of fresh compressed air. A sharp notch is further introduced at the tip of the pre-notch by tapping a razor blade. A mass (200 g) is dropped from a height of 18 cm through a guide fixed on the device. This runner impacts the razor blade placed in the throat of the pre-notch. The procedure results in reproducible sharp cracks, with a length at least four times larger than the radius of the pre-notch, in agreement with the recommendation of the ESIS-TC4 (Williams and Pavan, 2003). The procedure results in negligible stress induced by the notching operation.

The geometry of the sample used in the present analysis is depicted schematically in Fig. 1. A screw is used to prescribe a dominant opening mode. The displacement is increased by regular displacement increments. After each increment, an image of the region around the crack tip is recorded to allow the analysis of the displacement field.

The images were recorded using a telecentric lens mounted on a 2048 × 2048 pixel digital camera. The lens and the sensor resolution allow for magnification such that the physical size of a pixel is 0.5 µm. The sample surface is prepared by scratching with an ultra-fine glass paper and then by applying carbon black powder. As shown in Fig. 2,

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Schematic description of the sample and experimental setup used to prescribe a displacement controlled, mode I dominated loading. The dimensions of the specimen are $W=50$ mm, $h=20$ mm and the crack length $a=30$ mm, so that $W-a=20$ mm; the thickness is 10 mm. The diameter of the screw is 5 mm and its axis is located 10 mm from the edge of the sample. With this setup, an opening mode I is dominant but due to some a-symmetry in the displacement load, some mode II (shear mode) is also present.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Reference image obtained with a 2048 × 2048 pixel camera. The pixel size is 0.5 µm/pixel.</figcaption></figure>
a random grey-level pattern suitable for DIC analysis is obtained. Note that the area illuminated by the lens is smaller than the sensor resulting in a black region around pixels of the useful area.

### 3. Fracture mechanics analysis

#### 3.1. Displacement field decomposition and crack tip identification

For a 2D elastic body with a semi-infinite straight crack subjected to a mechanical load, an analytical solution for the displacement field was proposed by Williams (1957). The displacement is decomposed over a double series of $\phi_j^n$ solution fields that satisfy stress-free conditions along the crack faces. Two fracture modes (indexed by j) are usually distinguished: mode I (opening mode) and mode II (shear mode). These elementary solutions are indexed by n, they are written as follows:

$$ \phi_I^n(r, \theta) = r^{n/2} \left( \kappa e^{in\theta/2} - \frac{n}{2} e^{i(4-n)\theta/2} + \left(\frac{n}{2} + (-1)^n\right) e^{-in\theta/2} \right) \quad (1) $$

and

$$ \phi_{II}^n(r, \theta) = ir^{n/2} \left( \kappa e^{in\theta/2} + \frac{n}{2} e^{i(4-n)\theta/2} - \left(\frac{n}{2} - (-1)^n\right) e^{-in\theta/2} \right), \quad (2) $$

where i is the pure imaginary number, r the distance to the crack tip and θ the angle with respect to the co-ordinate system aligned with the crack. Hence, adopting complex notations, the position of a material point is defined by z which can be expressed as x+iy or by $re^{i\theta}$. κ is the Kolossov constant, namely, $κ = (3-v)/(1+v)$ for plane stress or $κ = (3-4v)$ for plane strain conditions, where v is the Poisson ratio.

This family of solutions is the appropriate basis to describe the displacement field around the tip of a crack. It contains:

* $n=0$ rigid body translations

* $n=1$ usual asymptotic terms that are weighted by stress intensity factors

* $n=2$ rigid body in-plane rotation and T-stress

* $n > 2$ sub-singular terms, i.e. the corresponding elastic stress field is not singular at the crack tip.

Usually, only $n \ge 0$ terms are considered because $n < 0$ terms exhibit an infinite strain energy density at the crack tip. They are thus not considered in the analysis. An actual displacement field $u$, determined over the previously introduced basis is

$$ u = \sum_{n \ge 0} K_j^n \phi_j^n(z). \quad (3) $$

One may be interested in the inverse problem, namely, from $u$, find the coefficients $K_j^n$. Due to the ill-posedness of this problem, one has to consider only a finite number of terms (in practice, not more than 5-7). Furthermore, not only the $K_j^n$ are unknown, but also the position of the crack tip and the orientation of the crack. This information cannot be accessed directly. However, finding the orientation of the crack from a displacement field (or an image) is usually easy. Eventually, the only remaining unknown is the position of the crack tip along a line defining the crack orientation. Following Hammam et al. (2006), Roux and Hild (2006), or Roux et al. (2009), a recursive formula can be obtained when Eq. (1) is differentiated with respect to the abscissa x:

$$ \frac{\partial \phi_j^n}{\partial x} = \frac{n}{2} \phi_j^{n-2}. \quad (4) $$

Suppose now that the estimated crack tip position is shifted by a small distance d. If the displacement $u$ is now expressed as a function of the Williams' series in the shifted co-ordinate system ($\tilde{z} = z+d$)

$$ u = \sum_{n \ge 0} K_j^n \phi_j^n(z) = \sum_{n,j} \tilde{K}_j^n \phi_j^n(\tilde{z}), \quad (5) $$

using a first order Taylor expansion of $\phi_j^n(\tilde{z})$ and Eq. (4) allows consideration of not only $n \ge 0$ terms but also $n < 0$ terms. Then, one can write the following by considering that the expression of $u$ from Eq. (3) has a zero term for $n = -1$ and mode I:

$$ -\frac{d}{2} \tilde{K}_I^1 + \tilde{K}_I^{-1} = 0. \quad (6) $$

The implication is that using the first super-singular ($n = -1$) term, one can identify the mis-positioning of the crack tip. In practice, the process is iterative: the inverse problem is solved in a least-squares sense, finding $K_j^I$ for mode I, mode II and for n between e.g. -3 to 5. Then, the mis-positioning d of the crack tip is estimated using Eq. (6). The co-ordinate system is progressively corrected until sufficient precision (for the crack tip position) is reached. Because the first order correction only is considered in order to obtain Eq. (6), the relation between the actual shift, d, and its estimation from $-2\tilde{K}_I^{-1}/\tilde{K}_I^1$ is linear. In practice, this linear dependence is observed to occur for quite a large range of crack tip shifts and convergence is usually obtained in 2 or 3 iterations. Note that data from the vicinity of the crack tip are omitted due to the super-singular nature
of $n < 0$ terms. This technique thus provides not only the stress intensity factors $K_I = K_I^1$ and $K_{II} = K_{II}^1$ but also the position of the crack tip along the crack faces.

When material non-linearities occur in the vicinity of the crack tip, the obtained crack tip location is the one generated by the elastic fields that best describe the actual displacement field $\mathbf{u}$, namely the equivalent elastic crack tip.

## 3.2. Extraction from digital images

During the experiment, digital images of the sample surface are recorded. Let $f$ be the grey-level function describing the initial image and $g$ the function for a generic deformed image. Between the two images, the sample has been loaded, resulting in a displacement field $\mathbf{u}$ at the observed surface. The grey-level distribution is thus advected by the displacement field $\mathbf{u}$. This can be written as

$$f(z) = g(z + \mathbf{u}). \qquad (7)$$

Following Roux and Hild (2006), instead of trying to solve the inverse problem of the previous section from a known displacement field, it is solved directly from the digital images. Given $f$ and $g$, $K_I^n$ are found such that

$$f(z) = g \left( z + \sum_{n,j} K_j^n \phi_j^n(z) \right). \qquad (8)$$

The fracture parameters and the crack tip position can thus be extracted directly from the digital images. The inverse problem is solved in a least-squares sense by minimizing the cost function $\eta$ defined by

$$\eta^2 = \int_{ROI} \left[ f(z) - g \left( z + \sum_{n,j} K_j^n \phi_j^n(z) \right) \right]^2 dz, \qquad (9)$$

ROI stands for the region of interest over which the problem is solved. A non-linear least-squares resolution is performed to obtain the best estimates of $K_I^n$ given a position of the crack tip. Then the crack tip shift $d$ is calculated using Eq. (6), the pixel co-ordinates $z$ are updated using the new crack tip position and a new estimate of $K_I^n$ is calculated. Once a sufficiently small value of $d$ is obtained (in practice 1 pixel), the process is stopped and the stress intensity factors and other fracture parameters of the elastic field corresponding to the equivalent elastic crack are stored.

## 3.3. Cohesive zone analysis

As proposed by Bazant and Beissel (1994), the field emanating from a crack with a cohesive zone can be seen as the superposition of the contribution of elastic crack fields whose tips are distributed along the cohesive segment. We will show in the following analysis how a cohesive crack can be described whether in terms of an equivalent elastic crack using the formalism described in the previous section or by using a smeared distribution of elastic cracks. Furthermore, a link between these two descriptions is established. First, the superposition of two cracks is considered. Each crack field has a $k$ singularity. The distance between two crack tips is denoted as $\ell$ and is considered infinitely small compared to the crack length. The resulting displacement field is

$$\mathbf{u} = k\phi_I^1(z) + k\phi_I^2(z+\ell). \qquad (10)$$

Using Eq. (4) and a first order expansion of $\phi_I^1(z+\ell)$, the resulting displacement field now involves the first super-singular function $\phi_I^{-1}$

$$\mathbf{u} = 2k\phi_I^1(z) + \frac{k\ell}{2}\phi_I^{-1}(z). \qquad (11)$$

The resulting stress intensity factor is $2k$ and the first super-singular coefficient is $k\ell/2$. From Eq. (6), the tip of the equivalent elastic crack is translated by $d = \ell/2$. Once the tip of the equivalent elastic crack is translated by the appropriate distance $d$, then the first order correction cancels out and the second order rest is accounted for

$$\mathbf{u} = 2k\phi_I^1\left(z+\frac{\ell}{2}\right) - \frac{k\ell^2}{4} \frac{\partial^2 \phi_I^1}{\partial x^2}\left(z+\frac{\ell}{2}\right). \qquad (12)$$

Again using Eq. (4), the third super-singular field must now be involved. Its amplitude is $-k\ell^2/16$. A remarkable property is that the distance from the actual equivalent crack tip over which the non-linear processes are distributed is

$$\frac{\ell}{2} = \sqrt{\frac{-8K_I^{-3}}{K_I^1}}, \qquad (13)$$

in agreement with Roux and Hild (2006) and Roux et al. (2009). Finally, the first and third super-singular amplitudes of the Williams' series provide a full description of the geometrical features of the non-linear processes at the tip of a cohesive crack. However, the evaluation of the third super-singular term $K_I^{-3}$ from experimental data is extremely difficult (Réthore et al., 2011). Indeed, it requires the estimation of the contribution of a field that varies as $r^{-3/2}$, resulting in a high
sensitivity to noise. From this simple example with two cracks, it is shown that the equivalent crack tip description of a smeared distribution of cracks has its tip shifted inside the cohesive zone. It does not coincide with the end of the cohesive zone. Thus, starting from an initially sharp crack, the formation of a cohesive zone will result in an advance of its equivalent elastic crack tip.

These results are extended by Roux et al. (2009) to an arbitrary smeared tip distribution. The moment of order q is linked to the 1-2q order super-singular coefficient in the Williams series

$$K_I^{1-2q} = \frac{(-1)^{q+1}(2q)!}{(2q-1)2^{2q}(q!)^2} \int s^q k(s) \, ds. \qquad (14)$$

Obviously, if a constant distribution of stress intensity k(s) is considered, the result is the same as for two discrete cracks. However this case is not admissible as the stress field remains singular at the tip of the cohesive zone. The stress intensity distribution must vanish at the tip of the cohesive zone. Then the crudest distribution to be considered is the linear one. In this case, the equivalent elastic crack tip is shifted by one third of the cohesive segment length from the usual tip.

### 3.4. Analysis of the experiment

The sequence of 21 images obtained during the experiment is analyzed with the procedure described above. A ROI of 250 × 200 µm² centered at the presupposed position of the crack tip is analyzed. A circular region of 15 µm in diameter is omitted around the crack tip as well as a band of 20 µm along the crack faces. As an example, Fig. 3 shows the two components of the displacement field for the last image before failure. One can observe that both the displacement orthogonal to the crack and parallel to the crack faces have discontinuities. This means that the experiment cannot be analyzed as a pure mode I experiment. This is confirmed by the variation of the identified stress intensity factors during loading (Fig. 4). At first, $K_{II}$ is positive whereas at the last image it is negative and the ratio $K_{II}/K_I$ is -0.25 which emphasizes that a non-negligible mode mixity is obtained. Furthermore, the maximum value for $K_I$ is 1.1 MPa√m which is in good agreement with reported values for fracture toughness in PMMA (Döll, 1983; Döll and Könczöl, 1990). Fig. 5 shows the variation of the equivalent elastic crack tip position $x_{tip}$ during the experiment. This can be plotted as a function of an equivalent stress intensity factor $K_{eq} = |K_I + iK_{II}|$ as presented in Fig. 6. This plot is equivalent to a resistance-curve and its variation during the loading shows that two phases can be distinguished:

1. Before image 15 (i.e. for $K_{eq} < 0.7$ MPa√m), the position of the equivalent elastic crack remains fixed and no change in $x_{tip}$ is measured.

2. Beyond image 15, a regular advance of the crack tip position is measured, and it increases with $K_{eq}$.

The value of $K_{eq}$ increases up to a maximum value that corresponds to the onset of crack propagation. When the crack starts to grow, the crack tip rapidly escapes the view of the camera so that this point corresponds to the last observation and record of the crack displacement fields. In the R-like representation of Fig. 6, three stages can be distinguished: (i) up to a load level $K_{eq} \approx 0.7$ MPa√m, the crack tip position remains unchanged but (ii) starts to increase for $K_{eq} > 0.7$ MPa√m. This is interpreted as the onset of crazing and subsequent craze widening related to a continuous shift in the equivalent elastic crack tip position up to the maximum $K_{eq} = 1.2$ MPa√m. This latter value is consistent with the measure of the toughness reported in PMMA (Döll, 1983; Döll and Könczöl, 1990; Saad-Gouider et al., 2006).

Furthermore, just before the onset of failure for $K_{eq} = 1.2$ MPa√m, the advance of the equivalent elastic crack tip position is about 10 µm. This variation has to be considered as an indicator for the onset and growth of the failure process zone. It does not correspond to the length of the fracture process zone. The estimation of the process zone length is postponed for the moment and will be derived from the cohesive zone model determined in the following section. In fact, there is no direct link between the advance of the equivalent elastic crack tip and the front of the cohesive zone process zone.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Horizontal and vertical displacement field in pixels (1 pixel is 0.5 µm) for the last image before failure.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Variation of the mode I and mode II stress intensity factors during the experiment.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Variation of equivalent crack tip position $x_{tip}$ during the experiment.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Variation of the equivalent elastic crack tip position with the equivalent stress intensity factor during the experiment.</figcaption></figure>
The results presented in this section allow for the analysis of crazing in PMMA using fracture mechanics concepts. The far-field is analyzed to extract averaged estimates (stress intensity factors, equivalent crack tip position, mode mixity) of the local mechanical state at the crack tip. These estimates will be used to identify a cohesive zone model that represents the failure mechanism.

## 4. Cohesive zone modeling and identification

The goal of this section is to identify a cohesive zone model that describes the mechanics of the failure process for the material under consideration. For this purpose, numerical simulations invoking cohesive zone concepts are carried out. Then, the results are analyzed within the methodology presented in the previous section and an inverse problem is formulated so that cohesive zone parameters can be identified. In this section, a refined model is used at the crack tip, where local information about the mechanical state at the crack tip is required.

### 4.1. Cohesive zone models

The basic idea of cohesive zone models is that the non-linearities related to the fracture mechanism are localized so that the non-linear processes are considered to be confined along an interface (a line in 2D and a surface in 3D). Hence, across a virtual extension of a crack beyond the crack tip, it is supposed that a displacement discontinuity $[\mathbf{u}]$ is allowed and also that a traction $t$ can be transmitted across the crack faces. This traction cancels the stress singularity at the crack tip and the initial boundary value problem is supplemented by a constitutive relation

$$C([\mathbf{u}], t) = 0. \qquad (15)$$

In order to account for mixed-mode situations, the cohesive zone model is formulated using an equivalent traction (e.g. Camacho and Ortiz, 1996)

$$t_{eq} = \sqrt{t_n^2 + \beta t_t^2}, \qquad (16)$$

where the subscripts $n$ and $t$ denote the normal and tangential components of the traction vector and $\beta$ is a parameter controlling the influence of the shearing mode. The procedure is the same for the displacement jump

$$u_{eq} = \sqrt{u_n^2 + u_t^2 / \beta}. \qquad (17)$$

Within this framework, the ratio between the stiffness of the cohesive zone in the tangential direction and in the normal direction is $1/\beta$. In the elastic regime, the ratio between shear and tensile bulk stiffness is $1/2(1+\nu)$. $\beta$ is thus set to $2(1+\nu)$ in order to get the same stiffness ratio in the cohesive zone and in the bulk. The constitutive cohesive zone model now reduces to a relationship between $t_{eq}$ and $u_{eq}$. A generic form for the law to be identified is presented in Fig. 7. The critical traction $t_{cr}$ corresponds to the onset of the failure mechanism (crazing in the case of PMMA), $u_{cr}$ being the critical opening for which a crack nucleates locally and $k_{cz}$ represents the slope of the traction-opening profile during crazing. This model appears similar to the Dugdale formulation (Dugdale, 1960) and it has also been used to describe failure in PMMA by Elices et al. (2002) for instance. It is versatile enough to capture almost any failure process. Unlike the derivations of cohesive zone models found in the literature (Abanto-Bueno and Lambros, 2005; Fedele et al., 2009), the formulation depicted in Fig. 7 is identified from the analysis of the displacement crack tip fields, as presented in the next section.

### 4.2. Identification strategy

As in the previous section, an inverse problem is formulated to extract the parameters of a cohesive law. The cost function that is chosen to be minimized is the quadratic error between the equivalent elastic crack tip position from Section 3.4 and from a numerical simulation of the experiment. For this purpose, the eXtended Finite Element Method

<figure><img src="image_9.png" /><figcaption>Fig. 7. Parametrization of the cohesive law to be identified.</figcaption></figure>
(X-FEM) (Black and Belytschko, 1999; Moës et al., 1999) is used as a convenient tool to introduce displacement discontinuity within a classical finite element discretization. The implementation of the cohesive zone model used within this framework is similar to Remmers et al. (2003). The mesh used for the numerical simulations covers approximately the same area of the image as the ROI in the previous section. It is shown in Fig. 8. The mesh is aligned with the crack orientation and it is refined in the vicinity of the cohesive zone. The elements in the refined region are 2 × 2 µm² so that from the result of Fig. 5, the crack process zone would propagate through several elements. To make the numerical simulations as realistic as possible, a displacement field obtained by finite element DIC (in Eq. (9) the displacement is determined using a finite element basis instead of the Williams' series) is applied along the entire boundary of the mesh. As an illustration, Fig. 9 shows the distribution of the equivalent strain $\epsilon_{eq} = \sqrt{\epsilon} : \epsilon$ on the deformed mesh for different steps in the numerical simulation.

With the proper kinematics being applied to the numerical model through appropriate boundary conditions, the only remaining unknown is the behavior of the cohesive zone. The bulk is assumed to be isotropic linear elastic with Young's modulus *E* and Poisson's ratio *v*, reported in Table 1 (Saad-Gouider et al., 2006).

<figure><img src="image_10.png" /><figcaption>Fig. 8. Finite element mesh for the numerical simulations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Distribution of the equivalent strain $\epsilon_{eq} = \sqrt{\epsilon} : \epsilon$ on the deformed mesh (amplification factor 5) for the steps 6, 11 (top), 16 and 21 (bottom).</figcaption></figure>
<table><caption>**Table 1**

Bulk parameters for PMMA (Saad-Gouider et al., 2006).</caption><thead><tr><td><i>E</i> (MPa)</td><td><i>ν</i></td></tr></thead><tbody><tr><td>5000</td><td>0.32</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Comparison of the evolution of the mode *l* and mode *ll* stress intensity factors for DIC and X-FEM with the identified cohesive law.</figcaption></figure>

After a numerical simulation is run for a given set of cohesive parameters, fracture parameters and crack tip position are extracted from the computed displacement field as described in Section 3.1. The evolution of the crack tip position $x_{tip}^{X-FEM}$ is obtained and compared to the values $x_{tip}$ plotted in Fig. 5. The cohesive zone parameters are then updated with a non-linear least-squares minimization of the following cost function:

$$ \eta_{tip}^2 = \sum_{step} [x_{tip} - x_{tip}^{X-FEM}]^2. \qquad (18) $$

For the present analysis, two parameters are sought: the traction $t_{cr}$ at which crazing initiates and the slope $k_{cz}$ of the curve describing the traction-separation law. In the experiment, the failure process by crazing is captured and the critical opening $u_{cr}$ is first set to infinity. The converged values are

$$ t_{cr} = 135 \text{ MPa} \quad \text{and} \quad k_{cz} = 20 \text{ MPa/µm}. \qquad (19) $$

For this set of parameters, Figs. 10 and 11 compare the stress intensity factors and the evolution of the position of the equivalent elastic crack tip for the numerical simulation and the analysis of Section 3.1. The two analyses agree of course in terms of equivalent crack tip position, because this was the cost function to be minimized but also in terms of stress intensity factors. These results show that with the identified set of cohesive parameters, the numerical model has the same behavior as the actual specimen, from a fracture mechanics point of view. As the cohesive surfaces have been modelled in the numerical calculation, a more detailed analysis of the mechanical state within the cohesive zone is possible. Fig. 12 shows the variation of the displacement jump within the cohesive zone for different steps during the loading. First, normal as well as tangential displacement jumps are obtained in agreement with the mode mixity that has already been measured. The maximum value of the normal opening is 2.4 µm and -0.55 µm for the tangential opening. These values correspond to the displacement jumps at the onset of crack propagation, or the last image recorded before failure.

With the cohesive zone analysis, one can manually identify the location of the cohesive zone tip, where the displacement jump vanishes. This gives the position of the equivalent kinematic crack tip. As shown in Fig. 11, the predicted position of the kinematic crack tip differs from the tip position of the equivalent crack defined in the foregoing section from the analysis of the far field and related stress intensity factors. Fig. 13 shows the traction distribution along the crack faces within the cohesive zone for different steps. Again, mode mixity is observed because normal as well as tangential tractions are obtained. The maximum value for the normal traction is about 130 MPa and -30 MPa in the tangential direction. The location of the maximum traction coincides with the kinematic crack tip and related cohesive zone tip. The variation of the kinematic crack tip position during loading corresponds to the extent of the crazing process zone, which is estimated here (see Fig. 11) to be close to twenty microns at the onset of crack propagation. This estimation is comparable to values reported by Döll and coworkers (≈ 30 µm) (Döll, 1983; Döll and Könczöl, 1990) that were based
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Comparison of the equivalent crack tip position $x_{tip}$ for DIC, X-FEM with the identified cohesive law from a fracture mechanics analysis and by manual pointing on the crack opening profile (kinematic).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Displacement jump distribution in the cohesive zone for different steps.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Traction distribution in the cohesive zone for different steps.</figcaption></figure>
on direct optical interferometry measurements. The present value (≈ 20 µm) is slightly smaller but this is consistent with the difference in critical opening: 2.4 µm compared to 3 µm reported in Döll (1983) and Döll and Könczöl (1990).

Thus, the craze critical opening and the extent of the craze process zone can be extracted from the mixed DIC/X-FEM analysis in addition to the identification of the traction-opening profile. The latter exhibits some decay of the traction with opening, which is similar to a Dugdale profile (Dugdale, 1960). In terms of fracture energy, an energy release rate of 269 J/m² is found from the traction opening profile, which is in good agreement with an estimate, 266 J/m², from the stress intensity factors using the Irwin formula.

**5. Discussion and concluding remarks**

In the literature, there are many studies reporting the identification of cohesive zone models from the analysis of displacement fields obtained by DIC (e.g. Fedele et al., 2009; Shen and Paulino, 2011; Abanto-Bueno and Lambros, 2005; Fuchs and Major, 2011). These examples focus on the traction-opening response of interfaces or adhesive layers and address fracture process zones that share common dimensions: a length of the process zone on the order of hundreds of microns and a critical opening of tens of microns. In such cases, a direct measurement of the displacement jump across the plane of separation is tractable. By calculating the strain fields and with a constitutive law for the surrounding bulk, the derivation of the tractions along the plane of separation is obtained. However, such a direct observation is *not* tractable for fracture processes at smaller scales such as the crazing in PMMA reported here. In this case, the cohesive zone opening is of the order of some microns, close to the limit of resolution of DIC optical techniques, and prohibiting a direct identification of the displacement jump. In addition, full account of the crack tip fields is necessary as the crack is subjected to a *K*-dominated loading. There are numerous cases in which these limiting conditions are present in polymer based materials but also in polymer-polymer assemblies or interfaces involving a polymer bulk. The size of the process zone is even smaller in ceramic/ceramic or metal/ceramic interfaces.

The proposed methodology tackles these configurations, through analysis of the displacement fields in the vicinity of the process zone within a fracture mechanics framework. The development of a fracture process zone induces cohesive tractions that result in a closure stress intensity factor, $K^{closure}$. This latter term is superposed to the one related to the external loading $K^{ext}$, and the total combined stress intensity factor is denoted *K*.

In the present study, a new approach to extract information about the process zone from the analysis of the displacement in the region where the *K*-fields are dominant is proposed. First, a characterization of the load present at the crack in terms of stress intensity factors (for the modes I and II) is necessary. This calculation is performed in accordance with the prescription of Roux and Hild (2006). From the calculated stress intensity factors, a linear elastic problem is considered for which the position of the equivalent elastic crack tip is derived (Hammam et al., 2006; Roux and Hild, 2006). During the loading, the position of the crack tip remains constant until the initiation of the fracture mechanism. Upon initiation, the position of the equivalent elastic crack tip is shifted progressively from its initial location towards the bulk. This shift of the equivalent elastic crack tip provides an experimental signature of the fracture process. The analysis of the local *K*-dominated displacement fields allows an estimation of the mode mixity and ultimately the identification of the normal *and* the tangential traction opening response (see Fig. 13).

The experimental crack tip displacement fields are employed for the identification of a cohesive zone model for PMMA. The traction-opening profile is similar to the one reported by Elices et al. (2002), although they considered test configurations only in mode I. Elices et al. (2002) identified the cohesive properties from the best fit comparison of the predicted macroscopic force-displacement obtained from a finite elements analysis and their experimental data. This methodology focuses on the crack propagation stage and misses local information regarding how realistic the cohesive model is. In the present work, local information such as the local displacement fields are analyzed within a fracture mechanics approach that captures the initiation of the failure mechanism, and the growth of the cohesive process zone up to the genesis of a crack. The proposed approach requires an initial natural crack, with careful machining.

It is worth noting that the present identification strategy does not require determination of the displacement fields at a length scale comparable to that of the cohesive opening profile but instead uses the information contained in the *K*-displacement fields in the vicinity of the crack. In the present case, the area observed by DIC is approximately 1 × 1 mm² and the cohesive zone characteristics are some microns in maximum opening and tens of microns long. The region along the crack plane is masked to ensure an optimized correlation of the 'far' displacement fields, 'far' from the failure process zone.

The methodology results in an identification of the cohesive zone model that underlies crazing in the case of PMMA with realistic parameter magnitudes for the tractions, traction-opening profile, and estimation of the maximum cohesive opening corresponding to the nucleation of a crack locally ($\delta^{cr} = 2.4$ µm). Three observations can be pointed out related to the three characteristic stages of crazing. Initiation is detected from the beginning in the shift of the position of the equivalent elastic crack during loading, with a local value of 130 MPa for the traction. The traction-opening profile identified during craze thickening is found to decay slightly (as in Elices et al., 2002) while a profile with a constant traction is suggested in Döll (1983), Döll and Könczöl (1990) and Saad-Gouider et al. (2006). There are also studies of fracture in PMMA that suggest some softening in the traction-separation profile (Murphy and Ivankovic, 2005) or a combination of softening followed by hardening (Seelig, 2008). Although no general agreement is found in the
traction-separation profile, the order of magnitude of the traction during crazing is comparable, the slight differences in the profile still suggest that an approximately constant or decaying traction-opening is realistic.

The third feature extracted from the present identification methodology concerns the estimation of the cohesive maximum opening $\delta^{cr} = 2.4 \ \mu\text{m}$ along the normal direction. This estimation is very close to that reported by Döll and co-workers (Döll, 1983; Döll and Könczöl, 1990): $2.5-3 \ \mu\text{m}$ for PMMA. To get this measurement, Döll and co-workers developed a sophisticated interferometry set-up which required steady-state crack growth for the observation of the interference patterns along the cohesive zone and the crack. The displacement jump was then derived from the analysis of the craze wedge and the interfringe distance. Although a direct observation of the displacement jump is not possible optically due to the dimensions of the opening profile under consideration, the identification of the cohesive model provides, not a direct, but a valuable estimation, and correct order of magnitude for the critical opening $\delta^{cr}$, in a simpler way. In addition, a tangential critical opening is extracted in the present identification by a tangential separation profile with a maximum traction at the onset of crazing followed by a continuous decay during crazing. This tangential response originates from the mixed-mode loading but also evidences some tangential stiffness of the web of craze fibrils. Such a stiffness has been addressed in detailed studies of the craze fibrils in Hui et al. (1992) and Sha et al. (1995) and within a cohesive surface methodology by Tijssens et al. (2000) who formulated a rate-dependent cohesive zone model with explicit contributions to the process zone from the normal and tangential components. The tangential component is not reported in the calibration of the cohesive zone model for crazing in PMMA by Saad-Gouider et al. (2006) who considered pure mode I, but also because mixed-mode is not easy to control experimentally. The present analysis shows that mode-mixity can nevertheless be evaluated. The present methodology may provide some insight regarding the influence of the tangential response in cohesive zone models. The contribution of the tangential mode to fracture and to the energy release rate is modest here, as the dominant loading is mode I. However, its influence could be more important for a larger mode mixity or in configurations where a compressive normal component is combined with a tangential one.

In conclusion, the methodology and results presented in this study for the identification of a cohesive zone model that describes the fracture mechanism in PMMA is summarized as follows:

* A fracture mechanics framework has been used to analyze the crack tip displacement fields of a natural crack loaded predominantly in mode I with some mode mixity. The crack tip fields contain information about the prescribed external load and the closure when a cohesive process zone appears. The loading is analyzed in terms of variation of stress intensity factors $K_I$ and $K_{II}$ with respect to the position of the equivalent elastic crack tip.

* The location of the equivalent crack tip during loading is used to detect the onset of the failure mechanism, and the growth of the cohesive process zone up to crack propagation.

* The identification of a cohesive zone model in a finite elements analysis is to capture the variation of the mode I and II stress intensity factors with respect to the position of the equivalent crack tip location.

* The parameters obtained from the identification procedure are consistent with available data reported in the literature for the normal components of the cohesive model. The methodology also provides information about the tangential response, however further experimental measurements are needed for validation.

The methodology for the identification of cohesive zone models can be used for other polymer materials for which cohesive model descriptions are not yet available. The methodology can also be applied to other materials provided that the displacement field information is accessible experimentally. The present methodology could be extended to interface problem. This is the topic of continuing studies.

## References

Abanto-Bueno, J., Lambros, J., 2002. Investigation of crack growth in functionally graded materials using digital image correlation. Eng. Fract. Mech. 69, 1695-1711.

Abanto-Bueno, J., Lambros, J., 2005. Experimental determination of cohesive failure properties of a photodegradable copolymer. Exp. Mech. 45, 144-152.

Abraham, F., 2001. The atomic dynamics of fracture. J. Mech. Phys. Solids 49, 2095-2111.

Argon, A., Hannoosh, J., 1977. Initiation of crazes in polystyrene. Philos. Mag. 36, 1116-1195.

Barenblatt, G., 1962. The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech. 7, 55-126.

Bazant, Z., Beissel, S., 1994. Smeared-tip superposition method for cohesive fracture with a rate effect and creep. Int. J. Fract. 65, 277-290.

Besnard, G., Hild, F., Roux, S., 2006. 'Finite-element' displacement fields analysis from digital images: application to portevin-le châtelier bands. Exp. Mech. 46, 789-803.

Black, T., Belytschko, T., 1999. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 45, 601-620.

Camacho, G., Ortiz, M., 1996. Computational modeling of impact damage in brittle materials. Int. J. Solids Struct. 33, 2899-2938.

Döll, W., 1983. Optical interference measurements and fracture mechanics analysis of crack tip craze zones. Adv. Polym. Sci. 52-53, 105-168.

Döll, W., Könczöl, L., 1990. Micromechanics of fracture: optical interferometry of crack tip craze zone. Adv. Polym. Sci. 91-92, 138-214.

Dugdale, D., 1960. Yielding of steel sheets containing slits. J. Mech. Phys. Solids 8, 100-104.

Elices, M., Guinea, G., Gomez, J., Planas, J., 2002. The cohesive zone model: advantages, limitations and challenges. Eng. Fract. Mech. 69, 137-163.

Estevez, R., Tijssens, M., Giessen, E.V.D., 2000. Modelling of the competition between shear yielding and crazing in glassy polymers. J. Mech. Phy. Solids 48, 2585-2617.

Fedele, R., Hild, B.R.F., Roux, S., 2009. Identification of adhesive properties in glare assemblies using digital image correlation. J. Mech. Phys. Solids 57, 1003-1016.

Fuchs, P., Major, Z., 2011. Experimental determination of cohesive zone models for epoxy composites. Exp. Mech. 51, 779-786.
Gall, K., Horstemeyer, M., Schilfgaarde, M.V., Baskes, M., 2000. Atomistic simulations on the tensile debonding of an aluminum-silicon interface. J. Mech. Phys. Solids 48, 2183-2212.

Grégoire, D., Maigre, H., Morestin, F., et al., 2009. New experimental techniques for dynamic crack localization. Eur. J. Comput. Mech. 18, 255-283.

Hammam, R., Hild, F., Roux, S., 2006. Stress intensity factor gauging by digital image correlation: application in cyclic fatigue. Strain 43, 181-192.

Hong, S., Kim, K., 2003. Extraction of cohesive-zone laws from elastic far-fields of a cohesive crack tip: a field projection method. J. Mech. Phys. Solids 51, 1267-1286.

Hui, C.Y., Ruina, A., Creton, C., Kramer, E.J., 1992. Micromechanics of crack growth into a craze in a polymer glass. Macromolecules 25, 3948-3955.

Kausch, H. (Ed.), 1983. Crazing in Polymers, vol. 1, Advances in Polymer Science. Springer, pp. 52-53.

Kausch, H. (Ed.), 1990. Crazing in Polymers, vol. 2, Advances in Polymer Science. Springer Verlag, pp. 91-92.

Kramer, E., 1983. Microscopic and molecular fundamentals of crazing. Adv. Polym. Sci. 52-53, 1-56.

Kramer, E., Berger, L., 1990. Fundamental processes of craze growth and fracture. Adv. Polym. Sci. 91/92, 1-68.

Kubair, D., Ciacchi, D.C., Spearing, S., 2009. Multiscale mechanics modeling of direct silicon wafer bonding. Scr. Mater. 60, 1125-1128.

McNeill, S., Peters, W., Sutton, M., 1987. Estimation of stress intensity factor by digital image correlation. Eng. Fract. Mech. 28, 101-112.

Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 46, 133-150.

Murphy, N., Ivankovic, A., 2005. The prediction of dynamic fracture evolution in pmma using a cohesive zone model. Eng. Fract. Mech. 72, 861-875.

Needleman, A., 1987. A continuum model for void nucleation by inclusion debonding. J. Appl. Mech. 54, 525.

Nguyen, T., Hall, S., Vacher, P., Viggiani, G., 2011. Fracture mechanisms in soft rock: identification and quantification of evolving displacement discontinuities by extended digital image correlation. Tectonophysics 503, 117-128.

Oxborough, R., Bowden, P., 1973. A general critical-strain criterion for crazing in amorphous glassy polymers. Philos. Mag. 28, 547-559.

Poissant, J., Barthelat, F., 2010. A novel subset splitting procedure for digital image correlation on discontinuous displacement fields. Exp. Mech. 50, 353-364.

Remmers, J., de Borst, R., Needleman, A., 2003. A cohesive segments method for the simulation of crack growth. Comput. Mech. 31, 69-77.

Réthoré, J., Gravouil, A., Morestin, F., Combescure, A., 2005. Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral. Int. J. Fract. 132, 65-79.

Réthoré, J., Hild, F., Roux, S., 2007a. Extended digital image correlation with crack shape optimization. Int. J. Numer. Methods Eng. 73, 248-272.

Réthoré, J., Hild, F., Roux, S., 2007b. Shear-band capturing using a multiscale extended digital image correlation technique. Comput. Methods Appl. Mech. Eng. 196, 5016-5030.

Réthoré, J., Roux, S., Hild, F., 2011. Optimal and noise-robust extraction of fracture mechanics parameters from kinematic measurements. Eng. Fract. Mech. 78, 1827-1845.

Rice, J., 1980. The Mechanics of Earthquake Rupture. Italian Physical Society and North- Holland, Amsterdam, The Netherlands, pp. 555-649.

Roux, S., Hild, F., 2006. Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches. Int. J. Fract. 140, 141-157.

Roux, S., Réthoré, J., Hild, F., 2009. Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks. J. Phys. D Appl. Phys. 42.

Saad-Gouider, N., Estevez, R., Olagnon, C., Séguéla, R., 2006. Calibration of a viscoplastic cohesive zone for crazing in pmma. Eng. Fract. Mech. 73, 2503-2522.

Seelig, T., 2008. Computational modeling of deformation mechanisms and failure in thermoplastic multilayer composites. Compos. Sci. Technol. 68, 1198-1208.

Sha, Y., Hui, C.Y., Ruina, A., Kramer, E., 1995. Continuum and discrete modeling of craze failure at a crack tip in a glassy polymer. Macromolecules 28, 2450-2459.

Shen, B., Paulino, G., 2011. Direct extraction of cohesive fracture properties from digital image correlation: a hybrid inverse technique. Exp. Mech. 51, 143-163.

Sternstein, S., Myers, F., 1973. Yielding of glassy polymers in the second quadrant of principal stress space. J. Macromol. Sci. Phys. B 8, 539-571.

Sternstein, S., Ongchin, L., Silverman, A., 1968. Inhomogeneous deformation and yielding of glasslike high polymers. Appl. Polym. Symp. 7, 175-199.

Sun, Y., Pang, J., Wong, C., Su, F., 2005. Finite-element formulation for a digital image correlation method. Appl. Opt. 44, 7357-7363.

Sutton, M., Cheng, M., Peters, W., Chao, Y., McNeill, S., 1986. Application of an optimized digital image correlation method to planar deformation analysis. Image Vision Comput. 4, 143-150.

Tijssens, M., der Giessen, E.V., Sluys, L., 2000. Modelling of crazing using a cohesive surface methodology. Mech. Mater. 32, 19-35.

Williams, J., Pavan, A. (Eds.), 2003. Fracture of Polymers, Composites and Adhesives, vol. II. Elsevier.

Williams, M., 1957. On the stress distribution at the base of a stationary crack. ASME J. Appl. Mech. 24, 109-114.

Xu, X., Needleman, A., 1994. Numerical simulations of fast crack growth in brittle solids. J. Mech. Phys. Solids 42, 1397-1434.

Yamakov, V., Saether, E., Phillips, D., Glaessgen, E., 2006. Molecular-dynamics simulation-based cohesive zone representation of intergranular fracture processes in aluminum. J. Mech. Phys. Solids 54, 1899-1928.