INVESTIGATION OF TRANSLAMINAR FRACTURE IN FIBREREINFORCED COMPOSITE LAMINATES—APPLICABILITY OF LINEAR
ELASTIC FRACTURE MECHANICS AND COHESIVE-ZONE MODEL
By
Fang Hou

A DISSERTATION
Submitted to
Michigan State University
In partial fulfillment of the requirements
for the degree of
Engineering Mechanics - Doctor of Philosophy
2014

ABSTRACT
INVESTIGATION OF TRANSLAMINAR FRACTURE IN FIBREREINFORCED COMPOSITE LAMINATES—APPLICABILITY OF LINEAR
ELASTIC FRACTURE MECHANICS AND COHESIVE-ZONE MODEL
By
Fang Hou
With the extensive application of fiber-reinforced composite laminates in industry, research on
the fracture mechanisms of this type of materials have drawn more and more attentions. A
variety of fracture theories and models have been developed. Among them, the linear elastic
fracture mechanics (LEFM) and cohesive-zone model (CZM) are two widely-accepted fracture
models, which have already shown applicability in the fracture analysis of fiber-reinforced
composite laminates. However, there remain challenges which prevent further applications of the
two fracture models, such as the experimental measurement of fracture resistance.
This dissertation primarily focused on the study of the applicability of LEFM and CZM for the
fracture analysis of translaminar fracture in fibre-reinforced composite laminates. The research
for each fracture model consisted of two sections: the analytical characterization of crack-tip
fields and the experimental measurement of fracture resistance parameters.
In the study of LEFM, an experimental investigation based on full-field crack-tip displacement
measurements was carried out as a way to characterize the subcritical and steady-state crack
advances in translaminar fracture of fiber-reinforced composite laminates. Here, the fiberreinforced composite laminates were approximated as anisotropic solids. The experimental
investigation relied on the LEFM theory with a modification with respect to the material
anisotropy. Firstly, the full-field crack-tip displacement fields were measured by Digital Image
Correlation (DIC). Then two methods, separately based on the stress intensity approach and the

energy approach, were developed to measure the crack-tip field parameters from crack-tip
displacement fields. The studied crack-tip field parameters included the stress intensity factor,
energy release rate and effective crack length. Moreover, the crack-growth resistance curves (Rcurves) were constructed with the measured crack-tip field parameters. In addition, an error
analysis was carried out with an emphasis on the influence of out-of-plane rotation of specimen.
In the study of CZM, two analytical inverse methods, namely the field projection method
(FPM) and the separable nonlinear least-squares method, were developed for the extraction of
cohesive fracture properties from crack-tip full-field displacements. Firstly, analytical
characterizations of the elastic fields around a crack-tip cohesive zone and the cohesive variables
within the cohesive zone were derived in terms of an eigenfunction expansion. Then both of the
inverse methods were developed based on the analytical characterization. With the analytical
inverse methods, the cohesive-zone law (CZL), cohesive-zone size and position can be inversely
computed from the cohesive-crack-tip displacement fields. In the study, comprehensive
numerical tests were carried out to investigate the applicability and robustness of two inverse
methods. From the numerical tests, it was found that the field projection method was very
sensitive to noise and thus had limited applicability in practice. On the other hand, the separable
nonlinear least-squares method was found to be more noise-resistant and less ill-conditioned.
Subsequently, the applicability of separable nonlinear least-squares method was validated with
the same translaminar fracture experiment for the study of LEFM.
Eventually, it was found that the experimental measurements of R-curves and CZL showed a
great agreement, in both of the fracture energy and the predicted load carrying capability. It thus
demonstrated the validity of present research for the translaminar fracture of fiber-reinforced
composite laminates.

To my parents, for their love, support and encouragement

iv

ACKNOWLEDGMENTS

I would like to take this opportunity to express my deepest gratitude to all the people who
have guided and supported me for my Ph.D. program.
First and foremost, I would like to thank my parents, Meiyu Tian and Yichuan Hou, who have
raised me and given me the foundation to be the person who I am today. Without their support I
would not be able to come to the United States to pursue my Ph.D. They have always been
supportive and showing confidence in me. Their encouragement is the biggest support for me to
go through all the struggling time during the past five years.
I would like to express my gratitude to my advisor, Dr. Soonsung Hong. I really appreciate the
opportunity that he provided for me to pursue my Ph.D. degree at Michigan State University. I
also thank him for his insightful guidance, which helped me develop the independent thinking
capability and research skills.
I express my deep appreciation to Professor Xinran Xiao, who kindly served as the chair of
my committee after Dr. Hong left Michigan State University. It is through her generous support,
invaluable suggestions and warm encouragement that I could finally accomplish this work.
Appreciations are given to Professor Dahsin Liu, Professor Arjun Tekalur and Professor
Parviz Soroushian, who kindly served in my dissertation committee, for their time and valuable
suggestions.
Special thanks are given to the Department of Mechanical Engineering for providing
continuous funding support for my education and research. Appreciation extends to my
officemates, classmates and friends who contribute to this work, for their support and help.

v

TABLE OF CONTENTS

LIST OF TABLES ...................................................................................................................... .viii
LIST OF FIGURES ....................................................................................................................... ix
Chapter 1............. ............................................................................................................................ 1
Introduction......... ............................................................................................................................ 1
1.1 Overview of fracture models for fiber-reinforced composite liminates ........................ 1
1.2 Objectives and scopes .................................................................................................... 6
1.3 Structure of dissertation ................................................................................................. 7
Chapter 2............. ............................................................................................................................ 9
Literature Review............................................................................................................................ 9
2.1 Application of linear elastic fracture mechanics in fiber-reinforced composite
laminates ................................................................................................................................ 9
2.2 Cohesive-zone model ................................................................................................... 16
2.3 Digital Image Correlation (DIC) .................................................................................. 24
Chapter 3............. .......................................................................................................................... 30
R-curve Characterization of Translaminar Crack Growth in Cross-ply Composite Laminates
using Digital Image Correlation.................................................................................................... 30
3.1 Introduction .................................................................................................................. 30
3.2 Estimation of the crack-tip field parameters in anisotropic solids ............................... 31
3.2.1 Asymptotic expansion of crack-tip fields in anisotropic solids ....................... 31
3.2.2 Least-squares method....................................................................................... 35
3.2.3 Conservation integrals ..................................................................................... 42
3.3 Experimental investigation ......................................................................................... 44
3.3.1 Experimental procedure ................................................................................... 44
3.3.2 Experimental results......................................................................................... 47
3.3.2.1 Load-CMOD curve and full-field displacements ................................... 47
3.3.2.2 Results of least-squares method .............................................................. 50
3.3.2.3 Results of conservation integrals ............................................................ 58
3.3.2.4 Comparison of the results of least-squares and conservation integrals .. 63
3.3.2.5 R-curve behaviors ................................................................................... 65
3.4 Discussion .................................................................................................................... 68
3.5 Summary ...................................................................................................................... 71
Chapter 4............. .......................................................................................................................... 72
Inverse Extraction of Cohesive-zone Law for Fiber-reinforced Composite Laminates Part I:
Numerical Analysis ....................................................................................................................... 72
4.1 Introduction .................................................................................................................. 72
4.2 An eigenfunction expansion of the cohesive crack-tip fields in anisotropic solids ..... 73
4.3 Extraction of CZLs in anisotropic elastic solids .......................................................... 75
vi

4.3.1 Field projection method ................................................................................... 76
4.3.1.1 Interaction J-integral ............................................................................... 76
4.3.1.2 Expansions of the elastic-fields in terms of orthogonal polynomial series
............................................................................................................................... 77
4.3.1.3 Traction probing and separation probing ................................................ 79
4.3.2 Least-squares method....................................................................................... 82
4.3.3 Determination of cohesive-zone size and position .......................................... 82
4.4 Numerical test with synthetic data ............................................................................... 86
4.4.1 Discription of numerical test ............................................................................ 87
4.4.1.1 Cohesive-zone laws and material property ............................................ 87
4.4.1.2 Procedure of numerical test ................................................................... 90
4.4.2 Results of test ................................................................................................... 95
4.4.2.1 Determination of CZLs with known cohesive-zone size and position .... 95
4.4.2.2 Determination of CZLs with unknown cohesive-zone size and position
............................................................................................................................. 105
4.4.3 Guidelines for experiments ............................................................................ 113
4.5 Summary .................................................................................................................... 115
Chapter 5............. ........................................................................................................................ 117
Inverse Extraction of Cohesive-zone Laws for Fiber-reinforced Composite Laminates Part II:
Experimental Validation ............................................................................................................. 117
5.1 Introduction ................................................................................................................ 117
5.2 Experimental procedure ............................................................................................. 118
5.2.1 Experimental setups ....................................................................................... 118
5.2.2 Data reduction scheme ................................................................................... 119
5.3 Experimental results................................................................................................... 121
5.3.1 Displacement fields ........................................................................................ 121
5.3.2 Cohesive-zone variables and cohesive-zone laws ......................................... 125
5.4 Verification of the experimental inverse solutions .................................................... 129
5.5 Summary .................................................................................................................... 133
Chapter 6............. ........................................................................................................................ 134
Conclusions and Future Work .................................................................................................... 134
6.1 Significance and contributions ................................................................................... 134
6.2 Recommendations for future work ........................................................................... 136
APPENDICES............. ............................................................................................................... 139
Appendix A Source Code of 2D Digital Image Correlation (DIC) Software .................... 140
Appendix B Derivation of the Asymptotic Expansion of Semi-infinite Crack Problem in
Anisotropic Solids............................................................................................................... 150
Appendix C Elastic Matrices for Orthotropic Solids under Plane Stress Condition ......... 152
Appendix D Additional Results of Translaminar Fracture Tests with Fiber-reinforced
Composite Laminates.......................................................................................................... 154
Appendix E Derivation of the Eigenfunction Expansion of Cohesive Crack-tip fields in
Anisotropic Solids............................................................................................................... 157
BIBLIOGRAPHY ....................................................................................................................... 160
vii

LIST OF TABLES

Table 3.1 Elastic properties of the composite laminate, Cyply 1002 ....................................... 44
Table 4.1 Coefficients t n and b n (n  0,1, 2) for the construction of different CZLs ................. 89
Table 4.2 Normalized root-mean-square (NRMS) errors of the extracted CZLs for different
inverse distances and noise levels .......................................................................................... 112
Table D.1 Elastic properties of composite laminates used in the translaminar fracture tests . 154

viii

LIST OF FIGURES

Figure 1.1 Micro-scale fracture mechanisms of fiber-reinforced composites .......................... 2
Figure 2.1 Cohesive-zone model ............................................................................................ 16
Figure 2.2 Cohesive-zone law................................................................................................. 17
Figure 2.3 A typical image used in DIC with random gray intensity distribution ................. 25
Figure 2.4 A schematic illustration of the basic principles of DIC ........................................ 26
Figure 3.1 Crack-tip coordinate system of a semi-infinite crack problem. ............................ 33
Figure 3.2 An arbitrary Cartesian coordinate system x  and y  with the crack-tip location
( x0 , y0 ) ............. ......................................................................................................................... 37
Figure 3.3 Flowchart of the separable nonlinear least-squares method for the determination of
coefficients of asymptotic expansion and effective crack-tip position. .................................... 39
Figure 3.4 Schematic diagrams of the in-plane rigid body motions of specimen: (a) in-plane
rigid body translation (b) in-plane rigid body rotation. ............................................................ 40
Figure 3.5 Schematic diagrams of the out-of-plane rigid body motions of specimen: (a) outof-plane translation (b) out-of-plane rotation with respect to y axis (c) out-of-plane rotation
with respect to x axis................................................................................................................. 41
Figure 3.6 Integral paths of conservation integrals. ................................................................ 42
Figure 3.7 (a) Configuration of ECT (ESET) specimen (b) Speckle patterns for DIC
measurement (c) Unpainted side of the specimen after test ..................................................... 46
Figure 3.8

Load versus Crack Mouth Opening Displacement (CMOD) Plot........................ 48

Figure 3.9 An example set of displacement fields measured by DIC: (a) U field (b) V field 49
Figure 3.10 Variation of the (a) stress intensity factor and (b) effective crack length obtained
by the least-squares method with respect to the number of terms in crack-tip asymptotic
expansion........ .......................................................................................................................... 51
Figure 3.11 Variation of the (a) stress intensity factors and (b) effective crack length obtained
by the least-squares method with respect to the half length of interior deleted area in the data
domain, which is represented by W .......................................................................................... 52
Figure 3.12 Comparison of the V-field contours measured by DIC and regenerated from leastsquares method at the moment when CMOD is 0.30 mm. The results of least-squares methods
ix

are obtained (a) without the out-of-plane rotation term and (b) with the out-of-plane rotation
term................. .......................................................................................................................... 54
Figure 3.13 Variation of the (a) Mode I stress intensity factor KI and (b) effective crack length
estimated by the least-squares method with respect to CMOD ................................................ 57
Figure 3.14 Effect of radius of integral path on the estimation of the (a) J-integral and (b)
effective crack length obtained by conservation integrals ........................................................ 60
Figure 3.15 Variation of the (a) J-integral and (b) effective crack length obtained by
conservation integrals with respect to CMOD .......................................................................... 62
Figure 3.16 Comparison of the (a) mode I stress intensity factor and (b) effective crack length
estimated by different methods ................................................................................................. 63
Figure 3.17 R-curves in terms of stress intensity factor, determined by the least-squares
method and conservation integrals............................................................................................ 66
Figure 3.18 R-curves in terms of strain energy release rate, determined by the least-squares
method and conservation integrals............................................................................................ 67
Figure 3.19 Fracture mechanisms in the translaminar fracture of fiber-reinforced composites 68
Figure 3.20 Comparison of (a) mode I and (b) mode II stress intensity factors estimated by
least-squares method and interaction J-integral ........................................................................ 70
Figure 4.1 Schematic diagram of semi-infinite crack with a cohesive zone of size 2c ............ 74
Figure 4.2 Interaction J-integral paths surrounding a cohesive zone at the far field  , and
along the cohesive zone surface 0 .......................................................................................... 77
Figure 4.3 Flowchart of the separable nonlinear least-squares method for the determination of
cohesive-zone law, cohesive-zone size and position ................................................................ 84
Figure 4.4 Three representative CZLs used in the numerical test ............................................ 87
Figure 4.5 Characteristic parameters of input data fields and displacement measurement noise
for the convex upward CZL (Law I): (a) Horizontal analytical displacement field U
(perpendicular to the loading direction) (b) Vertical analytical displacement field V (along the
loading direction) (c) Horizontal analytical displacement field U with the NSTD of noise
being 1 10 3 (d) Vertical analytical displacement field V with the NSTD of noise being 1 10 3
................................................................................................................................................... 93
Figure 4.6 Comparison between the inverse solutions and exact values of the cohesive-zone
variables and cohesive-zone laws: (a) traction (b) separation gradient (c) separation (d)
cohesive-zone law.. ................................................................................................................... 96

x

Figure 4.7 Influence of inverse distance on the accuracy of inverse solutions: (a) the CZLs
extracted by the field projection method (FPM) (b) the CZLs extracted by the linear leastsquares method (LSM) .............................................................................................................. 97
Figure 4.8 Relative errors of the eigenfunction expansion coefficients with respect to the
inverse distance ......................................................................................................................... 98
Figure 4.9 Relative errors of the eigenfunction expansion coefficients extracted from various
data sets with different intervals between neighbouring data points ...................................... 100
Figure 4.10 Relative errors of the eigenfunction expansion coefficients extracted from various
data sets with different noise levels: (a) interval between two neighbouring points is 5 pixels
(b) interval between two neighbouring points is 20 pixels ..................................................... 103
Figure 4.11 Searching paths of the separable nonlinear least-square method on the contour of
objective function: (a) for convex upward CZL (Law I) (b) for linear CZL (Law II) (c) for
concave upward CZL (Law III) .............................................................................................. 108
Figure 4.12 Estimated cohesive-zone size and positions, and extracted CZLs with different
noise levels.......... ................................................................................................................... 110
Figure 5.1 Sketch of the cohesive zone behavior during fracture process: the cohesive traction
and separation distribution profiles, the corresponding equivalent LEFM crack-tip opening
profile and stress distribution; and the different stages of crack growth ................................ 120
Figure 5.2 Load versus Crack Mouth Opening Displacement (CMOD) ............................... 123
Figure 5.3 R-curves in terms of stress intensity factor, determined by the nonlinear leastsquares method........................................................................................................................ 123
Figure 5.4 Displacement fields measured at different moments shown in Figure 5.2 and 5.3
during the steady-state crack growth: (a) at the moment I (b) at the moment II (c) at the
moment III....... ....................................................................................................................... 124
Figure 5.5 Cohesive-zone traction t, separation gradient b, separation δ, and CZL (relationship
between the cohesive-zone traction and separation) extracted from the data set obtained at the
moment II........ ........................................................................................................................ 126
Figure 5.6 Cohesive-zone sizes and positions extracted from the data set obtained at the
moment II........ ........................................................................................................................ 126
Figure 5.7 Cohesive-zone laws extracted at different time.................................................... 128
Figure 5.8 Cohesive-zone sizes and positions estimated at different time. ........................... 128
Figure 5.9 Finite element model of ECT specimen with inserted cohesive elements ........... 131

xi

Figure 5.10 Comparison of the load versus crack mouth opening displacement (CMOD)
curves between the experiment and finite element simulation ............................................... 131
Figure D.1 Images of specimens captured at the end of translaminar fracture tests: (a) UD00
(b) UD45 (c) UD90 (d) CP00 (e) CP45 (f) CP90 ................................................................... 155
Figure D.2 R-curves in terms of stress intensity factor, for specimen UD90 ........................ 156

xii

Chapter 1
Introduction
1.1 Overview of fracture models for fiber-reinforced composite
laminates
Fiber-reinforced composite laminates have already emerged as an important category of
structural materials, which have been widely applied in automotive industry, aerospace,
transportation, infrastructure and civil engineering. Like other structural materials, the fracture
analysis and safe-operating life prediction of fiber-reinforced composite laminates are critical in
its industrial applications, particularly when there are crack-like defects in the structures (from
where the fracture processes usually emanate). It thus necessitates the research efforts to
understand the fracture mechanisms of this type of materials. In the past decades, plenty of
research attempts have been made in this field. However, the overall progress in the development
of fracture models and experimental techniques for the study of fracture mechanisms of fiberreinforced composite laminates is still in a relatively early stage.
The lag in the research progress of the fracture mechanisms is mainly due to the complexity of
material characteristics of fiber-reinforced composite laminates, as compared to conventional
industrial materials (e.g. metallic materials). Most conventional industrial materials can be
approximated to be homogeneous or even isotropic. Hence the fracture behaviors of these
materials can be characterized through the linear elastic fracture mechanics or elastic plastic
fracture mechanics, both of which have considerably simple theoretical frames and fracture
criteria. Conversely, the heterogeneous structures of fiber-reinforced composite laminates
1

introduce a variety of mechanisms into their fracture processes, which cannot be easily
characterized with simple fracture models. A schematic diagram of the micro-scale fracture
mechanisms around the crack-tip in fiber-reinforced composite is shown in Figure 1.1. Various
micro-scale fracture mechanisms can be seen, including fiber bridging and pullout, micro-cracks
in the matrix and fiber/matrix debonding. Moreover, the ply-level fracture characteristics of
fiber-reinforced composites laminates can be classified as interlaminar, intralaminar and
translaminar, each of which involves one or several micro-scale fracture mechanisms.

Figure 1.1 Micro-scale fracture mechanisms of fiber-reinforced composites

Various fracture mechanisms of fiber-reinforced composites laminates present a fracture
mechanics researcher with a very complicated situation. Instead of a widely-accepted fracture
theory, numerous fracture models have been proposed for fiber-reinforced composite laminates
since the 1970s. From the multi-scale point of view, these fracture models can be categorized

2

into three groups: the micromechanical models [1]–[7], the cohesive-zone model (CZM) [8]–[10]
and the linear elastic fracture mechanics (LEFM) theory with anisotropic modification [11]–[13].
Each group of fracture models will be briefly introduced below, with their advanteges and
disadvanteges.
Firstly, it is known that the micro-scale fracture mechanisms within a crack-tip fracture process
zone can cause local load redistributions and stress concentrations. The micromechanical fracture
models incorporate these effects into the fracture analysis. In the past several decades, a variety
of studies with different types of micromechanical models have been carried out for fiberreinforced composite laminates [1]. Recently, a comprehensive review paper regarding these
work has been published [14]. Generally, the micromechanical models are attractive for their
capability to predict the micro-scale fracture behaviours of fiber-reinforced composite laminates
based on the properties of fibres and matrix. However, the lag in the development of
computationally-efficient multi-scale modeling techniques prevents the micromechanical
fracture models from being used in the fracture analysis of macro-scale structures.
In recent years, the cohesive zone model (CZM) [8]–[10] has drawn more and more attentions.
It has been widely applied for the fracture analysis of fiber-reinforced composite laminates, as
well as other quasi-brittle materials such as concrete. Compared to the micromechanical fracture
models, the CZM does not consider the separate behaviour of fibre and matrix during the fracture
process. Instead, the bulk material is assumed to be homogeneous, while the crack-tip fracture
process zone is simplified as a cohesive zone. The local fracture process within the cohesive
zone is characterized by the relationship between the cohesive-zone traction, which is in
equilibrium with the surrounding stress fields, and the cohesive-zone separation, which is
compatible with the surrounding deformation fields. This relationship constitutes the cohesive-

3

zone law (CZL). One important advantage of the CZM is that it is especially suitable for the
numerical modeling of fracture due to its balanced conceptual simplicity against fracture
predictive capability, as well as the computational efficiency. A feature of fracture simulation
with CZM is that the functional shape of the CZL can significantly influence the numerical
simulation results. Therefore, a good estimation of the CZL is critical to achieve truly predictive
simulations. However, although a variety of methods have been developed, the identification of
the details of CZL remains a challenge.
Both of the micromechanical models and cohesive-zone model are based on the consideration
of local failure mechanisms within a small area around the crack-tip (i.e. fracture process zone).
On the other hand, the linear elastic fracture mechanics (LEFM), which analyzes sharp-tip macro
cracks, has also been applied to fiber-reinforced composite laminates. The LEFM treats the fiberreinforced composite laminates as homogeneous materials, instead of considering heterogeneous
micro-structures. Thus, it has higher conceptual simplicity and computational efficiency
compared to the cohesive zone models and micromechanical models. However, the limit of
LEFM is also significant. In the classical LEFM, the fracture toughness of cracked body is
measured by a single parameter--the critical stress intensity factor, which is assumed to be a
material property that is independent of structural size and geometry. This single-parameter
fracture toughness is valid for most brittle materials, in which the size of crack-tip fracture
process zone is negligible compared to the characteristic structural dimensions (e.g. the original
crack length). However, usually the crack-tip fracture process zone in fiber-reinforced composite
laminates and other quasi-brittle materials is relatively large. Consequently, the fracture
toughness of cracked structures made of fiber-reinforced composite laminates cannot be
represented by a single size-independent constant. Instead, the crack growth resistance curve, i.e.

4

R-curve, is considered as an appropriate alternative to represent the fracture toughness of
materials. Moreover, when the crack-tip fracture process zone becomes too large to satisfy the
“small-scale yielding” condition, even the R-curve depends on the structure geometry and size.
For this “large-scale yielding” case, the LEFM cannot be used for the fracture analysis any more
[15]. The interlaminar fracture of some fiber-reinforce composite laminates falls into this
category of fracture problems.
Besides of the above fracture models, it is worth mentioning that a variety of fracture models
were developed as modifications to the LEFM by compensating the size effect caused by the
relatively large crack-tip fracture process zone in fiber-reinforced composite laminates. These
fracture models were developed in terms of empirical equations mainly for engineering purposes.
For instance, Waddoups et al [16] employed an empirical extension of LEFM to predict the
notched strength of composite panels. They argued that the damage zone adjacent to the cracktip effectively increases the crack length, similar to Irwin’s plastic zone correction in metals. The
fracture then occurs when the applied stress intensity factor determined by the load and the
modified crack length reaches a critical value. This fracture model is denoted as the Inherent
Flaw Model. Whitney and Nuismer [17], [18] proposed two stress-based criteria which are
similar to the Inherent Flaw Model, namely the Point Stress Criteria and Average Stress Criteria
respectively. The Point Stress Criteria assumes that fracture occurs when the stress at a
characteristic distance ahead of the notch tip reaches the un-notched tensile strength, whereas the
Average Stress Criteria assumes that fracture occurs when the average stress over a characteristic
distance ahead of the notch tip equals the un-notched tensile strength. Subsequent modifications
to the Whitney and Nuismer model include the work of Karlak [19] and Pipes et al [20]–[22],
which introduced more fitting parameters into the fracture criteria in addition to the un-notched

5

tensile strength and the characteristic distance. An comprehensive review of these models has
been published [23]. Even though these models worked well for specific materials, a basic
problem is that the characteristic distance or other parameters are empirically determined, so that
the universal applicability of these models to other materials or other cracked body
configurations with the same material cannot be guaranteed.

1.2 Objectives and scopes
As stated above, the heterogeneous microstructures in fiber-reinforced composite laminates
introduce various ply-level failure modes into their fracture processes, including interlaminar
fracture, intralaminar fracture and translaminar fracture. While a great number of published
studies have been dedicated to the investigation of interlaminar fracture of fiber-reinforced
composite laminates, their translaminar fracture mechanisms have received relatively little
attention from the fracture mechanics community. However, with the increasing use of large
composite primary structures in industry, it is envisaged that a comprehensive understanding of
the translaminar fracture of fiber-reinforced composite laminates will play an increasingly
important role over the coming years [24].
In this thesis, the mode I translaminar fracture behavior of a cross-ply glass fiber reinforced
composite laminates is investigated, in the framework of linear elastic fracture mechanics
(LEFM) and cohesive-zone model (CZM). The objective of this study is to investigate the
applicability of these two fracture models for the translaminar fracture analysis of any fibrereinforced composite laminates which can be approximated as anisotropic solids and show selfsimilar crack growth. Specifically, two sets of analytical characterizations of the elastic fields

6

near crack-tip in anisotropic solids are separately derived within the frameworks of LEFM and
CZM. Then two experimental methods are separately developed to measure the R-curve of
LEFM and the cohesive-zone law of CZM, based on the analytical characterization of crack-tip
fields and full-field crack-tip displacement measured with Digital Image Correlation (DIC)
technique. Eventually, the applicability of LEFM and CZM for the translaminar fracture analysis
of fibre-reinforced composite laminates are investigated.

1.3 Structure of dissertation
The remaining part of dissertation is organized as follows. In Chapter 2, the literature review
regarding the up-to-date progress of fracture mechanics studies with LEFM and CZM for fiberreinforced composite laminates is provided. Besides, the full-field optical measurement
technique in the study-Digital Image Correlation (DIC) is briefly introduced.
In Chapter 3, the translaminar fracture problem is studied in the context of LEFM. In this study,
the fiber-reinforced composite laminates are approximated as homogeneous anisotropic solids.
Firstly, an asymptotic expansion of crack-tip field in anisotropic solids is derived with the
anisotropic plane elasticity theory. Based on the asymptotic expansion, two methods are
developed to estimate important crack-tip field parameters from crack-tip displacement fields.
One method employs the nonlinear least-squares fitting algorithm, while the other method makes
use of the conservation integrals of fracture mechanics. The two methods respectively
correspond to the two major analysis approaches in LEFM, i.e. the stress intensity approach and
the energy approach. The investigated crack-tip field parameters include stress intensity factor,
energy release rate and effective crack length. An experimental validation of mode I translaminar
fracture test with a cross-ply glass fiber-reinforce composite laminates and the extended compact
7

tension (ECT) specimen geometry is also conducted. The crack-tip field parameters and Rcurves are obtained from full-field crack-tip displacement fields measured by 2-D Digital Image
Correlation. Special considerations are made to investigate the errors and limitations of the
methods. Then the R-curve behaviors for the translaminar fracture of fiber-reinforced composite
laminates are discussed.
In Chapter 4, the translaminar fracture problem of fiber-reinforced composite laminates is
investigated in the framework of cohesive-zone model (CZM). The homogeneous anisotropic
approximation of fiber-reinforced composite laminates retains in this chapter. Firstly, an
eigenfunction expansion of the planar elastic field around the crack-tip cohesive zone is derived
for anisotropic solids. Two analytical inverse methods, namely field projection method and leastsquares method, are developed based on the eigenfunction expansion to extract the details of
cohesive-zone law (CZL) from the anisotropic elastic fields around the crack-tip cohesive zone.
A sequence of numerical tests are carried out to assess and compare the accuracy and stability of
the two methods, through the extraction of cohesive-zone laws from synthetic displacement
fields that are generated with predefined analytical solutions. In the numerical tests, several
factors that can influence the accuracy of inverse solutions are investigated, including the shape
of cohesive-zone law, the noise level of inputs, the inverse distance of the data field and the
number of data points. Then some guidelines for the implementation of inverse methods in
practice are provided based on the results of numerical tests. In Chapter 5, the inverse methods
are applied to the same translaminar fracture test for the study in the context of LEFM. To verify
the accuracy of extracted CZLs, finite element simulation results which use the extracted CZLs
are compared with the experimental outputs. Finally, some concluding remarks and expectations
for the future work are provided in Chapter 6.

8

Chapter 2
Literature Review
In this chapter, a comprehensive literature review regarding to the backgrounds of this
dissertation is provided. The whole chapter is divided in to three sections. In Section 2.1 and
Section 2.2, the research advance of linear elastic fracture mechanics (LEFM) and cohesive-zone
model (CZM) in the fracture analysis of fibre-reinforced composite laminates are introduced,
respectively.To keep consistent with the research objective of this thesis, the emphasis of this
literature review is laid on the translaminar fracture of fibre-reinforced composite laminates.At
the end of Section 2.1 and Section 2.2, the advantages and disadvantages of LEFM and CZM are
respectively summarized to illustrate the motivations of the study. In Section 2.3, the principle
and current development of the optical measurement technique used in this study—Digital Image
Correlation (DIC) are briefly introduced.

2.1 Application of linear elastic fracture mechanics in fiberreinforced composite laminates
Historically, the great success of LEFM in analyzing the fracture behavior of homogeneous
isotropic materials has motivated researchers to explore its applicability in the fracture analysis
of fiber-reinforced composite laminates. It is well known that the macroscopic elastic properties
of fiber-reinforced composite laminates can be considered as anisotropic. Thus, the application
of LEFM in fiber-reinforced composite laminates requires additional considerations of the

9

material anisotropy. In 1965, Sih et al. [11] first demonstrated the feasibility of extending LEFM
analysis methods for isotropic materials to anisotropic solids. In their study, a complex variable
approach was used to obtain a mathematical understanding of the stress fields and displacement
fields around a sharp-crack tip in anisotropic elastic solids. Mathematical expressions of the
crack-tip fields were derived in the form of asymptotic expansions, which is similar to the
Williams expansions for isotropic materials. More importantly, they concluded that the concepts
of stress intensity factor and fracture toughness can be extended into cracked anisotropic solids,
whereas the relationship between the energy release rate and stress intensity factors should be
modified to involve the effect of anisotropic material property.
Based on the pioneering work of Sih et al, many attempts were made to apply the LEFM
analytical solutions or formulas which were originally derived for isotropic materials to study the
fracture problems of fiber-reinforced composite laminates. In these studies, the fiber-reinforced
composite laminates are approximated to be anisotropic materials. Many researchers have
studied the influence of the anisotropic material properties on the accuracy of stress intensity
factors estimated with the formulas of isotropic materials. Mandell etal. [25] used hybrid finite
element analysis to find that the degree of material anisotropy has a varying effect on the
calculation of stress intensity factors for different specimen configurations, whereas this effect of
anisotropy is relatively constant for varying crack lengths in a given configuration. As an
advanced study, Sweeney [26], [27] calculated the finite-width correction factors for singleedge-notched tension test of orthotropic materials, through a finite element analysis with fracture
mechanics J-integral. Later on, Bao et al. [12], [13] developed a rescaling technique to quantify
the effect of material orthotropy on the stress intensity factor formulas of notched bars,
delaminated beams, and hybrid sandwiches. All the studies demonstrated that the elastic

10

anisotropy may introduce a considerable effect into the measurement of stress intensity factor
and fracture toughness of fiber-reinforced composite laminates with analytical stress intensity
factor formulas.
The material anisotropy is not the only concern in the application of LEFM for fiberreinforced composite laminates. Many direct applications of LEFM to this type of materials have
suffered as a consequence of violating the “small-scale yielding” condition. As a prerequisite of
the validity of LEFM theories, the “small-scale yielding” condition requires that the nonlinear
deformation zone (which is also referred to as fracture process zone in some publications)
surrounding the crack tip must be small compared to the other characteristic lengths (e.g. crack
length) in the cracked geometry. As the size of crack-tip nonlinear deformation zone grows, the
fracture analysis with LEFM becomes increasingly inaccurate. Indeed, the fracture behaviors of
materials can be categorized based on the size of crack-tip nonlinear deformation zone.
Appropriate fracture models should be chosen accordingly. Firstly, when the crack-tip nonlinear
deformation zone is negligible compared to the crack length (this is the case for most brittle
materials such as glass), the instant of fracture can be predicted by a single-parameter fracture
toughness, i.e. the critical stress intensity factor Kc or the critical energy release rate Gc. The Kc
and Gc for brittle fracture are independent on the structure configurations and hence are material
constants. Secondly, for the fracture of quasi-brittle materials such as fiber-reinforced composite
laminates, the crack-tip nonlinear deformation is moderate. In this case, the stress intensity factor
and energy release rate keep increasing at the starting stage of crack extension, which is caused
by the increase of crack-tip nonlinear deformation zone. Consequently, a single-parameter
fracture toughness such as Kc or Gc is no longer suitable for fracture prediction. Nevertheless, the
validity of LEFM can be retained by introducing an alternative fracture resistance measure—the

11

crack resistance curve (R-curve), in which the material fracture resistance is plotted with respect
to the crack extension. The material fracture resistance corresponding to a certain crack
extension is usually equivalent to the stress intensity factor or energy release rate which can
cause stable crack growth at the same instant. For the quasi-brittle fracture problems, it has been
found that the R-curve is independent of structure configurations and hence is a material property
[15]. Moreover, for more extensive yielding cases which totally violate the “small-scale yielding”
condition (fracture of ductile materials such as steel), the R-curve becomes dependent on the
shape and dimension of cracked solids. In this case, ductile fracture models must be applied to
take into account the considerable nonlinear material behavior. Some two-parameter fracture
theories were developed to provide size-dependent fracture criteria, such as the J-Q theory [28],
[29] and the cleavage scaling model [30], [31]. The cohesive-zone model, which will be
introduced in the next section, can also be used to analyze the ductile fracture problems. The
cohesive-zone law has been proved to be a size-independent fracture criterion even under largescale yielding condition [15].
Besides the development of theoretical models to describe the complicated fracture
mechanisms of fiber-reinforced composite laminates, research efforts to understand such
complex fracture processes also heavily rely on experimental observations and measurements.
Up until recently, most of the experimental measurements of the single-parameter fracture
toughness are based on measuring far-field loads and overall deformations during the mechanical
test. Then, the crack-tip field parameters, i.e. the stress intensity factor or strain energy release
rate, are inferred through analytical formula or numerical calibrations. The fracture toughness, i.e.
the critical stress intensity factor or critical energy release rate, is estimated with the maximum
far-field load and original crack length. However, this method has two significant drawbacks.

12

Firstly, as stated above, the single-parameter fracture toughness has limited applicability for the
fracture analysis of fiber-reinforced composite laminates because the stress intensity factor or the
energy release rate keeps varying at the starting stage of crack extension. Secondly, even though
the "critical" moment is defined as that with the maximum far-field load, the original crack
should not be used in the estimation of fracture toughness. The crack can increase for a
considerable length through the growth of crack-tip nonlinear deformation zone before reaching
the instant with maximum far-field load.
Instead of the single-parameter fracture toughness, the translaminar fracture behavior of fiberreinforced composite laminates was found to be best described by R-curves [24]. However, the
experimental measurement methods of R-curve for translaminar fracture of fiber-reinforced
composite laminates have not been standardized yet even for mode I tensile failure. Recently, a
extensive review of the state-of-the-art experimental measurement methods of translaminar
fracture toughness and R-curve was provided by Laffan et al. [24], [32] . They investigated the
commonly-used specimen configurations used for the fracture toughness and R-curve
measurements, including compact tension, three or four point bending, double edge notched
tension, extended compact tension, center notched tension and single edge notched tension.
Among these specimen configurations, only compact tension, three/four point bend and extended
compact tension specimens could exhibit stable crack growth and are applicable for measuring
the R-curves. In addition, several data reduction methods for the R-curve measurements were
compared. The methods under investigation included the area method, compliance calibration
method, modified compliance calibration method, and finite element analysis using the J-integral
computation or virtual crack closure technique. Among these methods, the “Modified
Compliance Calibration method” (MCC) was found to provide the most reliable measurements.

13

Because this method did not rely on the optical measurement of the crack length in which the
measurement accuracy is highly dependent on the experience and cautiousness of the
experimenter. However, the MCC method still requires much experimental work by measuring
the elastic compliances for various specimens with machined cracks of different lengths [32] .
On the other hand, owing to the advances in full-field optical measurement techniques, the
crack-tip field parameters can be directly determined from the deformation fields around the
crack-tip which can be measured by the optical techniques. The applications of optical
techniques to characterize the fracture toughness of engineering materials have already become
routine tasks at least for isotropic materials [34]–[38]. In one commonly used data reduction
scheme, the crack-tip field parameters such as the stress intensity factors and effective crack
length were inversely extracted by fitting the Williams’ expansion of isotropic crack-tip elastic
fields [39] to the measurement data [34], [37] in a least-squares sense. Some other methods were
developed based on the conservational integrals, such as the J-integral [40] and interaction
integrals [33]. The conservational integrals were calculated from the deformation data around the
crack-tip. The energy release rate can be directly obtained by calculating the J-integral, and the
mode I and mode II stress intensity factors can be determined by using the interactional integral
[36].
In contrast to the fruitful application of optical techniques for the measurement of fracture
resistances of isotropic materials, a relatively fewer number of investigations based on optical
techniques were reported for fiber-reinforced composite laminates [42]–[45]. The major reason is
that the extraction and interpretation of the crack-tip field parameters in composite laminates
require special considerations on the material heterogeneity and elastic anisotropy. Among the
published literature, most of the studies only measured the stress intensity factors, which implied

14

that the critical stress intensity factor is used as the fracture toughness. As stated above, this
single-parameter fracture toughness is not suitable for the fracture prediction of most fiberreinforced composite laminates. Recently, an attempt has been made to measure the R-curve of
translaminar fracture of composite laminates based on the energy release rate consideration. The
energy release rate was obtained by the calculation of J-integral [44]. However, the accuracy of
the measured R-curve in this work is somewhat questionable because the crack-tip position was
identified by searching the displacement discontinuity in the crack-tip displacement fields.
Generally, the optical measurement of the displacement fields very close to the crack-tip is
unreliable due to the high displacement gradients within this area.
Moreover, a significant drawback in most of the above studies based on the crack-tip full-field
optical measurement is that the small-scale yielding condition, which is the prerequisite for the
validity of linear elastic fracture mechanics, were not validated. Therefore, there was lack of
sufficient evidence that the applicability of the measured single fracture toughness or R-curve in
these works was a geometry-independent material property.
To sum up, currently there are two major challenges in the application of LEFM to fiberreinforced composite laminates. The first challenge is the accurate measure of R-curve without
much experimental complexity. The second challenge is the verification of the geometryindependence of R-curve, which is, equivalently, the verification of small-scale yielding so that
this LEFM fracture resistances is still valid.

15

2.2 Cohesive-zone model
Recently, a nonlinear fracture model named cohesive-zone model (CZM) has become more
and more popular in the analysis of quasi-brittle and ductile fracture behaviors. Particularly, it
has been widely used to study the nonlinear fracture behaviors of fiber-reinforced composite
laminates. In the CZM, a cohesive zone is used to simplify the nonlinear fracture process zone
around the crack-tip. Instead of considering the complicated micron-mechanical fracture
mechanisms within the crack-tip fracture process zone, the CZM characterizes the crack-tip
fracture process and the energy dissipation during crack propagation with the relationship
between the cohesive tractions and separations. This relationship constitutes the cohesive-zone
law (CZL), a material-dependent constitutive law. The illustrations of the CZM and cohesivezone law are shown in Figure 2.1 and Figure 2.2, respectively.

Figure 2.1 Cohesive-zone model

16

Figure 2.2 Cohesive-zone law

The concept of CZM was first introduced by Barenblatt [8]. He originally postulated the
existence of cohesive forces in the vicinity of crack tip to eliminate the unphysical stress
singularity integrated in the linear elastic fracture mechanics (LEFM). Next, Dugdale [9] used a
CZM to represent the fracture in perfectly plastic materials. In his model, the cohesive tractions
within the crack-tip plastic zone were assumed to be equal to the yielding strength. Hillerborg et
al. [10] also analyzed the crack growth in quasi-brittle concrete using the CZM (the model was
named as “fictitious crack model”). Since then the CZM has gained great success in analyzing a
variety of nonlinear fracture processes, such as the interface decohesion between dissimilar
materials [46], the void-growth in ductile metals [47], the crazing in glassy polymers
(Tijssens et al. 2000), the void nucleation at an interface between a particle and matrix [49], and
the fracture mechanisms of quasi-brittle materials like fiber-reinforced composite laminates [50].
Compared to other fracture models, the CZM has significant advantages, especially in the
analysis of nonlinear fracture problems. Firstly, unlike the fracture resistance measures in the
17

context of LEFM (such as critical stress intensity factor and R-curves), the CZL is an inherent
material property that is valid even under the large-scale yielding condition. Thus the CZM can
be applied to not only brittle and quasi-brittle, but also ductile fracture problems. Secondly, the
CZM is particularly well-suited for the finite element simulation of nonlinear fracture behaviors,
since it not only predicts the fracture initiation but also the crack propagation process. The
cohesive-zone modeling algorithms have already been integrated into some of the commercial
FEA software such as ABAQUS and ANSYS.
For the finite element simulation with CZM, it has been found that some global fracture
behaviours of solid structures can be approximately predicted by the CZLs with simple
functional shapes. For instance, the linear CZL in which the cohesive strength and fracture
energy are the sole important factors to define its shape, were widely used in simulations.
However, some recent investigations have shown that the functional shape of CZL could
influence the simulation results considerably. Gu [51] showed that, under the large-scale
bridging condition, the shape of the CZL influenced the load carrying capacity of a plate with an
elliptical hole. Suo et al. [52] and Rice et al. [53] found that the stability of elastic solids with
interfaces, including not only the bifurcation of de-cohesion processes but also frictional slip
instability, was sensitive to the details of CZLs. Chandra et al. [54] investigated the fracture
processes in diverse material systems by utilizing several CZLs and showed that the shape of the
CZL played a critical role in determining the macroscopic response of a composite system.
Volokh [55] compared the bilinear, parabolic, sinusoidal and exponential CZLs using a blockpeel test. They showed that the shape of CZL significantly affected the result. Alfano [56]
reported the differences caused by the use of different CZLs such as bilinear, linear-parabolic,
exponential and trapezoidal laws. They concluded that the effect of the shape of CZL on

18

numerical results depended on the boundary conditions of the problem and the ratio between
interface toughness and stiffness of the bulk material. Song et al. [57] also studied the influence
of the shape of CZLs on the fracture modeling of asphalt concrete. A nonlinearly decaying CZL
was found to be more appropriate for simulation than a linear CZL. In general, the influence of
the shape of CZL became considerable when the size of crack-tip fracture process zone was
relatively large compared to other structure dimensions, which was relevant to most of the quasibrittle and ductile materials.
Therefore, an accurate estimate of CZL is essential to achieve truly predictive simulations for
quasi-brittle materials such as fiber-reinforced composite materials. Several experimental
approaches have already been developed to extract cohesive-zone properties. Among them, the
uniaxial tension test is considered to be the most direct method to determine the mode I CZL. In
this experimental method, the curve which represents the relationship between the cohesive
traction and separation, i.e. the CZL, was directly measured [58]–[60]. However, some
theoretical and experimental issues prevented the wide application of this method. Firstly, the
measured cohesive traction was approximated with the far-field tensile stress, which was against
the basis of fracture mechanics that the tensile stress is highly concentrated around the crack-tip.
Furthermore, the CZL was assumed to be of simple functional shape, either because the number
of measurable outputs was limited [58] or because the significant noise made the details of
measurement unreliable so that the measured CZL had to be fitted to a simple functional shape
[60]. In addition, there were also a lot of difficulties in the setup of experiments [59].
There are several other “direct methods” which try to directly measure the characteristic
parameters that describe the shape of CZL. In these methods, the shape of CZL is usually
presumed to be linear or bilinear, which can be characterized with a few characteristic

19

parameters such as the tensile strength, fracture energy and critical crack opening displacement.
Then the characteristic parameters can be determined through experiments. For instance, Park et
al. [61], [62] assumed that the kink point in a bilinear CZL corresponded to the instant when the
critical crack-tip opening displacement was reached. Then the bilinear CZL could be determined
based on a single fracture test. The defects of these “direct methods” are obvious. Firstly, the
estimated CZL is dependent on the assumption of its shape. Secondly, these methods cannot
obtain the details of CZL.
It is worth mentioning that there is another method which utilizes the mathematical fact that
the relationship between the cohesive-zone traction and separation, i.e. the CZL, is a
differentiation of the energy release rate (J) with respect to the crack-tip opening displacement
(CTOD) [63]–[65]. The theoretical design of this method is valuable, but the application of this
method suffers from the difficulty in the experiment setups. The noise in the measurement of
crack-tip opening displacement is significant, since the crack-tip cohesive zone is small and
difficult to recognize. Moreover, the numerical error of differentiation is considerable because of
the noisy measurement data of J and CTOD. Thus the measured CZLs in these studies were
scattered. In addition, this method is only applicable for mode I problems.
Besides the “direct methods”, inverse methods have also been developed based on the
parametric fitting of experimental data with optimization algorithms. The common procedure of
these methods was summarized by Elices et al. [66]. Firstly, a functional form of CZL with a few
characteristic parameters was assumed, depending on the material system and corresponding
fracture processes. Then the characteristic parameters in the functional form were computed by
iterative optimization procedure, i.e. the parameters were adjusted at each iteration by comparing
the computational and experimental results. The most important advantage of these inverse

20

methods is that the details of CZL are obtainable by assuming flexible functional forms.
However, these methods are also accompanied with defects. In most of the studies, only the
global responses of the experimental data, i.e., load-displacement or load-crack mouth opening
displacement (CMOD), were used as the basis of the parametric fitting procedure. This was due
to the fact that the global responses were the only obtainable experimental outputs in most of the
studies. The usage of global response would exacerbate the non-uniqueness problem of the
inverse solution, i.e. the CZL, which is a highly localized material property.
Recently, the application of full-field optical techniques like Digital Image Correlation (DIC)
in the studies of material fracture behaviors has become more and more popular. The full-field
optical technique can provide rich measurement data within a small area surrounding the crack
tip. Thus it is suitable for the investigation of the localized material fracture properties. In
Section 2.1, the application of full-field optical techniques within the framework of linear elastic
fracture mechanics has been reviewed. The application of these techniques in the estimate of
CZL is also fruitful.
The methods of measuring CZLs with the full-field optical techniques can be categorized as
direct methods and inverse methods, too. One representative direct method is to directly measure
the cohesive tractions from the deformation fields of bulk material surrounding the cohesive
zone. The method assumed that the cohesive tractions were in equilibrium with the surrounding
stress fields in the bulk material which could be obtained with the stress-strain relationship and
bulk material properties. Then the CZL could be determined by correlating the directly measured
crack opening displacement with the cohesive traction [67], [68]. Obviously, the accuracy of the
measured CZL is dependent on the definition of cohesive zone within the deformation fields,
which is highly subjective. Another work used part of the whole-field displacement field, i.e. the

21

displacement profiles along edges of the specimen, to solve for the characteristic parameters in a
simple functional form of CZL [69]. However, the result was found to be sensitive to the location
of the measurement and the measurement errors.
As for the inverse method, Guo et al. [70] determined a linear CZL by iteratively fitting the
finite element simulation result of crack opening displacement to the experimental result. As a
further step, a hybrid inverse technique combining finite element simulation and full-field
measurement has been developed to determine CZLs with flexible functional forms [71]–[73] .
First, a flexible shape parameterization of the CZL was used without making any assumption
about its shape. Then the parameters used to characterize the shape of CZL were determined
through the unconstrained, derivative-free Nelder-Mead optimization method. The objective
function used in the optimization procedure was the norm of the difference between the
measured whole-field displacement data and the computed whole-field displacement field from
FEM. Several successful applications of this methods have already been reported [74]–[76].
However, this method still has drawbacks. The first one is its high computation cost, since a
complete numerical simulation must be carried out at each iteration in the optimization
procedure. Furthermore, the accuracy of the extracted CZL is reduced by using the displacement
data within the crack-tip high-displacement-gradient field in the finite element computation.
Recently, an analytical inverse method named the “field projection method” was developed by
Hong and Kim [75]. With this method, the CZLs can be inversely extracted from the elastic farfields from the crack tip. Therefore, it can reduce the error caused by the usage of unreliable
displacement data very close to the crack-tip. In their original work, a general form of cohesivecrack-tip fields was first expressed in terms of an eigenfunction expansion of the plane elastic
field using a complex variable representation. The coefficients of the eigenfunction expansion

22

within the cohesive zone were extracted by applying interaction J-integrals between the auxiliary
probing fields and the elastic far-fields from the crack-tip. Then the cohesive-zone variables, i.e.
the cohesive-zone tractions and separations, and thus the CZLs were separately determined by
using the coefficients of the eigenfunction expansion. An iterative method to determine the
position and size of a cohesive zone was also proposed in the paper. The applicability of the field
projection method for isotropic materials has been verified experimentally and numerically
[76], [77]. Furthermore, the field projection method has been extended to extract the CZLs of an
interface between two anisotropic solids at the nano scale [78]. The theoretical derivation of this
field projection method enhances the mathematical understanding of the CZM. Since this method
is developed based on the analytical representation of cohesive-crack-tip fields, it provides more
confidence to the accuracy of inverse solutions than the method with flexible shape
parameterization of CZL [74]–[76]. Moreover, the CZLs determined by the analytical methods
have nearly no restriction on their functional shapes, because the order of polynomial that is used
to describe the functional shape of CZL is flexible. Therefore, this analytical method has the
potential to overcome the defects of numerical inverse methods. However, the considerable illconditioning [75] and the high sensitivity to experimental noise [78], [79] are the two major
defects of this method.
In sum, the major challenge in the application of CZM for the fracture analysis of fiberreinforced composite laminates is the accurate measurement of the CZLs.

23

2.3 Digital Image Correlation (DIC)
In the present study, the crack-tip displacement fields on the surface of cracked specimen is
measured by Digital image correlation (DIC), a non-contacting full-field optical technique that is
originally developed in the 1980s [79]–[81]. Through 30 years' development, now DIC is one of
the most commonly-used full-field optical techniques, which is capable of measuring not only
the two-dimensional but also the three-dimensional surface deformations [82]. Countless
applications with DIC have been obtained in a vast variety of academic and industrial fields. In
the present study, the two-dimensional DIC is used for the measurement of displacements. A
brief introduction of this technique is provided in this section.
Compared to other optical measurement techniques such as the interferometry-based optical
techniques (e.g. Electric Speckle Pattern Interferometry (ESPI)), DIC has several special
advantages. First of all, the experimental setup and specimen preparation of DIC are simple. The
core of the experimental setup is an image acquisition system, which includes a light source for
illumination, a CCD (charge-coupled device) camera and an image storage device. In addition,
the only requirement of specimen is that the specimen surface in captured images must have a
random gray intensity distribution. It is easy to fulfill this requirement either by directly using the
natural texture of the specimen surface or by spraying black and white paints onto the surface to
generate speckle patterns. A typical image with speckle patterns is shown in Figure 2.3.
Secondly, the measurement environment of DIC has low limitation so that it is appropriate for
not only laboratory but also field applications. Here the "measurement environment" mainly
refers to the illumination conditions while capturing images with CCD camera. The laser source
is not required for DIC [83]. A while light source or natural light is good enough. Thirdly, the
range of measurement resolution of DIC is adjustable depending on the resolution of digital

24

image. For instance, with high-spatial-resolution image acquisition equipments such as optical
microscope, scanning electron microscope (SEM) and atomic force microscope (AFM), microscale and even nano-scale deformation can be measured with DIC [84]–[91].

Figure 2.3 A typical image used in DIC with random gray intensity distribution. The speckle
patterns are generated by spraying black and white paints on the specimen surface.

Essentially, DIC measures the surface displacement fields of specimen by tracking the
movements of points in the two images captured before and after deformation, which are referred
to as the reference image and the deformed image, respectively. The measurement procedure
consists of three successive steps. Firstly, the region of interest in the reference image is divided
into evenly spaced virtual grid. Then the displacement of each node in the grid is measured
separately. Eventually the full-field surface deformation can be obtained by measuring the
displacements for all the nodes.
The basic principles of DIC in measuring the displacements of each node are schematically
illustrated in Figure 2.4. To measure the displacement of node P in the reference image, a square
25

subset centered at P is used to determine the corresponding position of P in the deformed image,
i.e. the position of P'. A subset is used because it has a unique gray intensity distribution which
can distinguish itself from other subsets, so that it can be identified in the reference and
deformed images.

Figure 2.4 A schematic illustration of the basic principles of DIC

To identify the subset in the reference and deformed images, a criterion is needed to evaluate
the degree of similarity between two subsets. In DIC, the similarity between two subsets is
measured by a parameter named "correlation coefficient", which compares the gray level
distributions within two subsets. In other words, a combination of the grey values of all the
pixels in a subset is used as the comparable parameter in the correlation coefficient. There are
plenty of correlation coefficients in the literature, which can be categorized into two groups: the
26

cross correlation coefficients and the sum of squared differences coefficients [92], [93].
Although the formats of the coefficients within two groups seem different, it has been
mathematically proved that there are one-to-one correspondences between the coefficients within
two groups. Moreover, the zero-normalized cross-correlation coefficient or the zero-normalized
sum of squared differences coefficient is believed to be the most robust coefficient because of
their capabilities to offset the influence of linear variation of lighting [94].
Another important feature that should be considered in the identification of the deformed
subset is the shape change of subset during deformation. As shown in Figure 2.4, the shape of
subset in the reference image is a square, while the shape of corresponding subset in the
deformed image is distorted. This shape change can be represented mathematically. If the
coordinates of node

is ( , ), the coordinates of a random point

within the subset can be

represented as
.
After deformation, the point

moves to

(2.1)

, whose coordinates can be represented as
,

where and

(2.2)

are named shape functions [95] or displacement mapping functions [96].

Actually, the shape functions are derived based on the Taylor-series expansion of the
displacement fields in the vicinity of node . The form of shape functions is determined by the
type of shape change of subset. Firstly, if there is no shape change between the subsets before
and after deformation, i.e. there is only pure rigid body translation, the zero-order shape
functions are sufficient to describe the movement of square subset:
.
27

(2.3)

However, when there is a shape change in the deformed subset, the zero-order shape functions
are insufficient. In this case, the first-order shape functions can be used:

(2.4)

It can be seen that the translation, rotation, normal and shear strains of subsets are involved in the
first-order shape functions. Besides, the second-order shape functions which can represent more
complicated deformation situations are also available in the literature [96]:

.

(2.5)

However, advanced studies already showed that the first-order shape functions was a good
approximation to the shape change of subset in most cases. The application of second-order
shape functions is unnecessary, considering the much higher computational complexity and
limited improvement in the measurement accuracy [95].
Since the shape of deformed subset is changed and the point

may be located between

integer pixels, the grey values at the sub-pixel locations of deformed image are required to
compute the correlation coefficient. The sub-pixel grey values can be obtained by interpolation
methods, including the bilinear interpolation [81] and other higher-order interpolations such as
the bicubic spline interpolation [97] and biquintic spline interpolation [98]. The higher-order
interpolations can provide higher measurement accuracy compared to the bilinear interpolation.
The core of DIC implementation procedure is to search the position of the extreme of
correlation coefficient. In details, the procedure to determine the displacements of one node
includes two consecutive steps. Firstly, the integer-pixel-level displacements of the node are
obtained. Next, the sub-pixel-level displacements are computed starting from the integer-pixel28

level displacements. In the first step, usually the shape of subset is assumed to retain during
deformation, since only a rough estimate of displacement is needed. The integer-pixel-level
displacements can be obtained by simply applying a pixel-by-pixel searching scheme or the
coarse-to-fine searching scheme [99]. An alternative method is to calculate the correlation
between two subsets in Fourier's domain [100]. The computational speed of this method is
extremely fast, in spite of the fact that its measurement accuracy is sensitive to the displacement
gradients. In the second step, the measurement accuracy of displacements is enhanced to subpixel level, using the integer-pixel-level displacements as the initial guess. In the literature, there
are plenty of sub-pixel registration algorithms which can achieve sub-pixel measurement
accuracy as small as 0.01 pixels, such as the coarse-fine search method [101], curve fitting of the
correlation coefficient method [102], and Newton-Raphson iteration method [81], [97]. Among
all the algorithms, the Newton-Raphson iteration method is the most widely used one due to its
balanced measurement accuracy and computational cost.
In this thesis, a 2D-DIC software was custom-built for the measurement of crack-tip
displacements. As for the algorithms used in the DIC software, the zero-normalized sum of
squared differences coefficient and the first-order shape function shown in (2.4) are chosen; the
sub-pixel grey distribution of the digital image is obtained by the bicubic spline interpolation; the
measurement of integer-pixel-level displacements is implemented in Fourier's domain; the subpixel measurement is computed iteratively with the Levenberg-Marquardt algorithm [95]. The
measurement accuracy of the software is estimated to be 0.02 pixels, through a sequence of
numerical tests with synthetic speckle patterns. The DIC software was coded in MATLAB. The
source code is provided in Appendix A.

29

Chapter 3
R-curve Characterization of Translaminar Crack Growth in
Cross-ply

Composite

Laminates using

Digital

Image

Correlation
3.1 Introduction
This chapter presents experimental methods to estimate crack-tip field parameters and
characterize R-curve behaviors of translaminar fracture in cross-ply composite laminates, using
the experimental outputs of Digital Image Correlation (DIC). The investigated crack-tip field
parameters include stress intensity factor, energy release rate and effective crack length. The
nonlinear least-squares method and conservation integrals are used to determine these parameters
from elastic far-field displacements, based on the homogeneous approximation of anisotropic
solids for fiber-reinforced composite laminates. The sources of error of the two parameterestimation methods are investigated. Eventually, R-curves are constructed using the estimated
crack-tip field parameters.
The remainder of this chapter is organized as follows. In Section 3.2, firstly, the description of
asymptotic crack-tip fields in general anisotropic elastic solids is presented following the Stroh's
representation of anisotropic elasticity [103]–[106]. Then two methods, which are respectively
based on the stress intensity factor method and energy release rate consideration, are applied to
investigate this quasi-brittle fracture problem. A nonlinear least-squares method is proposed to

30

estimate the stress intensity factor and effective crack length from crack-tip deformation fields,
similar to the method proposed by Yoneyama et al [37]. In parallel, conservation integrals [40],
[107], [108] are also introduced to estimate the energy release rate and effective crack length. In
Section 3.3, both the nonlinear least-squares method and conservation integrals are applied to
process the experimental crack-tip displacement fields measured by Digital Image Correlation.
Special considerations are made to investigate the errors and limitations of two methods. Then
the results of estimated crack-tip field parameters and the R-curve behaviors are analyzed. In
Section 3.4, more extensive applications of these parameter-estimation methods are discussed.
Finally, some conclusion remarks are provided in Section 3.5.

3.2 Estimation of the crack-tip field parameters in anisotropic solids
3.2.1 Asymptotic expansion of crack-tip fields in anisotropic solids
In this section, a general form of crack-tip elastic fields in anisotropic solids is expressed in
terms of an asymptotic expansion of the complex functions in the Stroh formalism [104], [106].
In order to describe the present study in a self-contained manner, we firstly describe
anisotropic elasticity theory that has been well-documented in the work of Stroh and Suo [104],
[106]. To avoid confusion, no summation convention is assumed for repeated indices.
The generalized Hooke’s law relating stresses  ij to strains  kl for an anisotropic material can
be written as
ij 

3

 Cijkl kl

k , l 1

31

(3.1)

where Cijkl is the elastic stiffness tensor. Using the symmetry of the stiffness tensor and the
definition of infinitesimal strain tensor, the equilibrium equation can be written as
3



j , k , l 1

C ijkl

 2u k
 0.
xl x j

(3.2)

A general solution for the displacement u k in the two-dimensional plane ( x, y) is given
[103], [105] as
3

uk  2 Re  Ak f  ( z ) ,
 1


(3.3)

where f  ( z  ) are analytic functions of the complex arguments, z   x  p  y , and Ak is a 3 3
matrix depending on elastic constants. Each column of the matrix Ak and each of the
characteristic roots p are the eigenvector and the corresponding eigenvalue with positive
imaginary part of the following eigenvalue problem, respectively:
3

[Ci1k1  p (Ci1k 2  Ci 2k1 )  p 2Ci 2k 2 ] Ak  0 .

(3.4)

k 1

Then, the stresses  ij are expressed as
3



i 2  2 Re   Li f  ( z ) ,  i1  2 Re  Li p  f  ( z  ) ,
 1




(3.5)

where the matrix Li is given by
3

Li   (Ci 2 k1  pCi 2 k 2 ) Ak .

(3.6)

k 1

For a two-dimensional problem, i.e. with geometry and external loading invariant in the direction
normal to the plane ( x, y) , any solution of anisotropic elasticity problems can be represented with
an analytic function vector f (z ) which is defined as
f ( z )   f1 ( z ),

f 2 ( z ),
32

f3 ( z )

T

(3.7)

with the argument z  x  py in the generic form. Once the analytic solution of f (z ) is determined
for a boundary value problem, the argument z should be replaced with z  for each component
function to compute the field quantities from the equations (3.3) and (3.5).
Consider a semi-infinite traction-free crack problem in an anisotropic solid, as shown in
Figure 3.1, the crack tip is located at the origin, and the traction-free crack lies along the negative
x -axis (i.e. axis 1).

Figure 3.1 Crack-tip coordinate system of a semi-infinite crack problem

The boundary conditions consist of the continuity condition of traction along the entire x -axis
and the traction-free condition on the semi-infinite crack, which leads to a homogeneous Hilbert
arc problem [109]. Solving the problem, an asymptotic expansion of the crack-tip elastic fields
can be obtained as
Lf ( z ) 

N

N

n 0

n0


1  g( z )
 i h( z ) 

2 z


(3.8)

where g( z )   a n z n and h( z )   b n z n as N   ,and both a n and b n are column vectors with
three real-number elements. The detailed derivation of this asymptotic expansion is presented in
Appendix B. The expression in (3.8) can be rewritten in a compact form as
n
1 2 N 1
f ( z )   i[1( 1) ] / 2L1c n z
2 n0

33

n1
2

,

(3.9)

where c 2n  a n , and c 2n1  b n .The first term of coefficient series c 0 is related to stress intensity
factors as
k  2 c 0 ,

(3.10)

where the stress intensity vector k is defined as k  {K II , K I , K III }T and K I K II K III respectively
represents the mode-I, mode-II and mode-III stress intensity factor. Under small-scale
yield/bridging condition, the stress intensity vector k is related to strain energy release rate

G

[106] as
G

1 T
k Hk ,
4

(3.11)

where the matrix H is defined with a positive-definite hermitian matrix B  iAL1 as
HBB.

(3.12)

In the present study, the case of orthotropic solids under plane-stress condition is investigated.
Consider only in-plane deformation fields, the stress-strain relationship can be written as
 xx  c11 c12 c16    xx 
  


 yy   c12 c22 c26    yy  ,
   c


 xy   16 c26 c66  2 xy 

(3.13)

using the Voigt notation of elastic stiffness tensor c ij . Then, the eigenvalue problem in (3.4) is
reduced to [110]

Q  p(R  R

T



)  p 2T A  0 ,

(3.14)

where the 2x2 matrices Q, R and T are
c
Q   11
c16

c16 
c
, R   16

c 66 
c66

34

c12 
c
, T   66

c 26 
c 26

c 26 
.
c 22 

Now the matrices A, L, B and H are 2x2 matrices, and the field quantities in (3.3) and (3.5) can
be determined by two analytic functions f1 ( z1 ) and f 2 ( z 2 ) . By integrating (3.9), the analytic
functions of an asymptotic crack-tip field can be expressed as

f ( z ) 

n
2 N 1 [1( 1) ] / 2


n0

i

n 1

L

1 ( n )
1 1

c



 L12c2( n ) z

n1
2

(  1 or 2)

(3.15)

1
where L ( ,   1 or 2) are the elements of L1 , c( n) (  1 or 2) are the elements of c n .

Eventually, the in-plane displacement components near a crack tip in general anisotropic solids
can be written as
u

2 N 1



n 0

v

2 N 1



n 0

n 1
n 1
n 1
n 1
n

2


1
1
Re i[1 ( 1) ] / 2 ( A11L11
z1 2  A12L211 z2 2 )c1( n )  ( A11L12
z1 2  A12L221 z2 2 )c2( n ) 
n 1 


n 1
n 1
n 1
n 1
n

2


1
1
Re i[1 ( 1) ] / 2 ( A21L11
z1 2  A22L211 z2 2 )c1( n )  ( A21L12
z1 2  A22L221 z2 2 )c2( n ) 
n 1 



(3.16)

In (3.16), the crack-tip displacement fields including higher order terms are expressed in an
asymptotic expansion similar to the Williams-type expansion for isotropic cases.
In the present study, the fracture problem can be simplified as an orthotropic material under
plane-stress condition, with the crack surface lying on one plane of symmetry. In this case,
explicit expressions of matrices A, L, B and H are derivable and presented in Appendix C.

3.2.2 Least-squares method
Given a set of displacement data (u~i , v~i ) measured at points ( xi , yi ) with i  1, 2,, M , the
coefficients c1( n ) and c2( n ) in the asymptotic expansion in (3.16) can be determined using leastsquares methods [34], [37] . In linear least-squares method, the crack-tip location must be known
a priori [34]. However, for fiber-reinforced composites, identification of the crack-tip location is
35

extremely difficult due to the generation of relatively large crack-tip damage zone during
fracture process. Therefore, in the present study, the nonlinear least-squares method is used to
simultaneously determine not only stress intensity factors but also crack-tip location, higherorder terms in the asymptotic expansion of displacement fields and rigid-body displacement
components [37].
The computational procedure is briefly described below. Firstly, for convenience, the general
form of displacement components in (3.16) is rewritten in linear combinations as
u ( x, y ) 

4 N 3

 Pn ( x, y)dn

(3.17)

n 0

v ( x, y ) 

4 N 3

 Qn ( x, y)dn

(3.18)

n 0

where d 2n  c1( n) and d 2n1  c2( n ) . Then, the sum of squared

residuals for

M data

points( i  1, 2,, M ) is expressed using matrix notation as
2
2
R  [ Pin ]{d n }  {u~i }  [Qin ]{d n }  {v~i }

,

(3.19)

where the elements of matrices [ Pin ] and [Qin ] are defined as
Pin  Pn ( xi , yi ) , Qin  Qn ( xi , yi ) .

(3.20)

To minimize the sum of the squared residuals, let
R
0.
d n

(3.21)

Then a set of linear equations can be obtained as

[P

d   [P

T
T
im ] [ Pin ]  [Qim ] [Qin ]

n

T
im ]

u~i  [Qim ]T v~i 

(3.22)

As long as the number of the data points is greater than that of the unknown coefficients, the
coefficients can be determined by solving (3.22) in a least-squares sense.

36

In nonlinear least-squares method, the determination of crack-tip position ( x0 , y0 ) , which is
relevant to an arbitrary Cartesian coordinate system of x  and y  , as shown in Figure 3.2, is
included in the computational procedure as additional unknown parameters.

y

( x0 , y0 )

x

y

x

O

Figure 3.2 An arbitrary Cartesian coordinate system x  and y  with the crack-tip location ( x0 , y0 )

Since the coordinate variables

x

and y are relevant to a coordinate system with its origin

being at the crack-tip, then a coordinate transformation is made in (3.17) and (3.18) as

x  x  x0 ,

y  y  y0 .

(3.23)

As a result, the equations (3.17) and (3.18) can be represented as:

u ( x, y) 

4 N 3

 P ( x  x , y  y )d
n 0

v( x, y) 

n

0

0

n

(3.24)

4 N 3

 Q ( x  x , y  y )d
n0

n

0

37

0

n

(3.25)

where x0 , y0 and d n are unknowns. This is a nonlinear least-squares problem that can be solved
iteratively with a nonlinear optimization algorithm. In this paper, the problem is solved by a
separable nonlinear least-squares method, which separates the determinations of crack-tip
location and coefficient series. In details, the crack-tip location ( x0 , y0 ) is updated iteratively
using a nonlinear optimization algorithm, while at each iteration the coefficient series are
obtained by solving the linear equations in (3.22) in a linear least-squares sense with the
estimated crack-tip location. The flowchart of the separable nonlinear least-squares method is
shown in Figure 3.3. The nonlinear optimization method used to set new trial values of crack-tip
location ( x0 , y0 ) is Levenberg-Marquardt algorithm [111], [112]. Finally, the determined crack-tip
location is interpreted as the location of LEFM effective crack tip to which the analytic
continuation of elastic far field converges. This effective crack-tip position is equal to the
effective crack length in the case of self-similar crack growth if the coordinate origin is located at
the notch mouth.

38

Figure 3.3 Flowchart of the separable nonlinear least-squares method for the determination of
coefficients of asymptotic expansion and effective crack-tip position.

In the study, the crack-tip displacement fields are measured by a 2-D optical method—Digital
Image Correlation. The rigid body motions of specimen during the fracture test may influence
the accuracy of in-plane displacement measurement, and hence influence the estimated crack-tip
field parameters. Therefore, the terms representing the effects of rigid body motions must be
39

added to the analytical expressions of displacement fields in (3.24) and (3.25). The rigid body
motions can be categorized as in-plane and out-of-plane motions.
Schematic diagrams for in-plane rigid body translation and rotation of specimen are shown in
Figure 3.4 (a) and (b), respectively.

Tx
Tx

Rot

Ty
Tx

y'

y'
o

o

x'

x'

(a)
(b)
Figure 3.4 Schematic diagrams of the in-plane rigid body motions of specimen: (a) in-plane rigid
body translation (b) in-plane rigid body rotation
Firstly, the in-plane rigid body translations Tx and Ty must be added to the expressions of
displacement fields for least-squares method. On the other hand, the in-plane rigid body rotation

Rot can cause additional in-plane displacements as follows:
uRot   y  Rot

(3.26)

vRot  x  Rot

(3.27)

In the previous studies [37], these two additional displacements were added to the analytical
displacement fields. However, it is found that the term representing the effect of in-plane rigid
body rotation is already included in the asymptotic expansions of displacement fields in the
equation (3.16) as follows,









1
1
Rot   Re i A11L12
p1  A12L221 p2 c2(1)  Re i A21L12
 A22L221 c2(1)
1

40

(3.28)

Therefore, no additional terms corresponding to the in-plane rigid body rotation need to be added
to the expressions of displacement fields.
Schematic diagrams for the out-of-plane motions of specimen are shown in Figure 3.5 (a)-(c),
including out-of-plane translation, out-of-plane rotation with respect to x axis and rotation with
respect to y axis.


y'

y'

y'

x'
o
(a)

x'

x'
o

o
(b)

(c)

Figure 3.5 Schematic diagrams of the out-of-plane rigid body motions of specimen: (a) out-ofplane translation (b) out-of-plane rotation with respect to y axis (c) out-of-plane rotation with
respect to x axis. The original position of specimen is represented by solid lines, while the outof-plane motion is represented by dashed lines.

The effect of out-of-plane translation is negligible when the distance between the specimen
surface and the image-capturing device is large. Besides, it is found that the term which
represents the out-of-plane rotation with respect to y axis is also included in the asymptotic
expansions of displacement fields in the equation (3.16), similar to the case of in-plane rigid
body rotation. On the other hand, the effect of the out-of-plane rotation with respect to x axis
cannot be represented by any terms in the asymptotic expansions of displacements. Thus an
additional term must be added to the displacement fields. Eventually, the analytical expressions
of displacement fields used in the least-squares method are represented as follows,

41

u ( x, y) 

4 N 3

 P ( x' x , y' y )d
n0

v( x, y) 

n

0

4 N 3

 Q ( x' x , y' y )d
n0

n

0

 Tx

(3.29)

 Ty  y' sin 2 

(3.30)

0

0

n

n

where  denotes the out-of-plane rotation angle with respect to x axis, Tx and Ty respectively
represent the in-plane rigid body translations in x and y directions.

3.2.3 Conservation integrals
In parallel to the least-squares method, the stress intensity factors, energy release rate and cracktip location can also be extracted using conservation integrals. The conservation integrals are
path-independent, i.e. the integrals yield the same value even if different integral paths are used.
Figure 3.6 shows an arbitrary integral path of the conservation integrals. The integral path 
always starts from the lower side of crack faces and ends at the upper side of crack surface. The
conservation integrals used in this study are the J-integral and M-integral.



y

x

(x0, 0)

Figure 3.6 Integral paths of conservation integrals
42

The J-integral [40] in linear elastic solids is expressed as
1

J    u, n1  u,1n ds
2



where u ,  is u / x , n is the outward normal to



(3.31)

and ,   1 or 2 .Here the summation

convention is assumed for repeated indices. The J-integral is equal to the strain energy release
rate in elastic materials, i.e.
J G

1 T
k Hk .
4

(3.32)

For a pure mode-I case, the stress intensity factor can be determined by assuming k  0, K I T .
To determine the effective crack-tip position, the 2-dimensional M-integral [107], [108] is
introduced as follows
1

M     u, n  u,  n  x ds
2

 

(3.33)

where  ,  ,   1 or 2 . When evaluate M  based on a coordinate system with the origin arbitrarily
chosen on the crack plane, i.e., let the y-directional location of the crack tip be zero, the position
of the crack tip x0 can be determined by
x0 

M
,
J

(3.34)

with M  and J  being evaluated from far-field experimental data.
The conservation integrals can also be implemented by domain integral formulation [113].
Compared to the line integral, domain integral can reduce the random errors in the results by
using more data in the computation. Nevertheless, in this study the line integral is used because it
is able to provide accurate results with relatively simpler computational procedure.

43

3.3 Experimental investigation
3.3.1 Experimental procedure

The composite material investigated in this study was Cyply 1002 (formerly known as Scotchply
1002), a glass-fiber reinforced epoxy laminate, manufactured by Cytec Engineered Materials. A
total of six specimens with same geometries but different lay-up sequences were tested. All the
specimens are 3.3 mm-thick cured panels of laminates. The information about all the tested
specimens is provided in Appendix D. For the purpose of illustration, only the result of a crossply laminate [ (90 / 0)3 / 90 ]s was presented in this chapter. The orthotropic elastic properties of
the composite laminates are measured [114]–[116] and summarized in Table 3.1,which are
comparable to those available in literature [117].

Table 3.1 Elastic properties of the composite laminate, Cyply 1002
E1 [GPa]

E2 [GPa]

E45º [GPa]

12 [GPa]

12

Measured in this study

22.39

22.14

10.23

3.17

0.127

Ref. [117]

24.1

23.4

10.3

-

-

The translaminar fracture tests were conducted using the extended compact tension (ECT)
specimen geometry [118], also known as eccentrically loaded single-edge-notch tension (ESET)
specimen, following ASTM E1922 [119]. The ECT specimen is the only specimen configuration
that has become standardized for translaminar testing of composite laminates. Furthermore, the
ECT specimen exhibits stable crack growth, which is critical in the measurement of R-curves

44

[33]. The specimen was machined using an abrasive water jet cutter. The dimensions of the
specimen are shown in Figure 3.7 (a). The measured width of the specimen is 25.44 mm. The
width of initial notch is 0.20mm, narrow enough to be considered as a slim crack. The measured
original crack length is 12.93 mm.
To conduct full-field displacement measurements using DIC technique, a random speckle
pattern was spray-painted on the specimen. A monochrome CCD camera (1280 x 960 pixels)
was used to capture digital images of the speckle pattern. The size of the field-of-view was 28.8
mm x 21.6 mm (22.5 m/pixel) centered at the initial notch tip. An image of the speckle pattern
in un-deformed configuration is shown in Figure 3.7 (b).
The fracture test was conducted at a loading rate of 0.01 mm/sec using MTS servo-hydraulic
test frame. The extension rate is slow enough to cause a quasi-static crack growth, which is the
prerequisite of fracture toughness measurement. The loading direction lies along the y-axis (or
axis 2). The applied load and load-point displacement were recorded with a 10-Hz sampling rate.
A picture of the unpainted side of the specimen after the test is shown in Figure 3.7(c). The
picture confirms that the translaminar crack propagated in a self-similar manner. Moreover, no
extensive delamination or longitudinal splitting was found through the observation of fracture
surface.

45

Figure 3.7 (a) Configuration of ECT (ESET) specimen (b) Speckle patterns for DIC
measurement (c) Unpainted side of the specimen after test. In (a), W=25.44mm. The width of
notch is 0.20mm.The measured original crack length is 12.93 mm. In (b), the area surrounded by
the outer rectangle is the region of interest in which the displacements are measured by DIC, the
area within the polygon is the region from which data points used in least-squares methods are
sampled. In (c), it shows self-similar crack growth.

During the fracture test, a series of digital images of the speckle pattern was acquired with a
frame rate of 1 Hz. A total of 105 images were analyzed by using 2-D DIC software. The DIC
software was custom-built by implementing the Newton–Raphson iteration 2-D DIC algorithm
[97]. The quasi-Newton optimization method in the original paper was replaced by a more robust
optimization method--Levenberg-Marquardt method. The precision of the software was
estimated to be 0.02 pixels, through a numerical test using synthetic speckle patterns. A subsetsize of 41 x 41 pixels was used in the DIC analysis, and the spacing between two neighboring
subsets was 10 pixels.

46

3.3.2 Experimental results
3.3.2.1 Load-CMOD curve and full-field displacements

A load-crack mouth opening displacement (CMOD) curve obtained during the fracture test is
shown in Figure 3.8. The crack mouth opening displacement is measured by DIC. The loadCMOD relationship is linear up to approximate 70% of the peak load, and subsequently shows a
gradual deviation from the linear trend which represents the generation of a damage zone around
the crack tip. The deviation from the linear trend satisfies the requirement d n / d n0  0.3 , thus
the fracture toughness K IC can be obtained using the stress intensity factor formula provided in
ASTM standard E1922. After reaching the peak load 1185N, the applied load decreases
smoothly without discrete jumps, which means that the translaminar crack growth took place in a
macroscopically stable manner.

47

1400

d

dn-0

1200

n

Pmax=1185N

Load (N)

1000

800

600

400

200

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

CMOD (mm)

Figure 3.8 Load versus Crack Mouth Opening Displacement (CMOD) Plot

In this study, the frame of DIC measurement region was depicted by the outer white rectangle
in Figure 3.7(b). An example of the measured crack-tip displacement contour maps is shown in
Figure 3.9. This example set of displacement data corresponds to the moment when the CMOD
equals 0.30 mm. It can be seen that the symmetric U field and anti-symmetric V field with
respect to the crack plane are consistent with the nature of mode I fracture. Besides, because the
DIC subset near the crack surface contains displacement discontinuity so that the fundamental
assumption of DIC is violated, the displacement data near the crack surface was considered to be
invalid and discarded. In the estimation of the crack-tip field parameters, the displacement data
was used without additional smoothing.

48

-0.12

8

-0. 13

6
-0.14

4

-0. 15

y (mm)

2

-0.17
-0.18

-0.16

0
-0. 18

-2

-0. 17

-0. 16

-4

-0. 15

-6

-0.14
-0.13

-8
0

5

10

15

20

25

x (mm)
(a)

-0. 26

5
-0.
2

-0. 24

-0.22

-0. 20

-0.18

-0.17

-0.15

8
-0. 2

4

-0. 27

-0.11

6

-0.13

8

y (mm)

2
0
-2
27

6

-0
.

-0.
2

9
-0. 2

-0. 31

-0.32

-0.34

-0.36

-0.38

-0.39

-4
-6
-8
0

5

10

15

20

25

x (mm)
(b)
Figure 3.9 An example set of displacement fields measured by DIC: (a) U field (b) V field. It
was measured at the moment when the CMOD equals 0.30 mm. The unit of displacements is the
millimeter (mm).
49

3.3.2.2 Results of least-squares method

For least-squares method, the accuracy of estimated crack-tip field parameters is influenced by
several factors, including the number of terms in asymptotic expansion, the domain from which
the displacement data are sampled, and the out-of-plane rotation of the specimen. Before
analyzing the whole fracture test, two displacement data sets, measured at the moments when the
CMOD equals 0.30mm and 0.75 mm, were chosen to investigate the influence of these factors.
As was previously mentioned, the displacement data measured near the crack face and around
the crack tip is unreliable because of the displacement discontinuity and the large deformation
within the area, respectively. Therefore, the data points used in the estimation of crack-tip field
parameters were sampled from a square domain centered at the crack tip with an interior square
area being deleted, as shown in Figure 3.7 (b).
Firstly, the influence of the number of terms in asymptotic expansion on the determination of
stress intensity factor and effective crack length was investigated. For the data domain used in
the investigation, the half width of outer square domain was set to 350 pixels (7.88 mm), and the
half width of interior deleted square area was set to 100 pixels (2.25 mm) to avoid threedimensional effects [42]. The variations of determined stress intensity factor and effective crack
length with respect to the number of terms are respectively shown in Figure 3.10 (a) and (b). It is
observed that, the determined KI and effective crack length first fluctuate with increasing number
of terms and then stabilize after 4 terms. Therefore, in all the following investigations, 10 terms
of asymptotic expansion are used to ensure stable results.

50

(a)
18

Effective crack length (mm)

16
14
12
10
8
6
4
CMOD=0.75 mm
CMOD=0.30 mm

2

0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Number of terms
(b)
Figure 3.10 Variation of the (a) stress intensity factor and (b) effective crack length obtained by
the least-squares method with respect to the number of terms in crack-tip asymptotic expansion.
51

Secondly, the influence of the data domain used in the parameter estimation was investigated.
The investigation was conducted with the half width of outer square domain being fixed to 350
pixels, while the half width of interior deleted area was varied from 100 pixels to 300 pixels.
Figure 3.11(a) and (b) respectively shows the variations of stress intensity factor and effective
crack length with respect to W—the half width of interior deleted area. It is seen that the results
of KI and crack length remains stable before W gets close to 300 pixels. It means that the result is
not sensitive to the choice of data domain. Therefore, W is set to be 100 pixels for all the
following investigations to ensure that enough data points are used.

(a)
Figure 3.11 Variation of the (a) stress intensity factors and (b) effective crack length obtained by
the least-squares method with respect to the half length of interior deleted area in the data
domain, which is represented by W.

52

Figure 3.11 (cont’d)
18

Effective crack length (mm)

16
14
12
10
8
6
4
CMOD=0.75 mm
CMOD=0.30 mm

2
0
100

150

200

250

300

W (pixels)
(b)
As a confirmation of the least-squares fitting, Figure 3.12 (a) compares the experimental
displacement field in the tension direction with the displacement field that is regenerated with the
results of least-squares method. The data set is measured when CMOD equals 0.30 mm. It is
clearly observed that the regenerated displacement contours match very well with the
experimental data.

53

6
4

y (mm)

2
0
-2
-4
-6
5

10

15

20

15

20

x (mm)
(a)

6
4

y (mm)

2
0
-2
-4
-6
5

10

x (mm)
(b)

Figure 3.12 Comparison of the V-field contours measured by DIC and regenerated from leastsquares method at the moment when CMOD is 0.30 mm. The results of least-squares methods
are obtained (a) without the out-of-plane rotation term and (b) with the out-of-plane rotation term.
The DIC measurement is represented by solid lines, while the reconstructed field is represented
by dashed lines.
54

As previously stated in Section 3.2.1, the out-of-plane rotation of the specimen will influence
the accuracy of in-plane displacement measurement and thus influence the accuracy of crack-tip
field parameters. As an illustration, Figure 3.12 (b) shows the comparison of experimental and
regenerated displacement fields when the out-of-plane rotation is not considered, i.e., the term
containing  in the equation (3.30) is not used in the computation. In this case, it can be seen that
the regenerated field does not match well with the experimental field. Therefore, the inclusion of
out-of-plane rotation term in the least-squares fitting process improves the fitting results.
Then, the least-squares method was applied to process all the displacement data obtained
during fracture test. Figure 3.13 (a) and (b) respectively shows the evolution of stress intensity
factor KI and effective crack length as a function of CMOD. In Figure 3.13 (a), it is observed that
KI increases monotonically to a steady-state value of about 20 MPa m . As a confirmation, KI
estimated from the least-squares method was compared with that determined from the applied
load and specimen geometry. In isotropic cases, the stress intensity factor (SIF) formula for
extended compact tension (ECT) specimen [118], [119] is given as
KI 

  a / W , F ( )   1/ 2

P
F ( ) ,
BW 1 / 2

(3.35)

[1.4   ]
[3.97  10.88  26.25 2  38.9 3  30.15 4  9.27 5 ]
3/ 2
[1   ]

where P is the applied load, B is the specimen thickness, W is the width of the specimen, and a is
the crack length. Based on the work of Bao et al. [12], [13], the influence of material orthotropy
on KI is found to be negligible for this type of composite laminates. From Figure 3.13 (a), an
excellent agreement is found between the KI values estimated by least-squares method and SIF
formula with original crack length at the beginning stage of the fracture test, which validate the

55

applicability of least-squares method. However, during the steady-state crack propagation, the
fracture toughness KIC estimated by the equation (3.35) using the original crack length is about
16 MPa m , only 80 percent of the KI value estimated by least-squares method. The inconsistency
between two methods is due to the increase of effective crack length, which makes the SIF
formula with original crack length invalid. In Figure 3.13 (b), it is observed that the increase of
effective crack length can be divided into 3 different stages with different slopes, which are
caused by different fracture damage mechanisms. These fracture damage mechanisms will be
discussed in details in the characterization of R-curve behaviors.
Figure 3.13 (a) and (b) also shows the influence of out-of-plane rotation on the estimation of
stress intensity factor and effective crack length by the least-squares method. When the out-ofplane rotation is not considered, the discrepancy between KI obtained by SIF formula and K I
obtained by least-squares method becomes larger than the discrepancy when the out-of-plane
rotation is considered. Meanwhile, the effective crack length is significantly underestimated from
the original crack length at the beginning stage of the fracture test. Actually, the maximum of
out-of-plane rotation angle  estimated by least-squares method is as small as 1 degree.
Therefore, it can be seen that the least-squares method is very sensitive to the out-of-plane
rotation. Since it is very difficult, if not impossible, to eliminate the small out-of-plane rotation
experimentally, inclusion of the out-of-plane rotation term in the computational process is a
better choice.

56

(a)

20

Effective crack length (mm)

19
Stage III

18
17
16
Stage II
15
14

Stage I

13
12
With the out-of-plane rotation term
Without the out-of-plane rotation term

11
10
0

0.2

0.4

0.6

0.8

1

1.2

CMOD(mm)

(b)
Figure 3.13 Variation of the (a) Mode I stress intensity factor KI and (b) effective crack length
estimated by the least-squares method with respect to CMOD. The effect of out-of-plane rotation
is shown as well.

57

3.3.2.3 Results of conservation integrals

For the computation of conservation integrals, the displacement gradients were obtained from
the displacement data using a numerical differentiation method [120] with the computational
window size being 5 by 5. A circular contour containing a total of 600 data points was used as
the integral path. The data points on the integral path were obtained using the spline interpolation.
Path integral was calculated following the Simpson’s rule. In the computation, the displacement
gradients near the crack face were eliminated because the displacement data within this area are
unreliable. The error caused by the elimination of data is negligible because in mode I fracture
problems the displacement gradients near the crack-face are very close to zero.
The in-plane rigid body translation and rotation won’t cause any changes in the results of
conservation integrals. On the other hand, as same as for the least-squares method, the out-ofplane rotation of the specimen can influence the results of conservation integrals. However, the
effect of the out-of-plane rotation cannot be eliminated within the computational process of
conservation integrals. In this study, the effect of out-of-plane rotation is eliminated from the
experimental displacement data prior to the computation of conservation integrals, using the outof-plane rotation angle  estimated by least-squares method.
To validate the applicability of conservation integrals, the path-independency of integrals was
investigated before analyzing the whole fracture test. The two displacement data sets chosen for
least-squares method, which were obtained at the moments when the CMOD equals 0.30 mm
and 0.75 mm, were also chosen for the investigation of conservation integrals. For each data set,
the conservation integrals were computed with the radius of integral path varying from 100 to
300 pixels. Figure 3.14 (a) and (b) show the variations of J-integral value which equals energy
release rate and effective crack length with respect to the radius of integral path. In Figure 3.14
58

(a), when the CMOD equals 0.30 mm, J-integrals are stable in spite of the fluctuations which are
caused by the noise in displacement data. However, when the CMOD equals 0.75 mm, J-integral
first increases and then converges when the radius is larger than 200 pixels. The underestimation
of J-integral with small radii may be due to the relatively larger influence of the crack-tip
fracture damage zone. In Figure 3.14 (b), the estimated crack lengths are found to be stable with
respect to different radii. Based on the above investigation, the conservation integrals in the
following studies are first computed with 10 different radiuses from 200 to 300 pixels. Then the
average is used as the representative result.
The influence of the out-of-plane rotation on the computation of conservation integrals was
also shown in Figure 3.14 (a) and (b).When the effect of out-of-plane rotation is not eliminated,
the J-integral is overestimated by about 20 percent for the data set when the CMOD equals 0.30
mm. The overestimation is 5 percent when the CMOD equals 0.75 mm. Moreover, the estimated
crack length becomes path-dependent when the CMOD equals 0.75 mm. It can be seen that the
elimination of out-of-plane rotation effect before computation is necessary.

59

(a)

(b)
Figure 3.14 Effect of radius of integral path on the estimation of the (a) J-integral and (b)
effective crack length obtained by conservation integrals.
60

Analyzing all the experimental data for the fracture test, the variations of J-integral and
effective crack length were plotted as functions of CMOD in Figure 3.15 (a) and (b). In Figure
3.15 (a), it is observed that the J-integral is parabolic at the beginning stage of fracture test,
which is consistent with the relationship between J and K in the equation (3.32) since K varies
linearly at this stage. The J-integral increases to a steady-state value of about 25 KJ/m2, which is
close to 26 KJ/m2, the value predicted by using the steady-state value of KI—20 MPa m and the
equation (3.32). In Figure 3.15 (b), it can be seen that the effective crack growth also share the
same general trend of 3 stages as that was obtained by the least-squares method.
Figure 3.15 also shows the estimated J-integral and effective crack length when the effect of
out-of-plane rotation is not eliminated. In this case, the J-integral is significantly overestimated at
the early stage before reaching the steady state. The effective crack length is underestimated at
the beginning stage of fracture test, and is overestimated during the steady-state crack growth.

61

30

J-Integral (KJ/m2)

25

20

15

10

5
After correcting the out-of-plane rotation
Before correcting the out-of-plane rotation
0

0

0.2

0.4

0.6

0.8

1

1.2

CMOD (mm)

(a)

20

Effective crack length, M/J (mm)

19
18
17
16
15
14
13
12
After correcting the out-of-plane rotation
Before correcting the out-of-plane rotation

11
10

0

0.2

0.4

0.6

0.8

1

1.2

CMOD (mm)

(b)
Figure 3.15 Variation of the (a) J-integral and (b) effective crack length obtained by conservation
integrals with respect to CMOD. The effect of out-of-plane rotation is shown as well.
62

3.3.2.4 Comparison of the results of least-squares and conservation integrals

Now the results of the least-squares method and conservation integrals were compared to
validate the small-scale yielding condition. The variations of stress intensity factor and effective
crack length as functions of CMOD were plotted together in Figure 3.16 (a) and (b). The term “JIntegral” in Figure 3.16 (a) denotes KI obtained by using the equation (3.32) while assuming the
mode II stress intensity factor KII equals zero.

(a)
Figure 3.16 Comparison of the (a) mode I stress intensity factor and (b) effective crack length
estimated by different methods.

63

Figure 3.16 (cont’d)
20

Effective crack length (mm)

19
18
17
16
15
14
13
12
Least-squares method
Conservation integrals

11
10

0

0.2

0.4

0.6

0.8

1

1.2

CMOD(mm)

(b)

First of all, it is observed that the KI obtained by the least-squares method and that by the
conservation integrals show a good agreement when load increases linearly proportional to
CMOD. Then the difference between the two results becomes slightly larger. However, the
maximal relative difference between the results of two methods during the steady state is still
less than 3 percent. Furthermore, the effective crack lengths obtained by using the least-squares
method and conservation integrals are consistent only except at the early stage of the fracture test
when the signal-to-noise ratio is too low to get reliable results. The consistency of the estimated
crack-tip field parameters by the least-squares method and conservation integrals represents the
equivalence of the stress intensity factor approach and the energy release rate consideration.
Thus, the small-scale yielding condition and the applicability of LEFM are validated.

64

As an additional confirmation, KI was calculated by using the SIF formula in the equation
(3.35), with the original crack length a being replaced by the effective crack length obtained by
least-squares method. The result is shown as the dotted line in Figure 3.16 (a). In this case, the K I
estimated by using the equation (3.35) can obtain a good estimation of KI for the steady-state
crack propagation, if the fluctuation of the result is neglected. Therefore, it means that the SIF
formula provided in ASTM standard E1922 can be extended to estimate the toughness of steadystate crack propagation if the effective crack length is used instead of the original crack length.

3.3.2.5 R-curve behaviors
Since the small-scale yielding condition has been validated, the crack-growth resistance curves
(R-curves) can be applied to characterize the fracture behaviors of the composites. Based on the
variation history of stress intensity factor and effective crack length, the R-curves was
constructed in terms of stress intensity factors (KR) and shown in Figure 3.17. In engineering
design, the R-curves are usually used to predict load-carrying capability of the structure. For the
ECT specimen, a sequence of driving force curves (in terms of KI versus crack length) were
constructed corresponding to different loads by using the SIF formula in the equation (3.35). The
driving force curve that is tangent to the KR-curve was plotted together with the R-curve in
Figure 3.17. It is found that the load corresponding to this curve is 1160 N. This predicted
maximal load is very close to the experimental maximal load—1185 N. It is a proof of the
accuracy of estimated R-curves.

65

Figure 3.17 R-curves in terms of stress intensity factor, determined by the least-squares method
and conservation integrals.

R-curve was also plotted in terms of energy release rate in Figure 3.18 (JR). While KR curve
shows good promise for materials where small-scale yielding occurs in front of the crack tip, the
JR curve can be applied for more general cases, including the large-scale yielding conditions. The
obtained KR and JR curves demonstrate two distinct stages of rising R-curves with two different
slopes before they reaches a steady state. The effective crack length increases by 2 mm before
reaching the steady-state, which means that the length of fracture damage zone is larger than 2
mm.

66

30

25

JR (KJ/m2)

20

15

10

5
Least-squares method
Conservation integrals
0
12

13

14

15

16

17

18

19

20

Effective crack length (mm)

Figure 3.18 R-curves in terms of strain energy release rate, determined by the least-squares
method and conservation integrals.

This experimental observation is consistent with the reported R-curve behaviors of
translaminar crack growth in composite laminates [121], [122]. It is believed that two stages of
rising R-curves before reaching the steady state correspond to different fracture mechanisms. In
the rapidly rising stage, the fracture process starts by forming a small diffuse damage zone that
consists of multiple micro-cracks due to matrix cracking and fiber/matrix debonding. In the
relatively slowly rising stage, the fracture process creates a macroscopic fracture damage zone
due to the fiber bridging and fiber pull-out. A schematic diagram of these fracture mechanisms is
shown in Figure 3.19. Moreover, for the quasi-static fracture problem, the steady-state crack
propagation can be considered as the steady-state propagation of fully-generated fracture damage
zone.
67

Figure 3.19 Fracture mechanisms in the translaminar fracture of fiber-reinforced composites

3.4 Discussion
In this study, only the mode-I problem was investigated. However, it is noteworthy that both
of the least-squares method and conservation integrals have the potential to be extended to study
the mode-II or mixed mode problems. In the least-squares method, the mode-I and mode-II stress
intensity factors are obtained simultaneously. In parallel, the interaction J-integral of
conservation integral can be used to separately determine the mode-I and mode-II stress intensity
factors [41]. The interaction J-integral of two elastic fields S and Sˆ is defined as

J int [ S , Sˆ ]  J [ S  Sˆ ]  J [ S ]  J [ Sˆ ]

    uˆ , n1    uˆ ,1n  ˆ  u ,1n ds


68

(3.36)

From the equation (3.32), the interaction J-integral of two linear elastic crack-tip fields with



stress intensity vectors k  K II , K I T and kˆ  Kˆ II , Kˆ I



T

can be expressed as

1
J int[ S, Sˆ ]  k T H kˆ
2

(3.37)

If the two auxiliary crack-tip fields Sˆ I and Sˆ II are analytically constructed with kˆ I  0, 1T and
T
kˆ II  1, 0 respectively, the stress intensity vector of a measured crack-tip field S can be

determined by
 J int[ S, SˆII ]
K 
k   II   2H 1  int

 J [ S, SˆI ] 
 KI 

(3.38)

The experimental validation of conservation integrals has been carried out for isotropic materials
[36], [123]. For fiber-reinforce composite laminates, the only difference is that the auxiliary
crack-tip fields Sˆ I and Sˆ II need to be analytically constructed using the asymptotic expansion of
crack-tip fields in anisotropic solids.
Figure 3.20 shows the mode I and mode II stress intensity factors of this problem estimated by
least-squares method and interaction J-integral. It can be seen that the two methods have
obtained nearly identical results, which also demonstrates the equivalence of stress intensity
factor approach and energy release rate approach, and hence the validity of small-scale yielding
condition. The mode II stress intensity factor is negligible compared to the mode I stress
intensity factor, which is consistent with the nature of pure mode-I problem.

69

(a)

(b)
Figure 3.20 Comparison of (a) mode I and (b) mode II stress intensity factors estimated by leastsquares method and interaction J-integral.

70

3.5 Summary
In this chapter, an experimental investigation based on full-field optical displacement
measurements was presented as a way to characterize the subcritical and steady-state crack
advances in translaminar fracture of cross-ply composite laminates. Under the assumption of
homogeneous approximation, the experimental investigation relied on the linear elastic fracture
mechanics theory in anisotropic solids. Two methods of least-squares method and conservation
integrals were proposed to measure the crack-tip field parameters (i.e. stress intensity factor,
energy release rate and effective crack length) and R-curves from the crack-tip displacement
fields measured by digital image correlation. The sources of error of two methods were analyzed,
with an emphasis on the influence of out-of-plane rotation of specimen. Eventually, accurate Rcurve measurement of translaminar fracture of composite laminates was achieved with
experimental convenience. This study provides a foundation for the top-down approach to multiscale analysis of translaminar fracture in fiber-reinforced composite laminates.

71

Chapter 4
Inverse Extraction of Cohesive-zone Law for Fiberreinforced Composite Laminates Part I: Numerical Analysis
4.1 Introduction
In this chapter, two analytical inverse methods are developed to extract the cohesive-zone laws
(CZLs) from elastic far-field displacements surrounding a crack-tip cohesive zone, for the
translaminar fracture of fiber-reinforced composite laminates. In the study, the bulk material
property of fiber-reinforced composite laminates is approximated as elastic anisotropy.
The chapter is organized as follows. In Section 4.2, all the analytical derivations for the
establishment of inverse methods are provided. Firstly, an eigenfunction expansion of the planar
elastic field around the crack-tip cohesive zone for anisotropic solids is developed. The
displacements, displacement gradients and stresses around and within the cohesive zone can be
represented in terms of the eigenfunction expansion. Thus, the CZL which represents the
relationship between cohesive traction and separation can also be represented with the
eigenfunction expansion. Next, two analytical inverse methods are developed. The first method
is named “field projection method (FPM)”, in which the coefficients of eigenfunction expansion
are extracted using the interaction J-integrals between experimental fields and auxiliary probing
fields. As stated in Chapter 2, the basis of this FPM method was proposed in the work of Hong
and Kim [75] for isotropic materials. In this study, the method is extended to the anisotropic
materials. In parallel, another method based on the nonlinear least-squares method is developed

72

to determine the coefficients, by fitting the analytical expressions of crack-tip displacements to
the experimental data in a least-squares sense.
After establishing the scheme of analytical inverse methods, a set of numerical tests are
conducted to assess the applicability and error sensitivity of the two inverse methods, as well as
to provide guidelines for experiments. In Section 4.3, the procedure of the numerical tests and
the meanings of the studied parameters are introduced. In the numerical tests, synthetic
displacement fields are firstly constructed with predefined eigenfunction expansions. Then the
two analytical inverse methods are applied to extract the CZL from synthetic displacement fields.
Several factors which can influence the accuracy of the inverse solutions are investigated,
including the shape of CZL, the noise level of inputs, the inverse distance of the data field and
the number of data points. Results of the numerical tests and further discussions about the results
are presented in Section 4.4. Finally, some concluding remarks are provided in Section 4.5.

4.2 An eigenfunction expansion of the cohesive crack-tip fields in
anisotropic solids
The derivation of the eigenfunction expansion of cohesive crack-tip fields in anisotropic solids
is also based on the anisotropic plane elasticity theory in terms of Stroh's formalism that was
described from the equation (3.1) to (3.7). However, a different boundary value problem is
solved in this chapter, which leads to an analytical solutions of f (z ) that is different from the
equation (3.8).
The boundary value problem considered in this chapter is a traction-free semi-infinite crack
with a cohesive zone of size 2c, where the center of cohesive zone is located at the origin, and
73

the crack plane is located on the negative x -axis. The schematic diagram of this problem is
shown in Figure 4.1.

Figure 4.1 Schematic diagram of semi-infinite crack with a cohesive zone of size 2c

Using traction-continuity conditions on the entire x -axis and a superposition of two stressbounded sharp-crack-tip elastic fields with respect to

x  c

and x  c , an eigenfunction

expansion of the elastic field around a cohesive zone can be derived from the solution of a
homogeneous Hilbert arc problem [109] as follows
Lf ( z ) 

1
2



z  cq( z )  z  cr( z )  i s( z )

N

N

N

n0

n0

n0



(4.1)

where q( z )   q nQn ( z ) , r ( z )   rn Rn ( z ) and s( z )   s n S n ( z ) as N   , q n , rn and s n are
column vectors with real-number elements, and Qn (z ) , Rn (z ) and Sn (z ) are entire functions.

74

Then the cohesive-zone variables can be expressed in terms of eigenfunction expansions. The
closing-traction distribution t (x) within the cohesive zone (  c  x  c ), is described as
t ( x) 

x  cq( x)

(4.2)

The jump of separation-gradients in the cohesive zone b(x) which indicates the difference of the
separation-gradients between the upper and lower faces of the cohesive zone is expressed as:
b( x)  u' ( x)  u' ( x)   H c  xr( x)

(4.3)

where H is a Stroh’s matrix defined in the equation (3.12).
Thus, the separation in the cohesive zone δ(x) can be obtained by integrating (4.3),
c

δ( x)  x b( ) d   H x c   r( ) d  .
C

(4.4)

When the closing traction t (x) and the separation δ(x) are determined, the cohesive-zone law
(CZL) t (δ) can be obtained by correlating the two parameters within the whole cohesive zone. In
addition, it is noteworthy that s n in the equation (4.1) does not contribute to the cohesive-zone
variables.
The detailed derivation of eigenfunction expansions of cohesive-crack-tip elastic fields and
cohesive-zone variables are provided in Appendix E.

4.3 Extraction of CZLs in anisotropic elastic solids
Based on the eigenfunction expansions of cohesive-crack-tip fields, two inverse methods,
namely field projection method and least-squares method, are developed to extract the CZLs
from elastic far-field displacements surrounding a crack-tip cohesive zone. More specifically, the

75

inverse methods are developed to extract the column-vector coefficients q n and rn in q( z ) and
r ( z ) . Then the cohesive traction and separation within the cohesive zone can be obtained with

the equation (4.2) and (4.4), respectively. Eventually, the CZL can be constructed by correlating
the cohesive traction with separation.

4.3.1 Field projection method
In 2003, Hong and Kim developed an inverse method named “field projection method (FPM)”
to extract the CZLs for isotropic materials [75]. The method was developed based on the Jorthogonal representation of an eigenfunction expansion of cohesive-crack-tip elastic fields,
using the far-field interaction J-integral [41] between the physical field of interest and a set of
auxiliary probing fields. This method is then extended to a nano-scale planar field projection of a
cohesive crack-tip on an interface between two anisotropic solids by Choi and Kim [78]. In the
present study, the FPM is extended to a macro-scale planar cohesive-zone problem of an
anisotropic solid.

4.3.1.1 Interaction J-integral
The FPM makes use of the interaction J-integral, which is also introduced in Section 3.4 for
the separate determination of the mode-I and mode-II stress intensity factor in mixed-mode
fracture problems. Compared to its application in the framework of LEFM, the application of
interaction J-integral in the FPM is more advanced.

76

Figure 4.2 shows the coordinate system for a cohesive crack-tip field, with S and Sˆ respectively
denoting the physical field of interest and auxiliary probing field, and  and 0 respectively
denoting the integral contour in the far field and along the cohesive-zone surface. The pathindependence of the interaction J-integral leads to the relationship:

J int[S , Sˆ ]  J int0 [S , Sˆ ]

(4.5)

.

From the equation (3.44), the interaction J-integral along the contour 0 can be represented as
J int0 [S , Sˆ ]  

c

c

t ( x)bˆ ( x)  b ( x)tˆ( x) dx   cc
T

T



c 2  x 2 qT ( x)Hrˆ ( x)  r T ( x)HTqˆ ( x)dx

.

(4.6)

Figure 4.2 Interaction J-integral paths surrounding a cohesive zone at the far field  , and along
the cohesive zone surface 0 .

4.3.1.2 Expansions of the elastic-fields in terms of orthogonal polynomial series
The FPM makes use of the orthogonality of polynominals to determine the coefficients in
q( z ) and r ( z ) . In the equation (4.6), it can be seen the interaction J-integral is expressed in terms

77

of the inner product of q(x) and r (x) with a weight function of

c 2  x 2 in  c  x  c . Thus, it

is convenient to expand q(x) and r (x) with Chebyshev polynomials of the second kind, a
complete set of orthogonal polynomials generated by the Schmidt orthogonalization of inner
product.
The Chebyshev polynomials of the second kind, denoted as U n (x) , are orthogonal in the
interval [1,1] with a weight function of 1  x 2 , i.e.

1
1



1  x 2 U m ( x)U n ( x) d x 


2

 mn

(4.7)

where  mn is the Kronecker delta.
Now, The functions q(x) and r (x) are rewritten in terms of the Chebyshev polynomials of
N
N
the second kind as q( x)   t nU n ( x) and r( x)   b nU n ( x) for ~
x  x / c and N   , where t n
n 0

n 0

and b n are column vectors with real-number elements. Thus the eigenfunction expansion of the
crack-tip elastic fields in (4.1) can be rewritten as:
Lf ( z ) 

N
N
1

z

1
t
U
(
z
)

z

1
b nU n ( z )  i s( z ) 


n n

2
n 0
n 0


(4.8)

z  z / c and N   .
for ~
Then, the distributions of the traction, separation-gradient and separation within the cohesive
zone are rewritten as:
N

t ( x)  1  x  t nU n ( x)

(4.9)

n 0

N

b( x)  H 1  x  b nU n ( x)
n 0

78

(4.10)

c

1

x

x

δ( x)   b( ) d   cH 

N

1    b nU n ( ) d 

(4.11)

n 0

for ~
x  1 and N   .
In terms of the orthogonal polynomials’ expansion, once t n and b n (n  0, 1,2,...,N) are
determined, the cohesive-zone variables t (x) , b(x) and δ(x) , and hence the CZL t (δ) can be
determined.

4.3.1.3 Traction probing and separation probing
In order to determine t n and b n (n  0, 1,2,...,N) , we will use far-field interaction J-integrals
J int[ S , Sˆ ] between a physical cohesive-crack field S and a series of auxiliary probing fields Sˆ ( n )
(n)
(n  0, 1,2,..., N) . The field S and auxiliary probing fields Sˆ are represented by the analytical





functions Lf (z )and Lfˆ ( z )

( n)

, respectively.

Firstly, the traction-probing method, i.e. the method to determine t n , is introduced as follows.





If we choose Lfˆ ( z )

( n)

such that the tractions of the auxiliary field vanish on the cohesive zone,

i.e. tˆ ( n ) ( x)  [0,0]T for  c  x  c , then from the equation (4.6) it can be seen that the interaction
J-integral has contributions only from the cohesive tractions of the physical cohesive-crack-tip
field. Furthermore, the contributions to the interaction J-integral from the cohesive-zone tractions
of the physical field S can be limited to that from a single component of the traction along x or y
direction, i.e. t1 ( x) or t2 ( x) , by changing the analytical auxiliary function. Then the coefficients t n
can be determined from the interaction J-integrals which can be calculated from the elastic farfields due to its path-independence.

79

The mathematical derivation of the traction-probing method can be described as follows. The
analytical function defining the auxiliary fields for probing the cohesive-zone tractions of the





physical cohesive-crack field S in the α-direction is denoted as Lfˆ ( z )

(n)

t

. If such functions are

chosen as

Lfˆ ( z)

(n)

t



1 ~
z  1U n (~
z )[ 1 ,  2 ]T
4

(4.12)

where n  0, 1,2,..., N , and α= 1 or 2, then the tractions and separation-gradients of the auxiliary
(n)
field Sˆt within the cohesive zone  c  x  c become



tˆ t(n ) ( x)  [0,0]T , bˆ t( n ) ( x)  1 H 1  ~
xU n ( ~
x )[ 1 ,  2 ]T .

2

(c  x  c)

(4.13)

In this case, the interaction J-integral is evaluated from the equation (4.6) as

c
J int0 [ S , Sˆt( n ) ] 
( H1 t1( n )  H 2 t2( n) )

2

t(n ) (

where n  0, 1,2,..., N , and α equals 1 or 2.

(4.14)

 1 or 2) represents the

th element

of the

column vector t n (n  0, 1,2,...,N) . H  represents the element of Stroh’s matrix H.
Because the interaction J-integral is path-independent, we can replace the integration contour
0 with a far-field contour  . Rearrange the equation (4.14), the column-vector coefficients t n can
T

be represented as follows,

t Tn  t1( n ) , t 2( n )  





T
2 1 int
H J  [ S , Sˆt(1n ) ], J int [ S , Sˆt(2n ) ] .
c

It is noteworthy that H is a symmetric matrix.

80

(4.15)





Similarly, b n can be determined with the separation-probing method. If we choose Lfˆ ( z )

( n)

such that the separation-gradient of the auxiliary field vanishes in the cohesive zone, i.e.
bˆ ( n ) ( x)  [0,0]T for  c  x  c , then the interaction J-integral has contributions only from the

cohesive zone separation-gradients of the physical cohesive-crack field. Denoting the auxiliary
fields for probing the cohesive zone separation-gradients of the physical cohesive-crack field in





the α-direction as Lfˆ ( z )

( n)

b

, we choose such functions as

Lfˆ ( z)

(n)

b



1 ~
z  1U n (~
z )[ 1 ,  2 ]T
4

(4.16)

where n  0, 1,2,..., N , and α=1 or 2, then the tractions and separation-gradients of the auxiliary
(n)
field Sˆb within the cohesive zone  c  x  c become



1
tˆ b( n ) ( x) 
1 ~
xU n ( ~
x )[ 1 ,  2 ]T , bˆ b( n ) ( x)  [0,0]T . (c  x  c)
2

(4.17)

Then the interaction J-integral is evaluated as

c
J int0 [S , Sˆb( n ) ] 
( H 1b1( n )  H 2b2( n ) )

2
where n  0, 1,2,..., N , and α equals 1 or 2.
column vector b n (n  0, 1,2,...,N) .

(4.18)

b( n ) (β=1 or 2) represents the th element of the

H  represents the element of Stroh’s matrix H.

Based on the path-independence of the interaction J-integral, the integration contour 0 can be
replaced with a far-field contour  . The equation (4.18) can be rearranged to represent the
determination method of b n as follows,

b n  b1( n ) , b2( n )  
T





T
2 1 int
H J  [ S , Sˆb(1n ) ], J int [ S , Sˆb(2n ) ]
c

81

(4.19)

After the determination of t n and b n (n  0, 1,2,...,N) , the cohesive-zone law can be obtained
through the use of the equations (4.9), (4.10) and (4.11).

4.3.2 Least-squares method
Based on the eigenfunction expansion of the cohesive crack-tip fields that was given in
Section 4.2, the coefficients t n and b n (n  0, 1,2,...,N) can also be determined using linear leastsquares method. Similar to the procedure described in Section 3.2.2, the least-squares method
extracts the coefficients from the full-field cohesive crack-tip displacements. However, in this
problem, it is difficult to represent the crack-tip displacement fields in a compact form as the
equation (3.16). Nevertheless, the computational implementation ( i.e. the computer
programming) of the data reduction scheme of linear least-squares method is still simple. The
matrix of the normal equation of least-squares method can be generated from the eigenfunction
expansion in the equation (4.8). Each element in the matrix of least-squares' normal equation
equals to the analytical displacement of one data point which is obtained with the corresponding
unknown coefficient in the eigenfunction expansion equaling unity and other coefficients
vanishing.

4.3.3 Determination of cohesive-zone size and position
In the field projection method and linear least-squares method, the size and the position of the
crack-tip cohesive zone, i.e. c and ( x0 , y0 ) ,are assumed to be known a priori. However,
identification of the cohesive-zone size and location in practical experiments is an extremely
82

hard task. Even though the orientation and the y-directional location of the cohesive zone y 0 are
visually obtainable from the images of cracked structures, determination of the x-directional
central location of cohesive zone x0 cannot rely on the visual inspection. In this study, the
determination of the size and x-directional central location of cohesive zone, i.e. c and x0, is also
involved in the procedures of inverse methods.
An iterative process to estimate c and x0 was developed by Hong and Kim [75]. The initial
guess of cohesive zone size for the iterative process was sufficiently large to ensure that the true
cohesive zone is involved in the estimated one. Then the process was iteratively updated by
searching the zero-cross-over points of the traction and separation-gradient distributions within
the estimated cohesive zone. Finally, the iterative process stopped when the traction and
separation-gradient distributions did not make zero-cross-over points any more. In their paper,
the initial guess of x0 was made to be the LEFM effective crack-tip location, which always lay
within the cohesive zone [124]. As stated in Chapter 3 of this dissertation, the LEFM effective
crack-tip location can be obtained by the conservation integrals and the nonlinear least-squares
method.
An alternative way to determine the size c and x-directional centre of the cohesive zone x0 is
the separable nonlinear least-squares method, which is similar to the method described in Section
3.2.1 for the determination of the x-directional and y-directional coordinates of the LEFM
effective crack-tip ( x0 , y0 ) . For the determination of cohesive zone size and positions, the xdirectional centre of the cohesive zone x0 and cohesive-zone size c are treated as two unknown
parameters in addition to the coefficients of eigenfunction expansion. Then, the corresponding
linear least-squares problem is solved iteratively to obtain the coefficients of eigenfunction

83

expansion while updating the size c and the x-directional central position of cohesive zone x0 ,
The iteration keeps going until the sum of residuals converges to a minimum. The flowchart of
separable nonlinear least squares method is shown in Figure 4.3.

Figure 4.3 Flowchart of the separable nonlinear least-squares method for the determination of
cohesive-zone law, cohesive-zone size and position.
84

In the nonlinear separable least-squares method, the nonlinear optimization algorithm used to
set new trial values of x0 and c is Nelder-Mead method, an unconstrained and derivative-free
optimization method [125], [126]. Compared to Newton-Raphson or other Newton-like
algorithms, Nelder-Mead method does not require the computation of the gradient or Hessian of
the objective function, which makes it a more robust optimization method when the solution
space is not convex or when it is hard to obtain the explicit form of the gradient of objective
function [73]. The initial guess of the x-directional central position of the cohesive zone x0 is also
set to be the x-directional LEFM effective crack-tip location, which can be determined by the
conservation integrals and the nonlinear least-squares method.
The initial guess of cohesive-zone size c for the separable nonlinear least-squares method is
relatively flexible. However, a rough estimate is still necessary since the nonlinear optimization
process may become unreliable if the starting position is far away from the true solution. To get a
rough estimate for the cohesive-zone size c, it is assumed that the full-field displacements around
the crack-tip cohesive zone can be represented with only one-term eigenfunction expansion as
follows,
Lf ( z ) 

t0
2





z  c  z  c  i P0

(4.20)

where t 0 and P0 are column vectors with real-number constants. The first-term eigenfunction
expansion can be further expanded into a new series of polynomial functions with respect to
large z ( z  c, c / z  1 ) which represents the elastic far-fields,
Lf ( z ) 

1
t0 2 
c
c
t
z  1   1    iP0  0
2 
z
z
2

85

9
 
  12 1 3  52
2
cz

c
z

O
(
z
)  iP0

8



(4.21)

Dividing the coefficient of term z



5
2

by the coefficient of term z



1
2

, the initial guess of the

cohesive-zone size c can be obtained. Similar to the nonlinear least-squares method in Chapter 3,
the coefficients in the above series of polynomial functions are obtainable with nonlinear leastsquares method, by fitting the analytical fields constructed with the equation (4.21) to the
experimental far-fields.

4.4 Numerical test with synthetic data
To assess the accuracy and stability of two inverse methods, a sequence of numerical tests
were carried out with synthetic cohesive crack-tip fields. Firstly, the synthetic crack-tip fields
were constructed with predefined eigenfunction expansions. Then, both of the field projection
method (FPM) and least-squares method (LSM) were applied to extract the cohesive tractions
and separations from the synthetic displacement fields. Finally, the inverse solutions were
compared with the predefined values. The details about the procedure of numerical tests were
introduced in Section 4.4.1.
In this study, only the mode-I problem was investigated. In other words, only the second
components of the cohesive traction t (x) , separation-gradient b(x) and separation δ(x) were
investigated. However, it is noteworthy that the eigenfunction expansion of the cohesive-cracktip fields developed in Section 4.2 provides a platform for the present work to be extended to
study the mode-II or mixed mode problems.

86

4.4.1 Description of numerical test
4.4.1.1 Cohesive-zone laws and material property
To validate the generality of the inverse methods, the influence of the shape of cohesive-zone
law (CZL) on the accuracy of inverse solutions was studied. Figure 4.4 shows three
representative CZLs with different shapes which were investigated in the numerical test: convex
upward softening law, linear softening law and concave upward softening law. The convex
upwards softening law (CZL I) was used to simulate the fracture behaviors of PMMA [74] or
HIPS [76] , the linear softening law (CZL II) was used to simulate the fracture process of high
explosives [67], while the concave upward softening law (CZL III) was used for some quasibrittle materials such as concrete, rocks or wood [66], [127], [128]. In the test, the CZLs were
normalized to extend the generality of test results. Consequently, widely-applicable guidelines
for experiments can be obtained from the numerical test results.

Normalized Traction

1

I
II
III

0
0

Normalized Separation

1

Figure 4.4 Three representative CZLs used in the numerical test

87

Each cohesive-zone law was constructed with predefined coefficients t n and b n in the

~) vanishes because it does not contribute to the
eigenfunction expansion (4.8). In addition, s(z
cohesive-zone variables. Thus, the eigenfunction expansion used in the numerical test was in the
form of
Lf ( z ) 

N
N
1

z  1 t nU n ( z )  z  1 b nU n ( z )  .

2
n 0
n 0


(4.22)

In this study, only the first three sets of coefficients t n and b n (n  0,1, 2) were predefined to
construct the synthetic displacement fields around crack-tip cohesive zone. The limit was set for
the reason that a three-term representation is accurate enough to depict most of the common
functional curves. Hence the details of CZL curves can be well captured. Table 4.1 shows the
values of coefficients t n and b n (n  0,1, 2) corresponding to each of the CZLs in Figure 4.4.
Here,  1 ,  2 and  3 are the incipient stresses of decohesion for the three CZLs. The relationship
among them is set to be 1 :  2 :  3  3.66 : 4.64 : 6.83 to ensure that the cohesive fracture
energies of the three CZLs are identical.

88

Table 4.1 Coefficients t n and b n (n  0,1, 2) for the construction of different CZLs
Coefficients

Convex upward law (I)

Linear law (II)

Concave upward law(III)

t0

σ1[0, 0.9871]T

σ2 [0, 0.5626]T

σ3 [0, 0.2405]T

b0

σ1[0, 0.6856]T

σ2 [0, 0.7571]T

σ3 [0, 0.6560]T

t1

σ1[0, -0.1756]T

σ2 [0, 0.1450]T

σ3 [0, 0.1664]T

b1

σ1[0, -0.1258]T

σ2 [0, 0.0495]T

σ3 [0, 0.2492]T

t2

σ1[0, 0.0249]T

σ2 [0, -0.0478]T

σ3 [0, 0.0414]T

b2

σ1[0, 0.0249]T

σ2 [0, -0.0478]T

σ3 [0, 0.0414]T

With the coefficients t n and b n (n  0,1, 2) specified, the analytical solutions of cohesive
tractions, separation-gradients, and separations can be obtained using the equations (4.9), (4.10)
and (4.11), while the synthetic cohesive crack-tip displacement fields can be obtained using the
equation (3.3). Besides, the material properties are also required in the generation of analytical
cohesive-zone variables and cohesive-crack-tip fields. The bulk material was assumed to be the
same as the fiber-reinforced composite laminates studied in Chapter 3. The orthotropic elastic
properties of the material are shown in Table 3.1.
Then the two inverse methods, namely field projection method and separable nonlinear leastsquares method, were employed to extract CZLs from the synthetic cohesive crack-tip
displacement fields. Moreover, to directly compare the robustness of two methods, the
interaction J-integral in the field projection method was implemented by domain integral
formulation [113], with the same input data points for separable nonlinear least-squares method.

89

4.4.1.2 Procedure of numerical test
The computational implementations of both field projection method and separable nonlinear
least-squares method can be divided into two sections: (1) the determination of the coefficients
t n and b n (n  0,1, 2) in the eigenfunction expansion (shown in the equation (4.22)) with known

cohesive-zone size and location; (2) the estimate of cohesive-zone size and location. Therefore,
numerical tests need to be separately carried out for each section to accomplish a comprehensive
investigation of the inverse methods. Moreover, for the estimate of cohesive-zone size and
location, both of the two methods introduced in Section 4.3.3, i.e. the method of searching zerocross-over points on the cohesive traction and separation-gradient distributions and separable
nonlinear least-squares method, are iterative method. At each iteration step, the estimation
process of cohesive-zone size and location requires the implementation of field projection
method or linear least-squares method to determine the coefficients t n and b n (n  0,1, 2) with the
estimated cohesive-zone size and location. Therefore, the accuracy of estimated cohesive zone
size and location is governed by the accuracy of extracted coefficients in the eigenfunction
expansion. In other words, the first section of numerical test is the foundation of the second one.
In the first section of numerical test, several factors which may influence the accuracy of
inverse solutions were investigated, including the shape of CZL (as shown in Figure 4.4), the
inverse distance of input data field, the number of the data points and the noise level of input
data.
As shown in Figure 4.5 (a), the inverse distance of the displacement data field is a normalized
parameter, defined as the ratio between the half width of inner deleted area WI and the size of
cohesive zone C . The inverse distance represents the distance between the data field used in the
computation and the crack-tip cohesive zone. The reason why the inverse distance is used and
90

how it is defined are explained as follows. On one hand, it is difficult to experimentally obtain
reliable displacement data close to the cohesive crack-tip because of the large deformation
gradients within this region. Hence, the displacement fields near the crack-tip need to be
eliminated from the inverse computation. On the other hand, the numerical errors in the inverse
computation process increase significantly when the distance between input data fields and
crack-tip cohesive zone increases. This issue was investigated in the original paper of field
projection method [75], which demonstrated that the inverse solutions of cohesive-zone variables
showed instability while using the data at a considerable distance away from the cohesive zone.
Therefore, there exists a maximum limit of the inverse distance, within which the CZL can be
extracted accurately. To determine the limit, the accuracy of inverse solutions according to
different inverse distance should be tested. In the numerical tests, the range of inverse distance
was varied from 2 to 10. Without loss of generality, the value of inverse distance was changed by
only changing C while fixing WI .
While the inverse distance represents one important feature of input data field for inverse
computation, the number of data points is another important feature. The influence of the number
of data points on the accuracy of inverse solutions was investigated. In the test, the number of
data points was changed with the size of data field fixed while changing the interval between
neighbouring data points.
The noise sensitivity is a very important concern in evaluating the robustness of inverse
methods because noise in the input data is inevitable in practice. To extend the generality of
numerical test results, the noise level was measured by a normalized standard deviation (NSTD),
which was defined as the standard deviation (STD) normalized with respect to the difference
between the maximum and minimum of displacement component along loading direction, i.e.

91

NSTD  STD of noise/ max(V)  min(V)  . The order of magnitude of NSTD in the
numerical test was 10 3 , which was chosen due to the fact that it is obtainable in practice. For
instance, if Digital Image Correlation (DIC) method is used, the standard deviation of
displacement measurement can easily reach 3 10 2 pixels or even smaller [71]. In such case,
the NSTD can be made to be smaller than 110 3 by adjusting the resolution of camera or the
input data field to make sure that max(V)  min(V)  30 pixels . In the numerical test, random
noise which was generated from a standard normal distribution with the mean value of noise
being zero was added to the synthetic displacement data. Illustrations of the noise in synthetic
displacement fields are shown in Figure 4.5 (c) and (d). Numerical tests for each noise level were
repeated 20 times with different sets of randomly-generated noise. The average of results was
then used as the representative result.

92

(a)

(b)

(c)

(d)

Figure 4.5 Characteristic parameters of input data fields and displacement measurement noise for
the convex upward CZL (Law I): (a) Horizontal analytical displacement field U (perpendicular
to the loading direction) (b) Vertical analytical displacement field V (along the loading direction)
(c) Horizontal analytical displacement field U with the NSTD of noise being 110 3 (d) Vertical
analytical displacement field V with the NSTD of noise being 110 3 .
After the inverse extraction of coefficients t n and b n (n  0,1, 2) in eigenfunction expansion, a
(n)
relative error defined as t2( n )  t20

(n)
t20
was used to investigate the accuracy of inverse solutions.

(n )
Here t2( n ) denotes the second component of t n , t20
represents the exact value. Besides, relative

93

errors were computed through comparisons between the exact values and the inverse solutions of
non-dimensionalized



~
traction t ( ~
x )  t( ~
x) /

,

non-dimensionalized

separation-gradient



~
b( ~
x )  H 1 /  b( ~
x ) (  is the incipient stress of decohesion which is used in the construction of
~
synthetic fields) and non-dimensionalized separation δ( ~
x ) denotes
x )  δ( ~
x ) / max( δ0 ( ~
x )) ( δ0 ( ~

the exact value of separation). Moreover, the normalized root-mean-square (NRMS) error was
used to quantify the difference between the inversely-extracted CZL and its exact value,



~
max(
~ ~0
NRMS  (1 / max( max ,  max
))0

~0
max , max )

~t (~)  ~t (~ ) d~
0

0

2

1/ 2

(4.23)

~
~ ~
where ~
t ( ) represents the extracted non-dimentionalized CZL, while t 0 ( 0 ) denotes the exact

value.
As stated previously, the second section of numerical test was carried out based on the result
of the first section. Because the iterative computational procedure determines that it is possible to
obtain accurate cohesive-zone size and position only when the coefficients in the eigenfunction
expansion can be extracted with acceptable accuracy at each iteration step. The same factors
which can influence the accuracy of the extracted coefficients were also investigated in the
second section, including the shape of the CZL, the inverse distance of the data field, the noise of
data field, etc. The ranges of tested factors in the second section were the subsets of those in the
first section. In addition, to verify the robustness of iterative methods for the determination of
cohesive-zone size and position, different initial guesses were tested.

94

4.4.2 Results of test
4.4.2.1 Determination of CZLs with known cohesive-zone size and position
In this section, numerical tests were conducted to determine the coefficients t n and b n
(n  0,1, 2) and hence the CZLs with known cohesive-zone size and location. The results are

presented as follows.
First of all, a baseline test was conducted to verify the theoretical correctness of two inverse
methods—the field projection method (FPM) and the linear least-squares method (LSM). The
synthetic displacement data without noise was used. The two characteristic lengths used to define
the range of data field were the inverse distance and the half width of outer square, as shown in
Figure 4.5 (a). In this baseline test, WI / C  2 , and WO / C  6 . The total number of input data
points was 18816. Figure 4.6 (a)-(d) respectively shows the comparison between the inverse
solutions and the predefined values (i.e. exact values) of tractions, separation-gradients,
separations and CZLs. It can be seen that the inverse solutions of 3 different CZLs all converge
to the exact values, which verifies the theoretical derivations and computational implementations
of two inverse methods.

95

1.5

exact
FPM
LSM

1

b2

t2

1

1.5

exact
FPM
LSM

0.5

0
-1

0.5

-0.5

0
x/c
(a)

0.5

0
-1

1

-0.5

0
x/c
(b)

0.5

1.5

1.5

exact
FPM
LSM

exact
FPM
LSM

1

2

t2

1

1

0.5

0.5

0
-1

-0.5

0
x/c
(c)

0.5

0
0

1

0.5

2

1

1.5

(d)

Figure 4.6 Comparison between the inverse solutions and exact values of the cohesive-zone
variables and cohesive-zone laws: (a) traction (b) separation gradient (c) separation (d) cohesivezone law.

Starting from the baseline test, a sequence of numerical tests were carried out to study the
influence of inverse distance, number of the data points and noise on the accuracy of inverse
solutions. Due to limited space, only the results of linear softening law (CZL II) are shown here.
The results of concave upward softening law (CZL III) and convex upward softening law (CZL
I) are analogous to that of linear softening law.
Firstly, the influence of inverse distance WI / C on the inverse solutions was investigated. As
mentioned above, the value of WI / C was changed from 2 to 10, by changing the cohesive-zone

96

size C while fixing the half width of inner square WI . The other factors which define the input
data fields are the same as that of the baseline test. Figure 4.7 (a) and (b) respectively exhibits the
CZLs extracted by the field projection method (FPM) and the linear least-squares method
(LSM), with respect to different inverse distances. It can be seen that the errors of CZLs
extracted by both methods increase when the inverse distance WI / C increases. For FPM, it has
been found that the increasing error of CZLs with respect to the inverse distance was due to the
singular behaviour of the inversion scheme, which was caused by the term 2 N (r / c) N1 / 2 in the
integrand of interaction J-integral [75]. For LSM, the increase of inverse distance makes the
matrix in the normal equation more ill-conditioned, so that the computational error becomes
larger.

1.2

1.2

exact
W /C=3

exact
W /C=3

I

1

I

1

W /C=6

W /C=6

W /C=10

W /C=10

I

I

t2

t2

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.4



0.6

0.8

I

0.8

I

0.8

0

1

2

(a)

0

0.2

0.4



0.6

0.8

1

2

(b)

Figure 4.7 Influence of inverse distance on the accuracy of inverse solutions: (a) the CZLs
extracted by the field projection method (FPM) (b) the CZLs extracted by the linear least-squares
method (LSM).

97

Figure 4.8 shows the relative errors of extracted coefficients corresponding to different-order
terms in the eigenfunction expansion. The relative errors are plotted with respect to the inverse
distances. The term t (2n ) in the legend represents the second component of eigenfunction
expansion coefficients t n (n  0,1, 2) (the first component is zero due to the mode-I problem).
The result of coefficients b n is similar to that of t n and hence is not presented. It can be seen that
the errors of higher-order-term coefficients are always larger than the errors of low-order-term
coefficients. In addition, the variation trends of errors with respect to the inverse distance for
both inverse methods are different. For FPM, when inverse distance WI / C increases, the relative
errors of the coefficients in higher-order terms increase faster than the coefficients in lower-order
terms. For LSM, the relative errors of all the coefficients increase dramatically at the first stage,
and reach plateaus when WI / C is larger than 6.

Figure 4.8 Relative errors of the eigenfunction expansion coefficients with respect to the inverse
distance.
98

In addition, when WI / C  10 , the relative errors of the 3rd-term coefficients t (2)
for both
2
inverse methods are around 100 percents. From Figure 4.7, it can be seen that the extracted CZLs
at this moment can barely represent the shapes of CZLs. So the value 100 percent in the relative
error of the 3rd-term coefficients t (2)
is set as a limit to judge if the accuracy of extracted CZL is
2
acceptable or not. In other words, the inverse solution of CZL is considered as a reliable result
only when the relative error of the 3rd term is smaller than 100 percent.
Another important observation which can be obtained from Figure 4.8 is the comparison of
accuracy between two methods. The relative errors in the coefficients extracted by FPM are
smaller than those obtained by FPM when the inverse distance is smaller than 4. In contrast, the
FPM can get more accurate result when the inverse distance increases from 5 to 9. These
observations are helpful for the selection of appropriate method, based on the inverse distance of
input data.
In all the above analyses, the number of data points was arbitrarily chosen. However, the
obtainable number of data points in practice relies on the experimental measurement techniques.
In this chapter, it is assumed that the displacement data used for the inverse computation is
measured by Digital Image Correlation (DIC). Thus, the domain of displacement data and the
interval between neighbouring data points were quantified in pixels. In all the numerical tests,
the data domain was fixed with the sizes WO  360 pixels and WI  120 pixels , which are
obtainable in practice.
The influence of the number of data points on the inverse solutions was investigated. In this
test, the inverse distance was set to be the same as the baseline test, WI / C  2 and WO / C  6 .
The number of data points was modified by changing the interval between neighbouring points.

99

For the baseline test, the interval between neighbouring data points was 5 pixels. In the new
tests, different input data sets were generated by setting the intervals between neighbouring data
points to be 5 pixels, 10 pixels, 15 pixels, 20 pixels and 30 pixels. The relative errors of the
extracted eigenfunction expansion coefficients t (2n ) are shown in Figure 4.9.

Figure 4.9 Relative errors of the eigenfunction expansion coefficients extracted from various
data sets with different intervals between neighbouring data points. WI / C  2 and WO / C  6 .

For LSM, it can be seen that the relative errors of the extracted coefficients t (2n ) decrease
slightly as the interval between neighbouring data points increases (i.e. the number of data points
decreases). It is due to the fact that the matrix of normal equation in LSM becomes less illconditioned with less number of input data points. Contrarily, the increase of interval between
neighbouring data points reduces the accuracy of FPM results. The relative error of the 3rd-

100

coefficient t (2)
becomes larger than 100 percent when the interval increases to 30 pixels. The
2
low accuracy of the FPM results under large intervals is attributed to the large numerical errors
in the calculated displacement gradients, which is required in the computation of interaction Jintegral. Specifically, the raw input data from DIC is the displacement. The displacement
gradients are determined by numerical differentiation of the displacements. For the test results
shown in Figure 4.9, there is no noise in the input data. Thus, the source of error in the calculated
displacement gradients is numerical truncation error, which can be illustrated as follows.
Assume that the numerical differentiation algorithm used to estimate the displacement
gradient is central difference method, the truncation error ET can be represented as:
ET 

u '
u ( x  x)  u ( x  x) '
 u ( x) 
 u ( x)
x
2x

where u represents the displacement,

x represents

(4.24)

the position or coordinate of data point. From

Taylor series expansion, the equation can be re-written as:.
1
ET  u ''' ( x)(x)2  O((x)4 )
6

(4.25)

It can be seen that the truncation error will increase when the interval between neighbouring data
points x increases.
There is no noise in the synthetic input data for all the previous numerical tests. However,
noise is inevitable in the practical data sets which are measured by experimental techniques. It
may significantly reduce the accuracy of inverse solutions for both methods. Take FPM as an
example, if noisy input data is used, there will be an additional term in the total error of
numerical differentiation of displacement. Denoting the error that is caused by the noise of input
data as EN , the total error of numerical differentiation can be represented as

101

u*
u* u
u
 1 '''
'
E  EN  ET 
 u ( x)  (
 )  (  u ' ( x)) 
 u ( x)(x)2  O((x) 4 ) (4.26)
x
x x
x
x 6

where u* ( x)  u(x)   ( x) , and  denotes the noise, i.e. the measurement error. It can be
concluded that the error caused by noise will increase as the interval between neighbouring data
points x decreases. It is opposite to the variation trend of numerical truncation error.
Furthermore, the equation (4.26) suggests that there exists an optimal interval x , with which the
displacement gradients and hence the inverse solutions of FPM calculated from noisy data will
be the most accurate.
A set of numerical tests were then conducted to investigate the influence of noise on the
accuracy of inverse solutions. In the tests, the inverse distance of the input data was WI / C  2 ,
as same as the baseline test. Various data sets were used as inputs, with the intervals between
neighboring data point equaling to 5, 10, 15 and 20 pixels. The maximal interval for the test, 20
pixels, was set according to the results in Figure 4.9, which shows that it is possible to get
inverse solutions with acceptable accuracy (i.e. the relative error of the 3rd-coefficient t (2)
is
2
smaller than 100 percent) only when the interval is smaller than 20 pixels, In addition, the
numerical tests for each noise level were repeated 20 times. Then the medians of all the 20 sets
of test results were used as the representative result. Limited by the space, only the results with
interval equaling to 5 pixels and 20 pixels were presented. Figure 4.10 (a)-(b) show the relative
errors of the eigenfunction expansion coefficients extracted from various data sets with different
noise levels.

102

(a)

(b)
Figure 4.10 Relative errors of the eigenfunction expansion coefficients extracted from various
data sets with different noise levels: (a) interval between two neighbouring points is 5 pixels (b)
interval between two neighbouring points is 20 pixels.
103

Figure 4.10 reveals several important features regarding the accuracy of inverse solutions.
Firstly, significant increases in the errors of the extracted coefficients can be observed for both
inverse methods when there is noise in the input data. When the noise level is 0.5e-3 (NSTD),
only barely visible effect can be seen in the displacement field. However, errors in the extracted
coefficients are several orders of magnitude larger than the errors obtained with no-noise inputs.
In addition, it can be seen that the increase of errors becomes slow after the first leap between the
test with no-noise input data and that with noisy input data.
As for the comparison of two methods, it can be seen that the errors of coefficients obtained
by LSM are always smaller than those obtained by FPM, with the same noise level and interval
in the input data. A more important observation is that the relative errors of the 3rd-term
coefficient extracted by FPM are always larger than 100 percent, even with the smallest noise
level, NSTD 0.5e-3, which generates nearly invisible fluctuations in the input data fields. It is
highly possible that the noise level in the experimental data is larger than 0.5e-3 (NSTD).
Therefore, it can be concluded that the FPM is not reliable to extract CZLs which requires at
least three terms in the eigenfunction expansion to depict their functional shapes. On the other
hand, the largest error of the 3rd-term coefficient extracted by LSM is always less than 50 percent
with all the noise levels. Therefore, the LSM is preferable for practical use if the input data has
considerable noise level.
In all the previous numerical tests, the cohesive-zone size and central position were assumed
to be known a priori, for the purpose of isolating the sources of error in the inverse solutions.
Generally, it is expected that the error in the final inverse solutions of CZL will become larger if
the determination of cohesive-zone size and position is also involved in the inverse scheme.

104

Therefore, it can be concluded that the FPM will never get inverse solutions that is accurate
enough to depict the details of CZLs (i.e., the 3rd term coefficient should be smaller than 100
percent) if noisy inputs are used. So the LSM is the only method that will be studied in the
following numerical tests, which involves the accuracy estimate of the simultaneously
determined cohesive-zone size, central position and CZL.

4.4.2.2 Determination of CZLs with unknown cohesive-zone size and position
In this section, the inverse extraction of CZL was implemented without pre-known cohesivezone size and position. The numerical tests in Section 4.4.2.1 were the bases for the design of
numerical tests in this section. Firstly, the FPM was eliminated from the tests because of its high
noise sensitivity. In addition, the same factors which can influence the accuracy of inverse
solutions in Section 4.4.2.1 were also investigated in this section, including the shape of the
CZL, the inverse distance of the data field, the noise in the data field, etc. The tested ranges of
these factors were the subsets of the ranges in the previous tests.
In Section 4.3.3, two methods were introduced for the estimate of cohesive-zone size and
central position (i.e., c and x0 ). In this section, the applicability of two methods were examined
with a set of baseline tests. Firstly, tests were conducted for the method developed by Hong and
Kim [75], which determines c and x0 by iteratively searching for the zero-cross-over points of the
traction and separation-gradient distributions within the estimated cohesive zone. The iterative
process is supposed to stop when the estimated traction and separation-gradient distributions do
not make zero-cross-over points any more. The estimates of c and x0 at the stopping moment are
the final solutions. In the test, the traction and separation-gradient distributions within the
estimated cohesive zone at each iteration step were obtained using linear least-squares method.
105

Although the theoretical design of this method seems reasonable, it was found that this method is
not robust in the numerical implementation. In most of the test cases, the solutions of c and x0
were different from the exact values, even with un-noisy synthetic data. In addition, different
initial guesses of c and x0 were tested, but could not improve the accuracy of solutions.
Moreover, it has been found that the decrease of inverse distance in the input data could reduce
the error in the solutions. Therefore, the inaccurate solutions of c and x0 should be due to the
errors in the extracted eigenfunction expansion coefficients during the iterative process.
Nevertheless, the overall performance of this method in the numerical tests did not demonstrate
its applicability, especially for the data with considerable noise level and inverse distance.
Numerical tests were then conducted to demonstrate the applicability of separable nonlinear
least-squares method, whose flowchart is shown in Figure 4.3. Essentially, the inverse
computation procedure of separable nonlinear least-squares method can be divided into two
coupled parts. One is the nonlinear optimization process for the estimate of cohesive-zone size
and position, while the other is the linear least-squares method for the extraction of eigenfunction
expansion coefficients. The accuracy of linear least-squares method with exact cohesive-zone
size and position have been investigated in the previous section. Therefore, the tests in this
section focused on the accuracy of cohesive-zone size and position estimated with the nonlinear
optimization process.
Firstly, tests with no-noise input data were conducted as the baseline. Here, the inverse
distance was chosen to be WI / C  2 , and the interval between neighbouring points was 10
pixels. The initial guess of cohesive-zone size and position was obtained from the method
introduced in the equations (4.20) and (4.21). Figure 4.8(a)-(c) show the searching paths of
separable nonlinear least-squares method on the contour of objective function. The x and y axis
106

respectively represents the cohesive-zone central position and cohesive-zone size which are
normalized with respect to the exact value. The objective function is defined as the square of
difference between the analytical displacement fields and the displacement fields regenerated
from the inverse solutions of eigenfunction expansion coefficients that are extracted with specific
c and x0 .Therefore, the minimum of objective function corresponds to the accurate inverse
solutions of eigenfunction expansion coefficients as well as the cohesive zone size and position.
Obviously, the minimum is at the origin.

107

(a)

(b)

(c)
Figure 4.11 Searching paths of the separable nonlinear least-square method on the contour of
objective function: (a) for convex upward CZL (Law I) (b) for linear CZL (Law II) (c) for
concave upward CZL (Law III). Inverse distance WI / C  2 .

First of all, high computation efficiency was observed in the test. Usually, the iterative process
took less than 10 steps to get the convergent solution. Secondly, it can be seen that the nonlinear
optimization processes for different CZLs all successfully converged to the correct solution. It
demonstrates the applicability of the nonlinear separable least-squares method, as well as the
108

method for the determination of initial guess. However, a nonlinear optimization method cannot
be considered as a robust method if it can only get correct solution from one specific initial guess.
In addition, it can be expected that the method for the determination of initial guess may not
work for the data with considerable noise level. Therefore, in the following tests, influence of
different initial guesses on the accuracy of inverse solutions were also studied, while using the
initial guess obtained from the equations (4.20) and (4.21) as a reference.
Next, the influence of noise and inverse distance on the accuracy of estimated CZL, cohesivezone size and position were investigated. The interval between neighbouring data points was set
to be 10 pixels. The numerical tests started from small inverse distance (WI / C  1.5) and small
noise level (NSTD of noise is 0.5e-3). Then, either the inverse distance or the noise level was
increased with the increasing interval 0.5 ( WI / C ) or 0.5e-3 (NSTD of noise), respectively. One
parameter is changed while keeping the other parameter fixed. For each CZL, numerical tests
were repeated 20 times for each combination of inverse distance and noise level. The initial
guess of cohesive-zone size and position was obtained from the method introduced in the
equations (4.20) and (4.21). The test stopped when the inverse distance or the noise level was
large enough to generate significant differences between the extracted CZLs and the exact values.
Figure 4.12 summarizes the estimated cohesive-zone sizes and positions, as well as the
extracted CZLs for different noise levels and different types of CZLs, with the inverse distance
WI / C =2 . In the first row, the estimated cohesive-zone sizes and positions are marked on the

contours of objective function. For each CZL, the green line with squares represents the
searching path of nonlinear optimization scheme when there is no noise; the red star marks (total
20 marks) represent the convergent location for each optimization process when the NSTD of
noise is 0.5e-3; the yellow circles represent the results when the NSTD of noise is 1e-3. The

109

corresponding inverse solutions of CZLs are also presented. Each set of results include 20
extracted CZL curves, which are represented by 20 curves with different colors and types.
Cohesive zone law I

Cohesive zone law II

Cohesive zone law III

Estimated
cohesivezone
sizes and
positions

Without
noise

NSTD of
noise
0.5e-3

NSTD of
noise
1e-3

Figure 4.12 Estimated cohesive-zone sizes and positions, and extracted CZLs with different
noise levels. The inverse distance WI / C  2 .
110

In addition, it is worth mentioning that the influence of initial guess of cohesive-zone size and
position on the estimated values was also investigated. Two different initial guesses have been
used. One initial guess was obtained from the method introduced in the equations (4.20) and
(4.21), while the other one was the origin of objective function contour which represents the
exact values of cohesive-zone size and position. In most cases, it has been found that the
nonlinear optimization process converged to the same solution, no matter which initial guess was
used. This observation demonstrated the applicability and robustness of the nonlinear
optimization searching process.
Several important findings were obtained from Figure 4.12. Firstly, it is found that the
accuracy of estimated cohesive-zone size and position is highly sensitive to noise. When there is
no noise in the input data, the determination of cohesive-zone size and position is very accurate,
which leads to an accurate inverse solution of CZL. However, the estimated cohesive-zone sizes
and positions become very scattered on the contours of objective function, even when the NSTD
of noise is as small as 0.5e-3. Furthermore, the error of extracted CZLs is found to be governed
by the error of estimated cohesive-zone size and position. The inverse solutions of CZLs
becomes more scattered when the noise level increases, which is due to the increasing error in
the determination of cohesive-zone size and position.
In this test, the normalized root-mean-square (NRMS) error was used to provide a quantified
measure for the accuracy of extracted CZLs. The expression of NRMS error was given in the
equation (4.23). As stated above, numerical tests were repeated 20 times for each combination of
noise level, inverse distance and type of CZL. Then the average (arithmetic mean) and median of
NRMS errors for the total 20 tests was used as the representative result. In addition, through a

111

visual inspection of extracted CZLs, the 5% of NRMS errors was set as the limit to judge if the
accuracy of extracted CZLs is acceptable or not. Because the shape of CZL can be barely
maintained with this amount of error. The results are summarized in Table 4.2.

Table 4.2 Normalized root-mean-square (NRMS) errors of the extracted CZLs for different
inverse distances and noise levels.
Inverse

Number of failed tests
Average

Distance

of

NRMS

NSTD of noise

Median of NRMS error (NRMS error>0.05)
error

WI/C

out of 20 tests
CZL I CZL II CZL III CZL I

CZL II CZL III CZL I CZL II CZL III

0.5e-3

0.0114 0.0063 0.0047 0.0095 0.0045 0.0031

0

0

0

1e-3

0.0259 0.0145 0.0082 0.0229 0.0069 0.0057

3

1

0

0.0428 0.0286 0.0133

0.037 0.0155 0.0076

5

6

0

2e-3

0.0513 0.0363 0.0298 0.0427 0.0145 0.0175

6

4

4

0.5e-3

0.0315 0.0185 0.0129 0.0338 0.0087 0.0059

3

2

1

0.0646 0.0343 0.0397 0.0474 0.0243

0.017

8

4

5

1.5
1.5e-3

2
1e-3
2.5

0.5e-3

0.0609 0.023 0.0268 0.0405 0.0098 0.0178

6

3

3

3

0.5e-3

0.1055 0.0758 0.0357 0.0531 0.0412 0.0270

12

9

4

Generally, it can be seen that the errors of extracted CZLs increase following the increase of
inverse distance and noise level. In other words, the combination of inverse distance and noise
level determines the accuracy of inverse solution. Another important testing item is the
difference in the accuracy of inverse solutions with different types of CZLs. The influence of
noise level on the accuracy of inverse solutions is found to be highly dependent on the types of
112

CZLs. With the same noise level and inverse distance, the error of inverse solutions for CZL I is
always the largest, while the inverse solutions of CZL III is the most accurate. Correlating this
observation to the estimated cohesive-zone sizes and positions shown in Figure 4.12, it is found
that the accuracy of estimated cohesive-zone size and position for different types of CZL has
different amount of influence on the accuracy of extracted CZLs. For instance, with the same
inverse distance and noise level, even though the error of estimated cohesive-zone size and
position for CZL III is the largest among the three CZLs, the error in the extracted CZL for CZL
III is the smallest.
The results presented in Figure 4.9 and Table 4.2 can provide guidelines for experimental
studies. Indeed, this is the most important goal of numerical tests.

4.4.3 Guidelines for experiments
From Figure 4.9 and Table 4.2, it is realized that the combination of the inverse distance and the
noise level of input data governs the accuracy of inverse solutions. The inverse distance and
noise should be modified accordingly to ensure that inverse solutions with acceptable accuracy
can be obtained. For instance, when the inverse distance WI / C is 1.5, it is possible to extract the
CZLs only when the NSTD of noise is equal to or smaller than 2e-3. On the other hand, when
the NSTD of noise is 2e-3, the inverse distance WI / C should be less than 1.5 to ensure that the
accuracy of extract CZLs is acceptable.
Moreover, since the shape of CZL usually cannot be known a priori in practice, the safest
combination of inverse distance and noise level should be followed as the guideline. From Figure
4.9 and Table 4.2, it can be realized that the CZL I is usually the most difficult one to achieve

113

accurate inverse solutions. Thus the setup of inverse distance and noise level need to be modified
based on the results of CZL I. Specifically, the inverse distance WI / C should be respectively set
to be less than 2.5, 2 and 1.5, when the NSTD of noise are 0.5e-3, 1e-3 and 2e-3, and vice versa.
Some additional guidelines are worth mentioning here if the displacement data are measured
by Digital Image Correlation (DIC). Recall that the parameter used to measure the noise level
NSTD is defined as NSTD  STD of noise /max( V )  min(V ). Since the measurement accuracy
of DIC, i.e. the STD of noise, is fixed and measured in pixels, the NSTD can be modified by
changing the range of displacements which is easily achieved by modifying the resolution of
captured images (physical length/pixel). However, it is noteworthy that the inverse distance will
be changed automatically while changing the resolution of image. Therefore, several factors need
to be considered simultaneously for the experimental setup.
At last, it should be pointed out that an important simplification was made in the numerical
tests, which is different from reality. For the numerical analysis, it was assumed that the elastic
fields around cohesive zone can be represented with a 3-term eigenfunction expansion. It is a
reasonable assumption when the region of interest is close to the crack-tip cohesive zone.
However, the higher-order terms in the eigenfunction expansion may contribute considerably to
the representation of cohesive crack-tip fields some distance away from the cohesive zone.
Moreover, it was found that the difference between the elastic fields defined with the cohesivezone eigenfunction expansion and the LEFM Williams expansion with identical J-integral values
become very small at the elastic far-fields. Thus, the accuracy and uniqueness of inverse solution
could be reduced if using the data fields far away from the cohesive zone in the inverse
computation. A data field that is as close as possible to the crack-tip cohesive-zone is preferable
in getting accurate inverse solutions.

114

4.5 Summary
In this chapter, two analytical inverse methods were proposed for the extraction of cohesive
fracture properties from full-field displacement data around the crack-tip, namely the field
projection method (FPM) and the separable nonlinear least-squares method. Both methods were
developed based on an eigenfunction expansion of the elastic fields around a crack-tip cohesive
zone. The coefficients of eigenfunction expansion, the cohesive-zone size and position can be
determined with the two inverse methods. Then the traction and separation distributions within
the cohesive zone can obtained with the extracted eigenfunction expansion coefficients, the
cohesive-zone size and position. Finally the cohesive-zone law can be generated as the
relationship between cohesive-zone traction and separation.
A set of numerical tests were carried out to investigate the applicability and robustness of the
two inverse methods. Several factors that can influence the accuracy of inverse solutions were
investigated. The tested factors involved the noise of input data, the number of data points, the
inverse distance of data field, and the shape of the cohesive-zone law. It has been found that the
field projection method had high noise sensitivity, and hence had limited applicability in
practice. On the other hand, the separable nonlinear least-squares method, which consists of the
extraction of eigenfunction expansion coefficients in a linear least-squares sense and the
nonlinear iterative searching process of cohesive-zone size and position, has been proven to be a
robust inverse method. It was reliable to extract accurate CZLs with practically obtainable noise
levels and inverse distances. Therefore, this method is recommended for experimental studies.
Guidelines for the experiment setup and data processing were provided based on the results of
numerical test. It was found that the noise level and the inverse distance worked together to
115

determine the accuracy of inverse solution. The detailed guidelines can be found in Table 4.2.
As a sequel of the analytical and numerical studies presented in this chapter, experimental
validation with DIC will be presented in Chapter 5.

116

Chapter 5
Inverse Extraction of Cohesive-zone Laws for Fiberreinforced Composite Laminates Part II: Experimental
Validation
5.1 Introduction
In Chapter 4, two analytical inverse methods, namely field projection method and separable
nonlinear least-squares method, have been developed to measure the cohesive-zone law (CZL) of
fiber-reinforced composite laminates. The methods were developed based on the eigenfunction
expansion of cohesive crack-tip fields for anisotropic solids. Numerical tests were conducted to
assess the applicability and robustness of the two methods. It was found that the separable
nonlinear least-squares method had the potential to be used in practice. Some guidelines for
experimental setup and data processing were provided.
In this chapter, an experimental validation of the separable nonlinear least-squares method is
presented. The inverse method is applied to measure the CZL for the translaminar fracture of a
laminated composite. In the experiment, the crack-tip displacement fields measured by 2-D
digital image correlation (DIC) are used for the inverse extraction of CZLs. The experimental
setups are based on the guidelines provided in Chapter 4.
The remainder of the chapter is organized as follows. In Section 5.2, the experimental setups
and data reduction scheme are introduced in details. Section 5.3 shows the experimental results,
117

including the experimental displacement fields, inverse solutions of CZLs, cohesive-zone size
and position. The accuracy of inverse solutions are discussed in Section 5.4, through the
computation of fracture energies and the finite element simulation with ABAQUS. Finally, the
research achievements are summarized in Section 5.5.

5.2 Experimental procedure
5.2.1 Experimental setups

The composite material and specimen geometry in the experimental study presented in
Chapter 3 are also used in this chapter for the experimental validation of inverse method. The
material properties are provided in Table 3.1. The specimen configurations are shown in Figure
3.7 (a).
The CZL of the fiber-reinforced laminated composite was measured based on a mode I
translaminar fracture test which was conducted using the extended compact tension (ECT)
specimen geometry [118] , following ASTM E1922 [119]. The ECT specimen can exhibit stable
crack growth, which is critical in the inverse extraction of CZLs. The experimental
displacements measured in Chapter 3 for the estimate of LEFM crack-tip field parameters and Rcurves were also used in this chapter for the inverse extraction of CZLs. However, only those
displacement measurements during the steady-state crack propagation were chosen for the
measurement of cohesive-zone laws. The reason will be explained in Section 5.2.2.

118

5.2.2 Data reduction scheme
The full-field crack-tip displacements measured by DIC which had been used in Chapter 3
were also used in this chapter. Instead of using the whole data from the beginning of fracture test,
it is noteworthy that only the displacements measured during the steady-state crack growth can
be used for the inverse extraction of cohesive-zone laws. It is due to the fact that the analytical
inverse method was developed based on an eigenfunction expansion of crack-tip fields
surrounding a complete cohesive zone. Essentially, the cohesive zone is completely generated at
the start of steady-state crack growth, and remains complete during the steady-state crack growth.
Figure 5.1 shows the evolutions of cohesive traction and separation distributions (solid curves)
and the corresponding equivalent LEFM crack-tip opening profile and stress distribution (dashed
curves) at different stages during the cohesive fracture process. In the first three figures, the red
curves represent the separation and crack-tip opening profile δ(x), while the blue curves
represent the cohesive traction and crack-tip stress distribution σ(x). The solid triangle and
hollow triangle in each figure denotes the equivalent LEFM crack-tip and the tip of cohesivezone, respectively. δcr and σcr denote the critical cohesive traction and separation, respectively.
The stages of fracture process corresponding to the first three figures are marked as A, B and C
with solid triangles on the LEFM R-curve plot. It can be seen that the cohesive zone is firstly
growing (Stage A) to the complete size (Stage B) when the R-curve reaches its plateau, i.e. at the
beginning of steady-state crack propagation, and thereafter propagates in a self-similar steadystate way (Stage C). The steady-state propagation of cohesive zone requires that the size of the
cohesive zone and the cohesive traction and separation distributions over the cohesive zone
remain constant [124]. Besides, another consequence of the steady-state regime is the

119

equivalence between the plateau value of R-curve, i.e. GRc, and the value of cohesive fracture
energy which is equivalent to the area under the curve of cohesive-zone law.

Figure 5.1 Sketch of the cohesive zone behavior during fracture process: the cohesive traction
and separation distribution profiles, the corresponding equivalent LEFM crack-tip opening
profile and stress distribution; and the different stages of crack growth.

120

The accuracy of extracted CZLs can be verified by examining the maximum traction and
cohesive fracture energy, which should be equal to the uniaxial tensile strength and the LEFM
energy release rate during steady-state crack propagation, respectively. Another verification
method is to compare the experimental global response with the results of finite element
simulation which uses the extracted CZLs. Since the CZLs are inversely extracted from the
experimental displacements, the "direct" finite element simulation with the CZLs should be able
to duplicate the experimental outcomes.

5.3 Experimental results
5.3.1 Displacement fields
The full-field displacements which had been obtained in Chapter 3 for the measurement of
LEFM R-curve were also used in this chapter for the measurement of cohesive-zone laws. As
stated above, only the data during steady-state crack growth were used.
Figure 5.2 and Figure 5.3 shows the Load-CMOD plot and the R-curves in terms of stress
intensity factor, respectively. The black dots on the curves in the two figures represent the
moments when the displacement data are chosen. It can be seen that there are total 11 data sets
used for the inverse computation. The corresponding time of data acquisition is from 90th second
to 100th second after the start of experiment. These moments are at the decreasing stage of LoadCMOD curve, in the meantime the crack is propagating in a steady-state manner.
In Figure 5.4, the measured displacement fields (U and V) at three moments (I, II and III)
during the steady-state crack propagation are presented. The evolution of displacement fields
121

shows a pure mode I stable crack growth. The choice of input data field for the inverse
computation follows the guideline provided by the numerical test results in Chapter 4. Generally,
for a fixed noise level, the inverse solutions were more accurate with smaller inverse distance. In
other words, the inclusion of data close to the crack-tip cohesive zone can help reduce the error
of inverse solution. Through an examination of the continuity and smoothness of experimental
displacements, it was found that the region within which the displacement measurement was
unreliable was very close to the crack free surface and crack tip. Therefore, a rectangular region
near the crack free surface and crack tip was eliminated from the input data fields, as shown in
Figure 5.4. In addition, a square data field is used, with the outer half width of the data field
equaling 360 pixels, the same as the numerical tests. More than 3000 data points were used in the
inverse computation.

122

1400
1200

Load (N)

1000

I
II

800

III

600
400
200
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

CMOD (mm)

Figure 5.2 Load versus Crack Mouth Opening Displacement (CMOD)

Figure 5.3 R-curves in terms of stress intensity factor, determined by the nonlinear least-squares
method
123

(a)

(b)

(c)
Figure 5.4 Displacement fields measured at different moments shown in Figure 5.2 and 5.3
during the steady-state crack growth: (a) at the moment I (b) at the moment II (c) at the moment
III. V is the displacement component along the loading direction.

124

5.3.2 Cohesive-zone variables and cohesive-zone laws
The experimental results of cohesive-zone properties are shown in this section. First of all, the
inverse solutions for one displacement data set are presented. The data set was obtained at the
moment II, as shown in Figure 5.2 and Figure 5.3. The corresponding displacement fields at this
moment are presented in Figure 5.4 (b). For the inverse computation, the LEFM equivalent
crack-tip was used as the initial guess of cohesive-zone central position. On the other hand,
several different initial guesses of cohesive-zone sizes were used. The consistency of inverse
solutions with different initial guesses was checked to examine if the nonlinear optimization
process successfully converged to the global minimum of objective function. The extracted
cohesive-zone traction t, separation gradient b, separation δ, and CZL are shown in Figure 5.5.
The estimates of cohesive-zone sizes and positions are shown in Figure 5.6. In both figures, the
cohesive-zone sizes and positions are presented in the unit of pixels, while the resolution of
image is 22.5 m/pixel.
From Figure 5.5 and Figure 5.6, it is found that there are only two different sets of inverse
solutions for total seven different initial guesses. Through a visual inspection of the extracted
cohesive-zone variables, it can be easily concluded that the set of inverse solutions
corresponding to a longer cohesive zone is incorrect. For instance, it is noticeable that the
extracted separation distribution vanishes in the middle of cohesive zone, which contradicts the
basis of CZM. The other set of inverse solution which has a shorter cohesive zone is more
reasonable. A further verification for the accuracy of inverse solutions will be presented in
Section 5.4.

125

Figure 5.5 Cohesive-zone traction t, separation gradient b, separation δ, and CZL (relationship
between the cohesive-zone traction and separation) extracted from the data set obtained at the
moment II. Different initial estimations of cohesive-zone size were used in the inverse
computations.
1000

Extracted cohesive zone size and location
LEFM effective crack-tip location (845.7 pixels)

Cohesive-zone size and location (pixel)

980
960
940

X =899.4 pixels
0

C=74.0 pixels

920
X =845.1 pixels

900

0

C=35.6 pixels

880
860
840
820
800
40

50

60

70
80
90
100
Initial estimation of C (pixel)

110

120

Figure 5.6 Cohesive-zone sizes and positions extracted from the data set obtained at the moment
II. Different initial estimations of cohesive-zone size were used in the inverse computations.
126

Following the above procedure, inverse computation was conducted for each input data set
Different initial guesses were used to ensure an accurate inverse solution. However, reasonable
results cannot be obtained for three cases among the total eleven cases, no matter what initial
guesses were used. The results of the other eight cases which were measured at different time are
presented in Figure 5.7 and Figure 5.8. The average of CZLs is also presented, which will be
used as the representative result in the following studies. Some of the CZLs in Figure 5.7
originally had small negative tails. These artificial features were trimmed away since they
violated the basis of CZM. Firstly, it can be seen that the eight sets of extracted CZLs are close
to each other, even when the estimates of cohesive-zone sizes and positions are scattering from
1.2 mm to 1.7 mm. Actually, this observation is consistent with the numerical test results in
Figure 4.12 and Table 4.2, which shows that the relative errors of the estimated cohesive-zone
size and position are not proportional to the relative error of extracted CZL. Secondly, all the
CZLs share a similar two-stage functional shape, including an initial increasing part and a
subsequent softening part. Another important observation is that the maximum of cohesive
tractions ranges from 365 MPa to 425 MPa, smaller than the manufacturer-provided uniaxial
tensile strength, 483 MPa. It may be due to the stress multi-axiality near the crack tip [76].

127

4.5

x 10

8

91s
92s
94s
95s
96s
97s
98s
99s
Average

Cohesive traction-t (Pa)

4
3.5
3
2.5
2
1.5
1
0.5
0
0

0.2

0.4

0.6
0.8
Separation- (m)

1

1.2

1.4
x 10

-4

Location of cohesive zone and LEFM effective crack-tip (pixel)

Figure 5.7 Cohesive-zone laws extracted at different time

940
920

Extracted cohesive zone size and location
LEFM effective crack-tip location

900
880
860
1.2mm 1.3mm 1.7mm
840

1.4mm
1.7mm

820
1.6mm

800
780
760
90

1.6mm
1.7mm
91

92

93

94

95
96
Time (s)

97

98

99

100

Figure 5.8 Cohesive-zone sizes and positions estimated at different time
128

5.4 Verification of the experimental inverse solutions
In this section, the accuracy of experimental measurement of CZLs is verified. Here, the
average of the total eight inverse solutions of CZLs is considered as the representative result.
Firstly, the accuracy of extracted CZL was verified by comparing the fracture energies within
the frameworks of LEFM and cohesive-zone model (CZM). In LEFM, the energy release rate is
a critical parameter which drives the fracture process. As stated in Chapter 3, the energy release
rate can be computed with J-integral, or from the area under the load versus load-line
displacement curves [129]. In CZM, the fracture process is driven by the CZL. The cohesive
fracture energy is equivalent to the area under the CZL curve. An important connection between
the LEFM and CZM is the equivalence between the LEFM energy release rate and the cohesive
fracture energy [63]–[65] during steady-state crack propagation. In Chapter 3, the steady-state
value of energy release rate estimated by J-integral is about 25 KJ/m2. In this chapter, the
cohesive fracture energy is estimated to be 25.4 KJ/m2. As expected, the two fracture energies
are consistent.
Moreover, another verification procedure was carried out through a finite element simulation.
Using the cohesive-zone modeling capability in ABAQUS [130], the mode I translaminar
fracture experiment with an ECT specimen was replicated with the extracted CZL. Then the
finite element solution of load versus crack mouth opening displacement (CMOD) was compared
to the experimental outcomes. The simulation procedure can be considered as a "direct" problem.
Since the CZL was inversely extracted from the same experiment, the solutions of "direct"
problem which uses the extracted CZL as the input should be consistent with the experimental
results. The finite element mesh and applied boundary conditions for the finite element
129

simulation are shown in Figure 5.9, while the comparison between finite element solution and
experimental result is shown in Figure 5.10.
From Figure 5.9, it can be seen that a set of cohesive elements was inserted into the specimen,
along the extension line of crack. This setup was consistent with the self-similar crack extension
which was shown in Figure 3.7 (c). The simulation was under the plane stress condition, with the
bulk material properties provided in Table 3.1. The loading condition applied on the model also
replicated the experiment: the upper loading point was fixed, while the lower one was loaded at a
constant rate of 0.01 mm/sec. In addition, the mesh dependence of the finite element simulation
results was investigated. The result shown in Figure 5.10 is the convergent result with fine mesh.

130

Figure 5.9 Finite element model of ECT specimen with inserted cohesive elements

Figure 5.10 Comparison of the load versus crack mouth opening displacement (CMOD) curves
between the experiment and finite element simulation.
131

From Figure 5.10, a good agreement can be seen between the experimental data and finite
element simulation result. The maximum load estimated by the finite element simulation is 1258
N, which is only 6 percent higher than the experimental maximum load, 1185 N. Moreover, the
increasing and decreasing parts in load-CMOD curves of the experiment and simulation are very
close. The consistency between the two sets of results provides a strong evidence for the
accuracy of extracted CZL.
It is noteworthy that the cohesive-zone size in the finite element simulation is estimated to be
1.04 mm, which is smaller than the inverse solutions shown in Figure 5.8. The discrepancy is
mainly due to the measurement noise in the inverse solutions. From Figure 4.12 and Figure 5.8,
it has been found that the inverse solution of cohesive-zone size is sensitive to the noise in input
data. Nevertheless, it was also found that the scattered values of estimated cohesive-zone size
and position did not cause significantly scattering of the extracted CZLs. Therefore, even though
the accuracy of estimated cohesive-zone size and position is unsatisfying, this analytical inverse
method is still considered to be reliable in the extraction of CZLs, which is the most important
feature of CZM.
Finally, recall that the maximum load predicted by the R-curves is 1160 N, as presented in
Section 3.3.2.5. It can be seen that the prediction of load carrying capability using LEFM Rcurve show a great agreement with that using finite element simulation with extracted CZL.

132

5.5 Summary
In this chapter, an experiment validation was conducted for the applicability of nonlinear
separable least-squares method, an analytical inverse method which had been developed and
numerically tested in Chapter 4. Following the guidelines provided by numerical tests, this
analytical inverse method was applied to measure the CZL for a translaminar fracture test of
cross-ply composite laminates. In the experiments, the crack-tip displacement fields measured by
2-D digital image correlation (DIC) were used as the inputs for the inverse extraction of CZLs.
Eleven displacement data sets which had been measured during the steady-state crack
propagation were chosen for the inverse computation. For each data set, different initial guesses
of cohesive-zone size and position were used to check the consistency of inverse solutions.
Consequently, eight of total eleven cases have successfully obtained reasonable and consistent
inverse solutions of CZLs, while the estimates of cohesive-zone sizes and positions show some
fluctuations. Then the accuracy of inverse solutions was discussed. Firstly, it was demonstrated
that the energy release rate measured in Chapter 3 was equivalent to the cohesive fracture energy
during steady-state crack propagation. Then an advanced verification work was carried out
through the comparison between the experimental load-CMOD curves with the finite element
simulation results which used the extracted CZL as input. A good agreement was found between
the two results, which verified the accuracy of extracted CZL. Therefore, the applicability of
separable nonlinear least-squares method is validated.

133

Chapter 6
Conclusions and Future Work
6.1 Significance and contributions
The primary focus of this dissertation is the investigation of translaminar fracture mechanisms
of fiber-reinforced composite laminates. The investigations were carried out within the
frameworks of linear elastic fracture mechanics (LEFM) and cohesive-zone model (CZM),
through the development of experimental methods which can provide accurate and convenience
measurements of R-curves and cohesive-zone law (CZL), respectively. Based on the bases of
LEFM and CZM, the fiber-reinforced composite laminates were approximated as anisotropic
solids. Generally, the studies for each fracture theory can be divided into two subsequent
sections. Firstly, the experimental methods were developed. Their applicability were
experimentally validated. Secondly, the applicability of LEFM and CZM for the fracture analysis
and fracture prediction of fibre-reinforced composite laminates were investigated based on the
experimental measurements.
The researches in this dissertation were motivated by the current challenges in the studies of
translaminar fracture of fiber-reinforced composite laminates within the frameworks of LEFM
and CZM. The challenges were explicitly indicated in Chapter 2. For each fracture theory, the
thrusts of research are separately denoted as follows.
The thrusts of the research within the context of LEFM can be summarized into two aspects.
Firstly, two experimental methods based on stress intensity approach and energy approach were
developed to measure the crack-tip field parameters from the crack-tip displacement fields of

134

anisotropic solids. The crack-tip displacement fields were experimentally measured with a fullfield optical method-Digital Image Correlation (DIC). The studied crack-tip field parameters
included the stress intensity factors, the energy release rate and the effective crack length. Rcurves can be generated with these crack-tip field parameters. Secondly, the geometryindependence of R-curve was verified by comparing the crack-tip field parameters separately
estimated by the stress intensity approach and energy approach.
Moreover, the intellectual merits of research with LEFM can be clarified as follows.
(1) This research is the first attempt to verify the geometry-independence of R-curve and
hence the applicability of LEFM for the fracture study of fiber-reinforced composite
laminates, with a measurement of crack-tip full-field displacements.
(2) The developed experimental methods which provide accurate and convenient R-curve
measurements can considerably enhance the fracture predictive capability of LEFM for
fiber-reinforced composite laminates.
(3) This study is expected to provide a foundation for top-down method, thus is helpful for
the development of multi-scale fracture models for fiber-reinforced composite laminates.
For the research within the framework of CZM, two analytical inverse methods were
developed to measure the cohesive-zone law with detailed functional shape for the translaminar
fracture of fiber-reinforced composite laminates. The inverse methods were developed based on
an analytical description of the elastic fields around a crack-tip cohesive zone in anisotropic
solids. With the inverse methods, the cohesive-zone properties such as the cohesive traction,
separation and cohesive-zone law can be inversely extracted from the crack-tip displacement
fields. The accuracy and reliability of two methods were studied with a sequence of numerical
tests, and a noise-resistant inverse method was identified between the two. Based on the

135

guidelines provided by numerical tests, the applicability of the chosen method was demonstrated
through an experimental validation study with DIC.
The intellectual merits of the research about CZM can be summarized as follows.
(1) This study is the first attempt to extract CZLs with analytical inverse methods from the
cohesive-crack-tip displacement fields in macro-scale anisotropic solids.
(2) As denoted in the dissertation, the analytical inverse methods can be used not only for
mode-I fracture problem, but also for mode-II and mixed mode fracture problems.
(3) Since the analytical inverse method can provide accurate measurement of CZLs, it opens
up an avenue for the improvement of micro-mechanical fracture models with the
extracted CZLs.
Furthermore, it was found that the experimental LEFM R-curves and CZL showed a great
agreement in terms of fracture energies and maximum load predictions. It is consistent with the
understanding that quasi-brittle fracture problems such as the case in this dissertation can be
successfully analyzed with both the LEFM and CZM, which also demonstrates the validity of the
research in this dissertation.

6.2 Recommendations for future work
The breadth and depth of the research in this dissertation were limited by time and resource.
There are several more potential research topics that can be conducted in the future, based on the
fundamental theories and ideas provided in the dissertation. In this section, these possibilities are
listed, for the purpose of inspiring any other researchers who would like to make use of this
dissertation.

136

Firstly, in this dissertation, only the mode I translaminar fracture problem has been
investigated. It is mainly because of the simplicity of mode I fracture problem, which is
convenience for the illustration and verification of ideas and methods. However, as separately
stated in Chapter 3 and Chapter 4, both of the experimental methods for the measurement of
LEFM crack-tip field parameters and CZLs have the potentials to be extended to study mode II
and even mixed-mode fracture problems. Moreover, it is also noteworthy that all the theoretical
derivations for mode II and mixed-mode fracture problems have been given in this dissertation.
Therefore, any future work about the experimental validation of current methods for mode II and
mixed-mode fracture problems can be easily guided.
Secondly, the material that has been studied in this dissertation is a cross-ply glass-fiber
reinforced epoxy laminate. This material was chosen due to the fact that its quasi-brittle fracture
behavior during self-similar translaminar crack growth was well-suited for the experimental
validation of the methods for the measurement of R-curve and CZL. However, it is noteworthy
that the developed methods are applicable for the translaminar fracture analysis of any fiberreinforced composite materials which can be approximated as anisotropic solids and show selfsimilar crack growth.
Thirdly, the separable nonlinear least-squares method for the inverse extraction of CZL should
not only be applicable for anisotropic materials like fiber-reinforced composite laminates, but
may also be applicable for isotropic materials. Essentially, the prerequisite of separable nonlinear
least-squares method is a analytical characterization (i.e., mathematical expression) of the
interested object. For isotropic materials, the eigenfunction expansion for the characterization of
crack-tip fields around a cohesive zone has already been derived [75] with the Muskhelishvili
formalism of isotropic plane elasticity [109] [131]. Although the eigenfunction expansions for

137

isotropic materials and anisotropic materials are significantly different in terms of appearance,
they have been derived from the same boundary value problem. Therefore, it is reasonable to
assume that the separable nonlinear least-squares method can also be used for isotropic materials.
As represented in this dissertation, the applicability of this inverse method for isotropic materials
can be verified with a sequence of numerical tests and experimental validations.
Finally, the analytical inverse methods for the measurement of CZL can be used as an
powerful tool to study the fracture predictive capability of CZM. Generally, in CZM, the CZL is
considered as a material constant. The material dependence of CZL can be interpreted with two
different ways. On one hand, for a specific material, the CZL extracted from one specimen
configuration and loading condition should be able to predict the fracture behaviours of other
specimen geometries and loading conditions. On the other hand, the CZLs extracted from
different experiments with various specimens and loading conditions should be identical.
Although the geometry-independent fracture predictive capability of CZM has been analytically
verified, there is little experimental evidences since the experimental measurement of CZL
remains challenging. With the analytical inverse methods developed in this dissertation, the
fracture predictive capability of CZM for different structural geometries and loading conditions
can be experimentally validated with a sequence of well-designed experiments. Based on the
experimental results, topics such as the improvement of analytical inverse methods or the
enhancement of the cohesive-zone modeling technique can be investigated as a further step.

138

APPENDICES

139

Appendix A
Source Code of 2D Digital Image Correlation (DIC) Software

% This DIC code can measure large deformation by accumulating successive deformations in a
sequence of images (Number of images>=2).
% All the images are captured successively during a deformation process. The format of image is
256 grey-scale bmp file.
% The sequence of images used for DIC measurement should be kept in the same folder.
% All the deformed images which are captured successively should be named in order (e.g.
imag000, imag001, imag002,...).
% The first image must be the reference image.
% This DIC code exports a sequence of comma-separated values (csv) files.
% Each csv file contains the coordinates of data points (x, y) and displacement results (u, v),
corresponding to one deformed image.
% The csv files are automatically ordered in terms of three-digit file names (e.g. 000.csv,
001.csv, ...).
% The 000.csv is to compare reference image with itself, so displacements are all zeros.
% This function is the main function that is used to load images and export displacement result
of DIC.
clear
close all
b=dir('*.bmp');
global a
global fx
global fy
global fxy
% Users' inputs. Define the region of interest, spacing between neighbouring data points and
sizes of subsets in DIC algorithm
m0start_x=220;
% lower limit of x-coordinate to define the region of interest (in the
coordinate system of image)
m0start_y=180;
% lower limit of y-coordinate to define the region of interest (in the
coordinate system of image)
m0end_x=1070;
% upper limit of x-coordinate to define the region of interest (in the
coordinate system of image)
m0end_y=780;
% upper limit of y-coordinate to define the region of interest (in the
coordinate system of image)
stepsize=10;
% spacing between neighbouring data points
w1=100;
% even number,size of the subset used in obtaining integer displacements
w2=41;
% odd number,size of the subset used in obtaining the sub-pixel displacements

140

m0=double(imread(b(1).name));
[M N]=size(m0);
y_start=M+1-m0start_y;
y_end=M+1-m0end_y;
[x0,y0]=meshgrid(m0start_x:stepsize:m0end_x,y_start:-stepsize:y_end);
[x1,y1]=meshgrid(m0start_x:stepsize:m0end_x,y_start:-stepsize:y_end);
d=dir('*.csv');
for i=(length(d)+1):length(b),
tic;
m1=double(imread(b(i).name));
a=zeros(M-1,N-1,16)*nan;
[fx,fy,fxy] = calderivative(m1);
k=floor(i-1);

% interpolation coefficients matrix
% determine the gradient of the grey value

if k>=1
dd=dir('*.csv');
data1=csvread(dd(k).name,1,0);
ui=data1(:,3);
vi=data1(:,4);
kkk=1;
for jjj=1:((m0end_x-m0start_x)/stepsize+1)
for iii=1:((m0end_y-m0start_y)/stepsize+1)
ui1(iii,jjj)=ui(kkk);
vi1(iii,jjj)=vi(kkk);
kkk=kkk+1;
end
end
x1=x0+round(ui1);
y1=y0+round(vi1);
end
[u(:,:,i),v(:,:,i)]=dic_full_lm(m0,m1,x0,y0,x1,y1,w1,w2); % measure displacement fields u and v.
% To export the DIC result into .csv files
un=u(:,:,i);
vn=v(:,:,i);
x00=x0(:);
y00=y0(:);
unn=un(:);
vnn=vn(:);
result=[x00 y00 unn vnn];
kk=int2str(k);
str0=num2str(0);
colnames={'x','y','u','v'};
if k<10
141

cellwrite(strcat(str0,str0,kk,'.csv'),colnames);
cellwrite(strcat(str0,str0,kk,'.csv'),num2cell(result));
figure;mesh(x0,y0,un);
filename=strcat(str0,str0,kk,'_u.eps');
print('-depsc2',filename);
close;
figure;mesh(x0,y0,vn);
filename=strcat(str0,str0,kk,'_v.eps');
print('-depsc2',filename);
close;
elseif k>=10&&k<=99
cellwrite(strcat(str0,kk,'.csv'),colnames);
cellwrite(strcat(str0,kk,'.csv'),num2cell(result));
figure;mesh(x0,y0,un);
filename=strcat(str0,kk,'_u.eps');
print('-depsc2',filename);
close;
figure;mesh(x0,y0,vn);
filename=strcat(str0,kk,'_v.eps');
print('-depsc2',filename);
close;
else
cellwrite(strcat(kk,'.csv'),colnames);
cellwrite(strcat(kk,'.csv'),num2cell(result));
figure;mesh(x0,y0,un);
filename=strcat(kk,'_u.eps');
print('-depsc2',filename);
close;
figure;mesh(x0,y0,vn);
filename=strcat(kk,'_v.eps');
print('-depsc2',filename);
close;
end
i
% output the order of image that has just been processed
toc;
% output the computation time
end

function [fx,fy,fxy] = calderivative (m1)
% This function is used to determine the gradient of the grey values for the whole image
fx=(circshift(m1,[0 -1])-circshift(m1,[0 1]))/2;
fy=(circshift(m1,1)-circshift(m1,-1))/2;
fxy=(circshift(m1,[1 -1])-circshift(m1,[-1 -1])-circshift(m1,[1 1])+circshift(m1,[-1 1]))/4;
end

142

function [u,v]=dic_full_lm(m0,m1,x0,y0,x1,y1,w1,w2)
% This function is used to measure the displacements with DIC algorithm
[up,vp]=dic_full_pixel(m0,m1,x0,y0,x1,y1,w1); % measure the integer-pixel-level displacements
[m n]=size(x0);
%Incomplete compensatory program to redefine the "nan" value in the integer displacements.
if find(isnan(up))~=0
[ni,nj]=find(isnan(up));
for i=1:length(ni)
if ni(i)~=1 && ni(i)~=m
up(ni(i),nj(i))=(up(ni(i)-1,nj(i))+up(ni(i)+1,nj(i)))/2;
vp(ni(i),nj(i))=(vp(ni(i)-1,nj(i))+vp(ni(i)+1,nj(i)))/2;
else
up(ni(i),nj(i))=(up(ni(i),nj(i)-1)+up(ni(i),nj(i)+1))/2;
vp(ni(i),nj(i))=(vp(ni(i),nj(i)+1)+vp(ni(i),nj(i)+1))/2;
end
end
end
u=zeros(m,n)*nan;
v=zeros(m,n)*nan;
h=waitbar(0,'Please wait...');
j=1;
for i=1:m,
waitbar(i/m,h);
if j==1,
for j=1:n,
if i==1&&j==1
P=[up(i),vp(j),0,0,0,0]';
end
[P,C]=LM(x0(i,j),y0(i,j),P,w2,m0,m1); % calculate the deformation vector P and
correlation coefficient C
u(i,j)=P(1);
v(i,j)=P(2);
end
else
for j=n:-1:1,
[P,C]=LM(x0(i,j),y0(i,j),P,w2,m0,m1); % calculate the deformation vector P and
correlation coefficient C
u(i,j)=P(1);
v(i,j)=P(2);
end
end
end
close(h);
end
143

function [u,v]=dic_full_pixel(m0,m1,x0,y0,x1,y1,w)
% This function is used to measure the integer-pixel-level displacements
[M N]=size(m0);[m n]=size(x0);
u=zeros(m,n)*nan;
v=zeros(m,n)*nan;
h=waitbar(0,'Please wait...');
for i=1:m,
waitbar(i/m,h);
for j=1:n,
i0=(M+1-y0(i,j)-w/2); i1=M+1-y0(i,j)+w/2-1;
j0=(x0(i,j)-w/2); j1=x0(i,j)+w/2-1;
p0=(M+1-y1(i,j)-w/2); p1=M+1-y1(i,j)+w/2-1;
q0=(x1(i,j)-w/2); q1=x1(i,j)+w/2-1;
if i0>=1 & i0<=M & j0>=1 & j0<=N & i1>=1 & i1<=M & j1>=1 & j1<=N &...
p0>=1 & p0<=M & q0>=1 & q0<=N & p1>=1 & p1<=M & q1>=1 & q1<=N,
dm0=m0(i0:i1,j0:j1);
dm1=m1(p0:p1,q0:q1);
[du,dv]=decorr(dm0,dm1);
% Measure the displacements in Fourier's domain
if (length(du)*length(dv)==1)
u(i,j)=du+x1(i,j)-x0(i,j);
v(i,j)=dv+y1(i,j)-y0(i,j);
else
disp('temp=0');
end
end
end
end
close(h);
end

function [du,dv]=decorr(dm0,dm1)
% This function is used to measure the integer-pixel-level displacements in Fourier's domain
[M,N]=size(dm0);
dm0=double(dm0); dm1=double(dm1);
c=abs(fftshift(ifft2(conj(fft2(dm0)).*fft2(dm1))));
cc=sum(sum((dm0-dm1).*(dm0-dm1)));
if cc==0
du=0;
dv=0;
else
[di,dj]=find(c==max(c(:)));
du=dj-(N/2+1);
dv=-(di-(M/2+1));
144

end
end
function [P,C]=LM(X0,Y0,P0,w,m0,m1)
% This function is used to calculate the deformation vector P and correlation coefficient C for
one point using Levenberg-Marquardt algorithm.
eps1=double(1e-6); % convergent limit 1
eps2=double(1e-6); % convergent limit 2
lamda=0.01;
% set an initial value of the damping factor for the Levenberg-Marquardt
algorithm
update=1; % to judge if we should wait until lamda changes to a suitable value
cnt=0;
MAXCNT=500; % maximal cycle index
while cnt<MAXCNT
% prevent infinite cycle
if update==1
% lamda is suitable
if cnt==0
[GC,H,C0]=grhessianZNSSD(P0,X0,Y0,m0,m1,w); % calculate correlation coefficient C0,
Hessian matrix H and coefficient gradient GC for next iterative step
else
[GC,H]=grhessianZNSSD1(P0,X0,Y0,m0,m1,w); % calculate Hessian matrix H and
coefficient gradient GC
end
end
Him=H+(lamda*diag(diag(H)));
P=P0-Him\GC;
% the iterative formula to calculate deformation vector
[C]=grhessianZNSSD2(P,X0,Y0,m0,m1,w); % calculate correlation coefficient C
if norm(GC)<=eps1 || norm(P-P0)<=eps2*(norm(P0)+eps2)
break;
end
if C<C0
P0=P;
C0=C;
lamda=lamda/10;
update=1;
else
lamda=lamda*10;
update=0;
end
cnt=cnt+1;
end
if cnt==MAXCNT
disp('divergence')
end
end

145

function [GC,H,C]=grhessianZNSSD(P,X0,Y0,m0,m1,w)
% This function is used to calculate Hessian matrix H, correlation coefficient C and coefficient
gradients GC of the subset with
% the center of (X0,Y0)---under x-y coordinate system
% The zero-normalized sum of squared differences (ZNSSD) is used
global a
H=zeros(6,6);
GC=[0 0 0 0 0 0]';
[M,N]=size(m0);
qq=1;
for ii=-(w-1)/2:(w-1)/2,
pp=1;
for jj=-(w-1)/2:(w-1)/2,
X=X0+ii+P(1)+P(3)*ii+P(5)*jj;
Y=Y0+jj+P(2)+P(4)*ii+P(6)*jj; %X0 and Y0 are the x-y coordinates of the center point of
this subset
i=floor(M+1-Y);
j=floor(X);
% to judge in which unit square the point is----unit(i,j)
if isnan(a(i,j,1))==1,
a(i,j,:)=bispcoeff(i,j,m1);
end
aa=[a(i,j,1) a(i,j,2) a(i,j,3) a(i,j,4)...
a(i,j,5) a(i,j,6) a(i,j,7) a(i,j,8)...
a(i,j,9) a(i,j,10) a(i,j,11) a(i,j,12)...
a(i,j,13) a(i,j,14) a(i,j,15) a(i,j,16)];
xx=X-j; yy=Y+i-M;
xx2=xx^2;xx3=xx^3;yy2=yy^2;yy3=yy^3;
g(w+1-pp,qq)=aa*[1 xx xx2 xx3 yy xx*yy yy*xx2 yy*xx3 yy2 xx*yy2 xx2*yy2 xx3*yy2 yy3
xx*yy3 xx2*yy3 xx3*yy3]';
gx(w+1-pp,qq)=aa*[0 1 2*xx 3*xx2 0 yy 2*xx*yy 3*xx2*yy 0 yy2 2*xx*yy2 3*xx2*yy2 0 yy3
2*xx*yy3 3*xx2*yy3]';
gy(w+1-pp,qq)=aa*[0 0 0 0 1 xx xx2 xx3 2*yy 2*xx*yy 2*xx2*yy 2*xx3*yy 3*yy2 3*xx*yy2
3*xx2*yy2 3*xx3*yy2]';
gp(w+1-pp,qq,1)=gx(w+1-pp,qq);
gp(w+1-pp,qq,2)=gy(w+1-pp,qq);
gp(w+1-pp,qq,3)=gx(w+1-pp,qq)*ii;
gp(w+1-pp,qq,4)=gy(w+1-pp,qq)*ii;
gp(w+1-pp,qq,5)=gx(w+1-pp,qq)*jj;
gp(w+1-pp,qq,6)=gy(w+1-pp,qq)*jj;
pp=pp+1;
end
qq=qq+1;
end
F=sum(sum(m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)));%scalar
Fm=F/w^2;
%scalar
F1=m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)-Fm; %%matrix
146

FF1=sum(sum(F1.^2)); %scalar
FF=sqrt(FF1); %scalar
gm=mean(g(:));
%scalar
g1=g-gm; %%matrix
gg1=sum(sum(g1.^2)); %scalar
gg=sqrt(gg1);
%scalar
Fs=F1/FF-g1/gg; %matrix
Fss=Fs/gg;
%matrix
for mm=1:6,
GC(mm)=-2*sum(sum(Fss.*gp(:,:,mm)));
end
gg2=2/gg1; %scalar
for mm=1:6,
for nn=1:6,
H(mm,nn)=gg2*sum(sum(gp(:,:,mm).*gp(:,:,nn)));
end
end
C=sum(sum(Fs.^2));
end

function [GC,H]=grhessianZNSSD1(P,X0,Y0,m0,m1,w)
% This funciton is used to calculate Hessian matrix H and coefficient gradients GC of the subset
with
% the center of (X0,Y0)---under x-y coordinate system
% The zero-normalized sum of squared differences (ZNSSD) is used
global a
H=zeros(6,6);
GC=[0 0 0 0 0 0]';
[M,N]=size(m0);
qq=1;
for ii=-(w-1)/2:(w-1)/2,
pp=1;
for jj=-(w-1)/2:(w-1)/2,
X=X0+ii+P(1)+P(3)*ii+P(5)*jj;
Y=Y0+jj+P(2)+P(4)*ii+P(6)*jj;%X0 and Y0 are the x-y coordinates of the center point of
this subset
i=floor(M+1-Y);
j=floor(X);
% to judge in which unit square the point is----unit(i,j)
if isnan(a(i,j,1))==1,
a(i,j,:)=bispcoeff(i,j,m1);
end
aa=[a(i,j,1) a(i,j,2) a(i,j,3) a(i,j,4)...
a(i,j,5) a(i,j,6) a(i,j,7) a(i,j,8)...
a(i,j,9) a(i,j,10) a(i,j,11) a(i,j,12)...
a(i,j,13) a(i,j,14) a(i,j,15) a(i,j,16)];
147

xx=X-j; yy=Y+i-M;
xx2=xx^2;xx3=xx^3;yy2=yy^2;yy3=yy^3;
g(w+1-pp,qq)=aa*[1 xx xx2 xx3 yy xx*yy yy*xx2 yy*xx3 yy2 xx*yy2 xx2*yy2 xx3*yy2 yy3
xx*yy3 xx2*yy3 xx3*yy3]';
gx(w+1-pp,qq)=aa*[0 1 2*xx 3*xx2 0 yy 2*xx*yy 3*xx2*yy 0 yy2 2*xx*yy2 3*xx2*yy2 0 yy3
2*xx*yy3 3*xx2*yy3]';
gy(w+1-pp,qq)=aa*[0 0 0 0 1 xx xx2 xx3 2*yy 2*xx*yy 2*xx2*yy 2*xx3*yy 3*yy2 3*xx*yy2
3*xx2*yy2 3*xx3*yy2]';
gp(w+1-pp,qq,1)=gx(w+1-pp,qq);
gp(w+1-pp,qq,2)=gy(w+1-pp,qq);
gp(w+1-pp,qq,3)=gx(w+1-pp,qq)*ii;
gp(w+1-pp,qq,4)=gy(w+1-pp,qq)*ii;
gp(w+1-pp,qq,5)=gx(w+1-pp,qq)*jj;
gp(w+1-pp,qq,6)=gy(w+1-pp,qq)*jj;
pp=pp+1;
end
qq=qq+1;
end
F=sum(sum(m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)));%scalar
Fm=F/w^2;
%scalar
F1=m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)-Fm; %matrix
FF1=sum(sum(F1.^2)); %scalar
FF=sqrt(FF1); %scalar
gm=mean(g(:));
%scalar
g1=g-gm; %%matrix
gg1=sum(sum(g1.^2)); %scalar
gg=sqrt(gg1);
%scalar
Fs=F1/FF-g1/gg; %matrix
Fss=Fs/gg;
%matrix
for mm=1:6,
GC(mm)=-2*sum(sum(Fss.*gp(:,:,mm)));
end
gg2=2/gg1; %scalar
for mm=1:6,
for nn=1:6,
H(mm,nn)=gg2*sum(sum(gp(:,:,mm).*gp(:,:,nn)));
end
end
end

function [C]=grhessianZNSSD2(P,X0,Y0,m0,m1,w)
% To calculate the correlation coefficient C of the subset with
% the center of (X0,Y0)---under x-y coordinate system
% The zero-normalized sum of squared differences (ZNSSD) is used
global a
148

[M,N]=size(m0);
qq=1;
for ii=-(w-1)/2:(w-1)/2,
pp=1;
for jj=-(w-1)/2:(w-1)/2,
X=X0+ii+P(1)+P(3)*ii+P(5)*jj;
Y=Y0+jj+P(2)+P(4)*ii+P(6)*jj;%X0 and Y0 are the x-y coordinates of the center point of
this subset
i=floor(M+1-Y);
j=floor(X);% to judge in which unit square the point is----unit(i,j)
if isnan(a(i,j,1))==1,
a(i,j,:)=bispcoeff(i,j,m1);
end
aa=[a(i,j,1) a(i,j,2) a(i,j,3) a(i,j,4)...
a(i,j,5) a(i,j,6) a(i,j,7) a(i,j,8)...
a(i,j,9) a(i,j,10) a(i,j,11) a(i,j,12)...
a(i,j,13) a(i,j,14) a(i,j,15) a(i,j,16)];
xx=X-j; yy=Y+i-M;
xx2=xx^2;xx3=xx^3;yy2=yy^2;yy3=yy^3;
g(w+1-pp,qq)=aa*[1 xx xx2 xx3 yy xx*yy yy*xx2 yy*xx3 yy2 xx*yy2 xx2*yy2 xx3*yy2 yy3
xx*yy3 xx2*yy3 xx3*yy3]';
pp=pp+1;
end
qq=qq+1;
end
F=sum(sum(m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)));%scalar
Fm=F/w^2;
%scalar
F1=m0(M+1-Y0-(w-1)/2:M+1-Y0+(w-1)/2,X0-(w-1)/2:X0+(w-1)/2)-Fm; %%matrix
FF1=sum(sum(F1.^2)); %scalar
FF=sqrt(FF1); %scalar
gm=mean(g(:));
%scalar
g1=g-gm; %matrix
gg1=sum(sum(g1.^2)); %scalar
gg=sqrt(gg1);
%scalar
Fs=F1/FF-g1/gg; %matrix
C=sum(sum(Fs.^2));
end

149

Appendix B
Derivation of the Asymptotic Expansion of Semi-infinite Crack Problem in
Anisotropic Solids

Consider a semi-infinite traction-free crack problem in an anisotropic solid, as shown in
Figure 3.1. An asymptotic expansion of the crack-tip elastic fields can be solved from the
boundary conditions as follows.
For the derivation, let
ti   i 2   Li f' ( z ) ,

(B.1)



which can be rewritten in the matrix notation as

t( z )  Lf ' ( z )  Lf ' ( z ) .

(B.2)

Firstly, the continuity condition of traction along the entire x -axis can be expressed as
t( x)  t( x)  0

(  x  )

(B.3)

where the superscript plus (  ) and minus (  ) signs denote the evaluation of the function value
at x  i as  is positive and approaches to zero, t(x )  and t(x ) can be represented as follows,

t( x)  Lf ' ( x)  Lf ' ( x)

(B.4)

t( x)  Lf ' ( x)  Lf ' ( x)

(B.5)

Substitute the equation (B.4) and (B.5) into (B.2), we can get

t( x)  t( x)  Lf ' ( x)  Lf ' ( x)  (Lf ' ( x)  Lf ' ( x) )  0 .

(B.6)

The above equation can be rearranged as

Lf ( x)  Lf ( x)  Lf ( x)  Lf ( x)  0 .
'

'



'

150

'



(B.7)

From the equation (B.7), a general solution can be derived as [109]

Lf ' ( z )  Lf ' ( z )  ig(z) ,

(B.8)

where g(z) is an entire function.
Secondly, the traction-free condition on the semi-infinite crack leads to a homogeneous
Hilbert arc problem [109] with a straight cut along the negative x -axis, which can be represented
as follows,
t ( x)   t ( x)   0

(  x  0)

(B.9)

Substitute the equation (B.4) and (B.5) into (B.9), then we can get

Lf (x)  Lf (x)  Lf (x)  Lf (x)  0 .
'



'

'



'

(B.10)

A general solution in terms of eigenfunction expansion can be represented as [109]
Lf ' ( z )  Lf ' ( z ) 

h( z )
X( z )

(B.11)

where X( z )  z for semi-infinite sharp crack with the crack-tip at z=0 (  x  0) .
Therefore, an asymptotic expansion of the crack-tip elastic fields can be obtained as
Lf ( z ) 
N

N

n 0

n0

1  g( z )

 i h( z ) 

2 z


(B.12)

where g( z )   a n z n and h( z )   b n z n as N   , and both a n and b n are column vectors with
three real-number elements.

151

Appendix C
Elastic Matrices for Orthotropic Solids under Plane Stress Condition

The characteristic equation for an orthotropic material under plane stress condition, of which
the principal axes align with the coordinate axes, is
1

λ p 4  2ρ λ 2 p 2  1  0

(C.1)

where the two non-dimensional parameters λ and ρ are expressed in terms of the Engineering
constants as

λ

1

E2
 21 ) .
, ρ  E1E2 (
E1
212 E2

(C.2)

The characteristic roots with positive imaginary parts p1 and p2 are explicitly given [106] as


p1 

i

2


p1 

1
4



1
4

2



( 1  ρ  ρ 1) , p2 

i

2


(i 1  ρ  1  ρ ) , p2 

1
4



1
4

2

( 1  ρ  ρ 1) ,

(i 1  ρ  1  ρ ) ,

for

for

1 ρ  

1  ρ  1

(C.3)

(C.4)

The elements of the matrices A , L , B and H are then given by
 p12  12
 E
1
A
2
1



21 p1
 E p
2 1


p22  12 
E1 
.
1  21 p22 
E2 p2 

  p  p2 
.
L 1
1 
 1

152

(C.5)

(C.6)

  i ( p1  p2 )

E1
B
i
(

p
 21 1 p2  1)
 E2 p1 p2

  2i ( p1  p2 )

E1
H

0


153

 i ( 12  p1 p2 ) 

E1

i ( p1  p2 )  .

E2 p1 p2




2i ( p1  p2 )  .
E2 p1 p2 

(C.7)

0

(C.8)

Appendix D
Additional Results of Translaminar Fracture Tests with Fiber-reinforced
Composite Laminates

The composite material studied in this dissertation was Cyply 1002 (formerly known as
Scotchply 1002), a glass-fiber reinforced epoxy laminate, manufactured by Cytec Engineered
Materials. In the translaminar fracture tests following ASTM E1922 [119], a total of six
specimens with same geometries but different lay-up sequences were tested. The material
properties of the six different specimens are shown in the following table.

Table D.1 Elastic properties of composite laminates used in the translaminar fracture tests
Specimen

Lay-up

E1 [GPa]

E2 [GPa]

12 [GPa]

12

UD00

[90/0/909/0/90]T

13.2

31.2

3

0.096

UD45

[45/-45/459/45/45]T

9.38

9.38

8.18

0.563

UD90

[0/90/09/90/0]T

31.2

13.2

3

0.226

CP00

[(90/0)6/90]T

22.14

22.39

3.17

0.126

CP45

[(45/-45)6/45]T

10.23

10.16

9.88

0.601

CP90

[ (90 / 0)3 / 90 ]T

22.39

22.14

3.17

0.127

In the above table, “UD” stands for “unidirectional”, and “CP” stands for “cross-ply”. The
loading direction is along the axis 2. The experimental setups and procedures for all the
specimens were the same. Figure A.1 shows the images of specimens which were captured at the
end of fracture tests. It can be seen that only UD90, CP00 and CP90 show self-similar crack
154

growth. Nevertheless, a significant delamination (i.e. interlaminar fracture) was observed at the
back surface of specimen CP00. Therefore, only UD90 and CP90 can be used for the study of
this dissertation, which requires translaminar fracture and self-similar crack propagation. The
results of CP90 were presented in Chapter 3.
Figure D.2 shows the R-curves of UD90 which was measured with the experimental methods
developed in this dissertation. It can be seen that the R-curves measured by the nonlinear leastsquares method and conservation integrals are consistent, in spite of the fluctuations at the
beginning stage of R-curve measured by the conservation integrals. These fluctuations were
caused by the low signal-to-noise ratio in the displacements. The conservation integrals are more
sensitive to noise because line integrals with limited data points were used in the computation.

(a)

(b)

(c)

(d)

(e)

(f)

Figure D.1 Images of specimens captured at the end of translaminar fracture tests: (a) UD00 (b)
UD45 (c) UD90 (d) CP00 (e) CP45 (f) CP90.

155

Figure D.2 R-curves in terms of stress intensity factor, for specimen UD90.

156

Appendix E
Derivation of the Eigenfunction Expansion of Cohesive Crack-tip fields in
Anisotropic Solids

Consider a traction-free semi-infinite crack with a cohesive zone in an anisotropic solid, as
shown in Figure 5.1. An eigenfunction expansion of the elastic fields around crack-tip cohesive
zone can be solved from the boundary conditions as follows.

For the derivation, let
ti   i 2   Li f' ( z ) ,

(E.1)



which can be rewritten in the matrix notation as

t( z )  Lf ' ( z )  Lf ' ( z ) .

(E.2)

Firstly, the continuity condition of traction along the entire x -axis can be expressed as
t( x)  t( x)  0

(  x  )

(E.3)

where the superscript plus (  ) and minus (  ) signs denote the evaluation of the function value
at x  i as  is positive and approaches to zero, t(x )  and t(x ) can be represented as follows,

t( x)  Lf ' ( x)  Lf ' ( x)

(E.4)

t( x)  Lf ' ( x)  Lf ' ( x)

(E.5)

Substitute the equation (E.4) and (E.5) into (E.2), we can get

t( x)  t( x)  Lf ' ( x)  Lf ' ( x)  (Lf ' ( x)  Lf ' ( x) )  0 .
The above equation can be rearranged as

157

(E.6)

Lf ' ( x)  Lf '( x)  Lf ' ( x)  Lf '( x)




 0.

(E.7)

From the equation (B.7), a general solution can be derived as [109]
Lf ' ( z )  Lf '( z )  is( z ) ,

(E.8)

where s(z ) is an entire function.
Secondly, the elastic fields around the cohesive zone can be described by the superposition of
two stress-bounded sharp-crack-tip elastic fields with respect to x  c and x  c , as stated by
Hong and Kim [75]. Here the stress-bounded means that there is no singular term in the
eigenfunctions. The eigenfunction expansion is, then, expressed as,
Lf ' ( z )  Lf '( z )  z  cq( z )  z  cr( z )

(E.9)

where q(z ) and r (z ) are entire functions.
Therefore, the eigenfunction expansion of the crack-tip elastic fields are derived from (B.8)
and (B.9) as
Lf ( z ) 

1
2



z  cq( z )  z  cr( z )  i s( z )

N

N

N

n0

n0

n0



(E.10)

where q( z )   q nQn ( z ) , r ( z )   rn Rn ( z ) and s( z )   s n S n ( z ) as N   , q n , rn and s n are
column vectors with real-number elements, and Qn (z ) , Rn (z ) and Sn (z ) are entire functions.
From (E.2) to (E.5), the closing-traction distribution t (x) within the cohesive zone (
 c  x  c ), can be expressed in terms of eigenfunction expansions as

t( x)  t( x)   t( x)   Lf ' ( x)   Lf '( x)   x  cq( x) .

(E.11)

It is noteworthy that the above derivation makes use of the fact that the line segment  c  x  c is
the branch cut for the function z  cr(z ) , but not the branch cut for the function z  cq(z ) .
In addition, the jump of separation-gradient in the cohesive zone (  c  x  c ) is defined as
158

b( x)  u' ( x)  u' ( x)  .

(E.12)

From the equation (3.3), the displacement gradients is expressed as
u' ( z )  Αf ' ( z )  Af ' ( z ) .

(E.13)

Therefore, the jump of separation-gradient in the cohesive zone b(x) can be derived as
b( x)  [(Af ' ( x)  Af '( x) )  (Af ' ( x)  Af '( x) )]  H c  xr( x)

(E.14)

where H is a Stroh’s matrix defined in the equation (3.12).
Finally, the separation in the cohesive zone δ(x) can be obtained by integrating (E.14) as
follows,
c

δ( x)  x b( ) d   H x c   r( ) d 
C

159

(E.15)

BIBLIOGRAPHY

160

BIBLIOGRAPHY

[1] H. L. Cox, “The elasticity and strength of paper and other fibrous materials,” British
Journal of Applied Physics, vol. 3, no. 3, pp. 72–79, Mar. 1952.
[2] H. E. Daniels, “The Statistical Theory of the Strength of Bundles of Threads. I,” Proc. R.
Soc. Lond. A, vol. 183, no. 995, pp. 405–435, Jun. 1945.
[3] F. Hild, A. Burr, and F. A. Leckie, “Fiber breakage and fiber pull-out of fiber-reinforced
ceramic-matrix composites,” European journal of mechanics. A. Solids, vol. 13, no. 6, pp.
731–749.
[4] F. Hild, A. Burr, and F. A. Leckie, “Matrix cracking and debonding of ceramic-matrix
composites,” International Journal of Solids and Structures, vol. 33, no. 8, pp. 1209–1220,
Mar. 1996.
[5] A. Burr, F. Hild, and F. A. Leckie, “Continuum description of damage in ceramic-matrix
composites,” European journal of mechanics. A. Solids, vol. 16, no. 1, pp. 53–78.
[6] J. Aveston, G. Cooper, and A. Kelly, “Single and Multiple Fracture in: The properties of
fiber composites,” presented at the Conference proceedings, Guildford, U.K., 1971, pp. 15–
24.
[7] L. Mishnaevsky, Computational Mesomechanics of Composites: Numerical Analysis of the
Effect of Microstructures of Composites on Their Strength and Damage Resistance. John
Wiley & Sons, 2007.
[8] G. Barenblatt, “The formation of equilibrium cracks during brittle fracture: general ideas
and hypothesis, axially symmetric cracks,” Appl. Math. Mech., vol. 23, pp. 622–636, 1959.
[9] D. S. Dugdale, “Yielding of steel sheets containing slits,” Journal of the Mechanics and
Physics of Solids, vol. 8, no. 2, pp. 100–104, May 1960.
[10] A. Hillerborg, M. Modeer, and P. E. Petersson, “Analysis of crack formation and crack
growth in concrete by means of fracture mechanics and finite elements,” Cement and
concrete research, vol. 6, no. 6, pp. 773–781, 1976.
[11] G. C. Sih, P. C. Paris, and G. R. Irwin, “On cracks in rectilinearly anisotropic bodies,”
International Journal of Fracture, vol. 1, no. 3, pp. 189–203, 1965.
[12] G. Bao, S. Ho, Z. Suo, and B. Fan, “The role of material orthotropy in fracture specimens
for composites,” International Journal of Solids and Structures, vol. 29, no. 9, pp. 1105–
1116, 1992.
[13] G. Bao, “Corrigenda,” Int. J. Solids Struct., vol. 29, no. 16, pp. 2115–2115, 1992.

161

[14] L. Mishnaevsky Jr. and P. Brøndsted, “Micromechanical modeling of damage and fracture
of unidirectional fiber reinforced composites: A review,” Computational Materials Science,
vol. 44, no. 4, pp. 1351–1359, Feb. 2009.
[15] G. Bao and Z. Suo, “Remarks on Crack-Bridging Concepts,” Applied Mechanics Reviews,
vol. 45, no. 8, pp. 355–366, 1992.
[16] M. E. Waddoups, J. R. Eisenmann, and B. E. Kaminski, “Macroscopic Fracture Mechanics
of Advanced Composite Materials,” Journal of Composite Materials, vol. 5, no. 4, pp. 446–
454, Oct. 1971.
[17] J. M. Whitney and R. J. Nuismer, “Stress Fracture Criteria for Laminated Composites
Containing Stress Concentrations,” Journal of Composite Materials, vol. 8, no. 3, pp. 253–
265, Jan. 1974.
[18] R. Nuismer and J. Whitney, “Uniaxial Failure of Composite Laminates Containing Stress
Concentrations,” in Fracture Mechanics of Composites, 1975.
[19] R. F. Karlak, “Hole Effects in a Related Series of Symmetrical Laminates,” Proceedings of
Failure Modes in Composites, vol. IV, The Metallurgical Society of AIME, no. Chicago, IL,
pp. 105–117, 1977.
[20] R. B. Pipes, R. C. Wetherhold, and J. W. Gillespie, “Notched Strength of Composite
Materials,” Journal of Composite Materials, vol. 13, no. 2, pp. 148–160, Apr. 1979.
[21] R. B. Pipes, J. W. Gillespie, and R. C. Wetherhold, “Superposition of the notched strength
of composite laminates,” Polymer Engineering & Science, vol. 19, no. 16, pp. 1151–1155,
1979.
[22] R. B. Pipes, R. C. Wetherhold, and J. W. J. Gillespie, “Macroscopic fracture of fibrous
composites,” Materials Science and Engineering, vol. 45, no. 3, pp. 247–253, Oct. 1980.
[23] J. Awerbuch and M. S. Madhukar, “Notched Strength of Composite Laminates: Predictions
and Experiments--A Review,” Journal of Reinforced Plastics and Composites, vol. 4, no. 1,
pp. 3–159, Jan. 1985.
[24] M. J. Laffan, S. T. Pinho, P. Robinson, and A. J. McMillan, “Translaminar fracture
toughness testing of composites: A review,” Polymer Testing, vol. 31, no. 3, pp. 481–489,
May 2012.
[25] J. F. Mandell, F. J. Mcgarry, S. S. Wang, and J. Im, “Stress Intensity Factors for
Anisotropic Fracture Test Specimens of Several Geometries,” Journal of Composite
Materials, vol. 8, no. 2, pp. 106–116, Apr. 1974.
[26] J. Sweeney, “Finite-Width correction factors for sen testing of orthotropic materials in
opening mode,” The Journal of Strain Analysis for Engineering Design, vol. 21, no. 2, pp.
99–107, Apr. 1986.

162

[27] J. Sweeney, “Revised values of finite width correction factors for sen testing of orthotropic
materials,” The Journal of Strain Analysis for Engineering Design, vol. 23, no. 4, pp. 227–
228, Oct. 1988.
[28] N. P. O’Dowd and C. F. Shih, “Family of crack-tip fields characterized by a triaxiality
parameter—I. Structure of fields,” Journal of the Mechanics and Physics of Solids, vol. 39,
no. 8, pp. 989–1015, 1991.
[29] N. P. O’Dowd and C. F. Shih, “Family of crack-tip fields characterized by a triaxiality
parameter—II. Fracture applications,” Journal of the Mechanics and Physics of Solids, vol.
40, no. 5, pp. 939–963, 1992.
[30] A. Wolfenden, T. Anderson, and R. Dodds, “Specimen Size Requirements for Fracture
Toughness Testing in the Transition Region,” Journal of Testing and Evaluation, vol. 19,
no. 2, p. 123, 1991.
[31] R. H. D. Jr, T. L. Anderson, and M. T. Kirk, “A framework to correlate a/W ratio effects on
elastic-plastic fracture toughness (J c ),” Int J Fract, vol. 48, no. 1, pp. 1–22, Mar. 1991.
[32] M. J. Laffan, S. T. Pinho, P. Robinson, and L. Iannucci, “Measurement of the in situ ply
fracture toughness associated with mode I fibre tensile failure in FRP. Part I: Data
reduction,” Composites Science and Technology, vol. 70, no. 4, pp. 606–613, Apr. 2010.
[33] M. J. Laffan, S. T. Pinho, P. Robinson, and A. J. McMillan, “Translaminar fracture
toughness testing of composites: A review,” Polymer Testing, vol. 31, no. 3, pp. 481–489,
May 2012.
[34] S. R. McNeill, W. H. Peters, and M. A. Sutton, “Estimation of stress intensity factor by
digital image correlation,” Engineering Fracture Mechanics, vol. 28, no. 1, pp. 101–112,
1987.
[35] R. Sanford, “Determining fracture parameters with full-field optical methods,”
Experimental Mechanics, vol. 29, no. 3, pp. 241–247, 1989.
[36] J. Réthoré, A. Gravouil, F. Morestin, and A. Combescure, “Estimation of mixed-mode
stress intensity factors using digital image correlation and an interaction integral,”
International Journal of Fracture, vol. 132, no. 1, pp. 65–79, Mar. 2005.
[37] S. Yoneyama, Y. Morimoto, and M. Takashi, “Automatic Evaluation of Mixed‐mode
Stress Intensity Factors Utilizing Digital Image Correlation,” Strain, vol. 42, no. 1, pp. 21–
29, Feb. 2006.
[38] P. Lopez-Crespo, A. Shterenlikht, E. A. Patterson, J. R. Yates, and P. J. Withers, “The
Stress Intensity of Mixed Mode Cracks Determined by Digital Image Correlation,” The
Journal of Strain Analysis for Engineering Design, vol. 43, no. 8, pp. 769–780, Aug. 2008.
[39] M. Williams, “On the stress distribution at the base of a stationary crack,” Journal of
Applied Mechanics, vol. 24, pp. 109–114, 1957.
163

[40] J. R. Rice, “A Path Independent Integral and the Approximate Analysis of Strain
Concentration by Notches and Cracks,” Journal of Applied Mechanics, vol. 35, no. 2, p.
379, 1968.
[41] F. H. K. Chen and R. T. Shield, “Conservation laws in elasticity of the J-integral type,”
Zeitschrift für Angewandte Mathematik und Physik (ZAMP), vol. 28, no. 1, pp. 1–22, 1977.
[42] C. Liu, A. J. Rosakis, R. W. Ellis, and M. G. Stout, “A Study of the Fracture Behavior of
Unindirectional Fiber- Reinforced Composite Using Coherent Gradient Sensing (CGS)
Interferometry,” International Journal of Fracture, vol. 90, no. 4, pp. 355–382, 1998.
[43] G. P. Mogadpalli and V. Parameswaran, “Determination of Stress Intensity Factor for
Cracks in Orthotropic Composite Materials using Digital Image Correlation,” Strain, vol.
44, no. 6, pp. 446–452, Dec. 2008.
[44] G. Catalanotti, P. P. Camanho, J. Xavier, C. G. Dávila, and A. T. Marques, “Measurement
of resistance curves in the longitudinal failure of composites using digital image
correlation,” Composites Science and Technology, vol. 70, no. 13, pp. 1986–1993, Nov.
2010.
[45] S. H. Ju, C. Y. Chiu, and B. J. Jhao, “Determination of SIFs, crack‐tip coordinates and
crack angle of anisotropic materials,” Fatigue & Fracture of Engineering Materials &
Structures, vol. 33, no. 1, pp. 43–53, Jan. 2010.
[46] A. Needleman, “An analysis of tensile decohesion along an interface,” Journal of the
Mechanics and Physics of Solids, vol. 38, no. 3, pp. 289–324, 1990.
[47] V. Tvergaard and J. Hutchinson, “The relation between crack growth resistance and
fracture process parameters in elastic-plastic solids,” Journal of the Mechanics and Physics
of Solids, vol. 40, no. 6, pp. 1377–1397, 1992.
[48] M. G. A. Tijssens, E. Giessen, and L. J. Sluys, “Simulation of mode I crack growth in
polymers by crazing,” International Journal of Solids and Structures, vol. 37, no. 48, pp.
7307–7327, 2000.
[49] X.-P. Xu and A. Needleman, “Void nucleation by inclusion debonding in a crystal matrix,”
Modelling and Simulation in Materials Science and Engineering, vol. 1, no. 2, pp. 111–132,
Jan. 1993.
[50] Q. Yang and B. Cox, “Cohesive models for damage evolution in laminated composites,”
International Journal of Fracture, vol. 133, no. 2, pp. 107–137, May 2005.
[51] P. Gu, “Notch sensitivity of fiber-reinforced ceramics,” International Journal of Fracture,
vol. 70, no. 3, pp. 253–266, 1993.
[52] Z. Suo, M. Ortiz, and A. Needleman, “Stability of solids with interfaces,” Journal of the
Mechanics and Physics of Solids, vol. 40, no. 3, pp. 613–640, 1992.

164

[53] J. R. Rice, N. Lapusta, and K. Ranjith, “Rate and state dependent friction and the stability
of sliding between elastically deformable solids,” Journal of the Mechanics and Physics of
Solids, vol. 49, no. 9, pp. 1865–1898, 2001.
[54] N. Chandra, H. Li, C. Shet, and H. Ghonem, “Some issues in the application of cohesive
zone models for metal–ceramic interfaces,” International Journal of Solids and Structures,
vol. 39, no. 10, pp. 2827–2855, 2002.
[55] K. Y. Volokh, “Comparison between cohesive zone models,” Communications in
Numerical Methods in Engineering, vol. 20, no. 11, pp. 845–856, 2004.
[56] G. Alfano, “On the influence of the shape of the interface law on the application of
cohesive-zone models,” Composites Science and Technology, vol. 66, no. 6, pp. 723–730,
May 2006.
[57] S. H. Song, G. H. Paulino, and W. G. Buttlar, “Influence of the Cohesive Zone Model
Shape Parameter on Asphalt Concrete Fracture Behavior,” AIP Conference Proceedings,
vol. 973, no. 1, pp. 730–735, Feb. 2008.
[58] P.-E. Petersson, Crack Growth and Development of Fracture Zones in Plain Concrete and
Similar Materials. Division of Building Materials, Lund Institute of Technology, 1981.
[59] J. G. M. Van Mier and M. R. A. Van Vliet, “Uniaxial tension test for the determination of
fracture parameters of concrete: state of the art,” Engineering fracture mechanics, vol. 69,
no. 2, pp. 235–247, 2002.
[60] J. Yang and G. Fischer, “Simulation of the tensile stress-strain behavior of strain hardening
cementitious composites,” Measuring, Monitoring and Modeling Concrete Properties, pp.
25–31, 2006.
[61] K. Park, G. H. Paulino, and J. R. Roesler, “Determination of the kink point in the bilinear
softening model for concrete,” Engineering Fracture Mechanics, vol. 75, no. 13, pp. 3806–
3818, Sep. 2008.
[62] K. Park, G. H. Paulino, and J. Roesler, “Cohesive fracture model for functionally graded
fiber reinforced concrete,” Cement and Concrete Research, vol. 40, no. 6, pp. 956–965, Jun.
2010.
[63] V. C. Li, C. M. Chan, and C. K. Y. Leung, “Experimental determination of the tensionsoftening relations for cementitious composites,” Cement and Concrete Research, vol. 17,
no. 3, pp. 441–452, 1987.
[64] B. F. Sørensen and T. K. Jacobsen, “Large-scale bridging in composites: R-curves and
bridging laws,” Composites Part A: Applied Science and Manufacturing, vol. 29, no. 11, pp.
1443–1451, Nov. 1998.
[65] P. F. Fuchs and Z. Major, “Experimental Determination of Cohesive Zone Models for
Epoxy Composites,” Experimental Mechanics, vol. 51, no. 5, pp. 779–786, May 2010.
165

[66] M. Elices, G. V. Guinea, J. Gomez, and J. Planas, “The cohesive zone model: advantages,
limitations and challenges,” Engineering fracture mechanics, vol. 69, no. 2, pp. 137–163,
2002.
[67] H. Tan, C. Liu, Y. Huang, and P. H. Geubelle, “The cohesive law for the particle/matrix
interfaces in high explosives,” Journal of the Mechanics and Physics of Solids, vol. 53, no.
8, pp. 1892–1917, 2005.
[68] J. Abanto-Bueno and J. Lambros, “Experimental determination of cohesive failure
properties of a photodegradable copolymer,” Experimental mechanics, vol. 45, no. 2, pp.
144–152, 2005.
[69] I. Mohammed and K. M. Liechti, “Cohesive zone modeling of crack nucleation at
bimaterial corners,” Journal of the Mechanics and Physics of Solids, vol. 48, no. 4, pp.
735–764, 2000.
[70] Z. K. Guo, A. S. Kobayashi, J. C. Hay, and K. W. White, “Fracture process zone modeling
of monolithic Al2O3,” Engineering fracture mechanics, vol. 63, no. 2, pp. 115–129, 1999.
[71] B. Shen, I. Stanciulescu, and G. H. Paulino, “Inverse computation of cohesive fracture
properties from displacement fields,” Inverse Problems in Science and Engineering, vol. 18,
no. 8, pp. 1103–1128, Dec. 2010.
[72] B. Shen and G. H. Paulino, “Direct Extraction of Cohesive Fracture Properties from Digital
Image Correlation: A Hybrid Inverse Technique,” Experimental Mechanics, vol. 51, no. 2,
pp. 143–163, Mar. 2010.
[73] B. Shen and G. H. Paulino, “Identification of cohesive zone model and elastic parameters of
fiber-reinforced cementitious composites using digital image correlation and a hybrid
inverse technique,” Cement and Concrete Composites, vol. 33, no. 5, pp. 572–585, May
2011.
[74] A. L. Gain, J. Carroll, G. H. Paulino, and J. Lambros, “A hybrid experimental/numerical
technique to extract cohesive fracture properties for mode-I fracture of quasi-brittle
materials,” International Journal of Fracture, vol. 169, no. 2, pp. 113–131, May 2011.
[75] S. Hong and K. S. Kim, “Extraction of cohesive-zone laws from elastic far-fields of a
cohesive crack tip: a field projection method,” Journal of the Mechanics and Physics of
Solids, vol. 51, no. 7, pp. 1267–1286, 2003.
[76] S. Hong, H. B. Chew, and K.-S. Kim, “Cohesive-zone laws for void growth — I.
Experimental field projection of crack-tip crazing in glassy polymers,” Journal of the
Mechanics and Physics of Solids, vol. 57, no. 8, pp. 1357–1373, Aug. 2009.
[77] H. B. Chew, S. Hong, and K. S. Kim, “Cohesive zone laws for void growth—II. Numerical
field projection of elasto-plastic fracture processes with vapor pressure,” Journal of the
Mechanics and Physics of Solids, vol. 57, no. 8, pp. 1374–1390, 2009.

166

[78] S. T. Choi and K.-S. Kim, “Nanoscale planar field projections of atomic decohesion and
slip in crystalline solids. Part I. A crack-tip cohesive zone,” Philosophical Magazine, vol.
87, no. 12, pp. 1889–1919, Apr. 2007.
[79] W. H. Peters and W. F. Ranson, “Digital Imaging Techniques In Experimental Stress
Analysis,” Opt. Eng, vol. 21, no. 3, pp. 213427–213427–, 1982.
[80] W. J. W. M. A. Sutton, “Determination of displacements using an improved digital
correlation method.,” Image Vision Comput., vol. 1, pp. 133–139, 1983.
[81] H. A. Bruck, S. R. McNeill, M. A. Sutton, and W. H. P. Iii, “Digital image correlation using
Newton-Raphson method of partial differential correction,” Experimental Mechanics, vol.
29, no. 3, pp. 261–267, Sep. 1989.
[82] M. A. Sutton, S. R. McNeill, J. D. Helm, and Y. J. Chao, “Advances in two-dimensional
and three-dimensional computer vision,” in Photomechanics, Springer, 2000, pp. 323–372.
[83] J. Brillaud and F. Lagattu, “Limits and possibilities of laser speckle and white-light imagecorrelation methods: theory and experiments,” Appl Opt, vol. 41, no. 31, pp. 6603–6613,
Nov. 2002.
[84] Z. Sun, J. S. Lyons, and S. R. McNeill, “Measuring Microscopic Deformations with Digital
Image Correlation,” Optics and Lasers in Engineering, vol. 27, no. 4, pp. 409–428, 1997.
[85] M. Luo, “Displacement/strain measurements using an optical microscope and digital image
correlation,” Optical Engineering, vol. 45, no. 3, p. 033605, Mar. 2006.
[86] D. Lei, F. Hou, and X. Gong, “Investigation of Deformation at the Grain Scale in
Polycrystalline Materials by Coupling Digital Image Correlation and Digital Microscopy,”
Experimental Techniques, vol. 36, no. 2, pp. 24–31, 2012.
[87] H. Wang, H. Xie, Y. Ju, and Q. Duan, “Error Analysis of Digital Speckle Correlation
Method Under Scanning Electron Microscope,” Experimental Techniques, vol. 30, no. 2, pp.
42–45, 2006.
[88] M. A. Sutton, N. Li, D. C. Joy, A. P. Reynolds, and X. Li, “Scanning Electron Microscopy
for Quantitative Small and Large Deformation Measurements Part I: SEM Imaging at
Magnifications from 200 to 10,000,” Exp Mech, vol. 47, no. 6, pp. 775–787, Dec. 2007.
[89] M. A. Sutton, N. Li, D. Garcia, N. Cornille, J. J. Orteu, S. R. McNeill, H. W. Schreier, X.
Li, and A. P. Reynolds, “Scanning Electron Microscopy for Quantitative Small and Large
Deformation Measurements Part II: Experimental Validation for Magnifications from 200
to 10,000,” Exp Mech, vol. 47, no. 6, pp. 789–804, Dec. 2007.
[90] S. Chang, C. S. Wang, C. Y. Xiong, and J. Fang, “Nanoscale in-plane displacement
evaluation by AFM scanning and digital image correlation processing,” Nanotechnology,
vol. 16, no. 4, p. 344, Apr. 2005.

167

[91] Y. Sun and J. H. L. Pang, “AFM image reconstruction for deformation measurements by
digital image correlation,” Nanotechnology, vol. 17, no. 4, p. 933, Feb. 2006.
[92] A. Giachetti, “Matching techniques to compute image motion,” Image and Vision
Computing, vol. 18, no. 3, pp. 247–260, 2000.
[93] W. Tong, “An Evaluation of Digital Image Correlation Criteria for Strain Mapping
Applications,” Strain, vol. 41, no. 4, pp. 167–175, 2005.
[94] B. Pan, H. Xie, Z. Guo, and T. Hua, “Full-field strain measurement using a twodimensional Savitzky-Golay digital differentiator in digital image correlation,” Opt. Eng,
vol. 46, no. 3, pp. 033601–033601–10, 2007.
[95] H. W. Schreier and M. A. Sutton, “Systematic errors in digital image correlation due to
undermatched subset shape functions,” Experimental Mechanics, vol. 42, no. 3, pp. 303–
310, Sep. 2002.
[96] H. Lu and P. D. Cary, “Deformation measurements by digital image correlation:
Implementation of a second-order displacement gradient,” Experimental Mechanics, vol. 40,
no. 4, pp. 393–400, Dec. 2000.
[97] G. Vendroux and W. G. Knauss, “Submicron deformation field measurements: Part 2.
Improved digital image correlation,” Experimental Mechanics, vol. 38, no. 2, pp. 86–92,
Jun. 1998.
[98] H. W. Schreier, J. R. Braasch, and M. A. Sutton, “Systematic errors in digital image
correlation caused by intensity interpolation,” Opt. Eng, vol. 39, no. 11, pp. 2915–2921,
2000.
[99] Z.-F. Zhang, Y.-L. Kang, H.-W. Wang, Q.-H. Qin, Y. Qiu, and X.-Q. Li, “A novel coarsefine search scheme for digital image correlation method,” Measurement, vol. 39, no. 8, pp.
710–718, 2006.
[100] D. J. Chen, F. P. Chiang, Y. S. Tan, and H. S. Don, “Digital speckle-displacement
measurement using a complex spectrum method,” Appl. Opt., vol. 32, no. 11, pp. 1839–
1849, 1993.
[101] M. A. Sutton, S. R. McNeill, J. Jang, and M. Babai, “Effects Of Subpixel Image
Restoration On Digital Correlation Error Estimates,” Opt. Eng, vol. 27, no. 10, pp. 271070–
271070–, 1988.
[102] P.-C. Hung and A. S. Voloshin, “In-plane strain measurement by digital image
correlation,” Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol.
25, no. 3, pp. 215–221, Sep. 2003.
[103] J. D. Eshelby, W. T. Read, and W. Shockley, “Anisotropic elasticity with applications to
dislocation theory,” Acta Metallurgica, vol. 1, no. 3, pp. 251–259, May 1953.

168

[104] A. N. Stroh, “Dislocations and Cracks in Anisotropic Elasticity,” Philosophical Magazine,
vol. 3, no. 30, pp. 625–646, 1958.
[105] S. G. Lekhnitskii, Theory of Elasticity of an Anisotropic Elastic Body, First Edition.
Holden-Day, 1963.
[106] Z. Suo, “Singularities, Interfaces and Cracks in Dissimilar Anisotropic Media,” Proc. R.
Soc. Lond. A, vol. 427, no. 1873, pp. 331–358, Feb. 1990.
[107] J. K. Knowles and E. Sternberg, “On a class of conservation laws in linearized and finite
elastostatics,” Archive for Rational Mechanics and Analysis, vol. 44, no. 3, pp. 187–211,
1972.
[108] J. R. Rice and B. Budiansky, “Conservation laws and energy-release rates,” J. of Appl.
Mech, vol. 40, pp. 201–203, 1973.
[109] J. R. Rice, “Mathematical analysis in the mechanics of fracture,” Fracture: an advanced
treatise, vol. 2, pp. 191–311, 1968.
[110] T. C. T. Ting, Anisotropic Elasticity: Theory and Applications. Oxford University Press,
USA, 1996.
[111] K. Levenberg, “A method for the solution of certain non-linear problems in least squares,”
Quarterly Journal of Applied Mathmatics, vol. 2, pp. 164–168, 1944.
[112] D. Marquardt, “An Algorithm for Least-Squares Estimation of Nonlinear Parameters,”
Journal of the Society for Industrial and Applied Mathematics, vol. 11, no. 2, pp. 431–441,
1963.
[113] C. F. Shih, B. Moran, and T. Nakamura, “Energy release rate along a three-dimensional
crack front in a thermally stressed body,” International Journal of Fracture, vol. 30, no. 2,
pp. 79–102, 1986.
[114] ASTM D3039, “Test Method for Tensile Properties of Polymer Matrix Composite
Materials,” ASTM International, 2008.
[115] ASTM D3518, “Test Method for In-Plane Shear Response of Polymer Matrix Composite
Materials by Tensile Test of a 45 Laminate,” ASTM International, 2013.
[116] C. C. Chamis and J. H. Sinclair, “Ten-deg off-axis test for shear properties in fiber
composites,” Experimental Mechanics, vol. 17, no. 9, pp. 339–346, Sep. 1977.
[117] D. J. DeRenzo, Advanced Composite Materials Products and Manufacturers, 1st ed.
William Andrew, 1989.
[118] R. S. Piascik, J. C. Newman Jr, and J. H. Underwood, “The extended compact tension
specimen,” Fatigue & Fracture of Engineering Materials & Structures, vol. 20, no. 4, pp.
559–563, Apr. 1997.
169

[119] ASTM Standard E1922, “Standard test method for translaminar fracture toughness of
laminated polymer matrix composites.” ASTM International, West Conshohocken,PA,
2004.
[120] F. Lanza di Scalea, S. Hong, and G. Cloud, “Whole-field strain measurement in a pinloaded plate by electronic speckle pattern interferometry and the finite element method,”
Experimental Mechanics, vol. 38, no. 1, pp. 55–60, 1998.
[121] S. T. Pinho, P. Robinson, and L. Iannucci, “Fracture toughness of the tensile and
compressive fibre failure modes in laminated composites,” Composites science and
technology, vol. 66, no. 13, pp. 2069–2079, 2006.
[122] C. Dávila, C. Rose, and P. Camanho, “A procedure for superposing linear cohesive laws
to represent multiple damage mechanisms in the fracture of composites,” International
Journal of Fracture, vol. 158, no. 2, pp. 211–223, 2009.
[123] J. Réthoré, S. Roux, and F. Hild, “Noise-robust stress intensity factor determination from
kinematic field measurements,” Engineering Fracture Mechanics, vol. 75, no. 13, pp.
3763–3781, Sep. 2008.
[124] S. Morel, C. Lespine, J.-L. Coureau, J. Planas, and N. Dourado, “Bilinear softening
parameters and equivalent LEFM R-curve in quasibrittle failure,” International Journal of
Solids and Structures, vol. 47, no. 6, pp. 837–850, Mar. 2010.
[125] J. A. Nelder and R. Mead, “A Simplex Method for Function Minimization,” The
Computer Journal, vol. 7, no. 4, pp. 308–313, Jan. 1965.
[126] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence Properties of
the Nelder-Mead Simplex Method in Low Dimensions,” SIAM Journal of Optimization, vol.
9, pp. 112–147, 1998.
[127] Z. P. Bažant, “Concrete fracture models: testing and practice,” Engineering Fracture
Mechanics, vol. 69, no. 2, pp. 165–205, Jan. 2002.
[128] J. Planas, M. Elices, G. Guinea, F. Gómez, D. Cendón, and I. Arbilla, “Generalizations
and specializations of cohesive crack models,” Engineering Fracture Mechanics, vol. 70,
no. 14, pp. 1759–1776, Sep. 2003.
[129] ASTM E1820, “Test Method for Measurement of Fracture Toughness,” ASTM
International, 2013.
[130] C. G. D. Pedro P. Camanho, “Mixed-Mode Decohesion Finite Elements for the
Simulation of Delamination in Composite Materials,” 2002.
[131] N. I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity,
Softcover reprint of hardcover 1st ed. 1977 edition. Springer, 2010.

170