Direct modelling technique is used to study stresses at the tip of a mode I crack where domain of the material microstructure is explicitly discretised and represented by finite elements. Two cases of internal (at crack faces) and external uniform loading are considered. A multiscale asymptotic model accounting for both the external boundaries and the non-singular stresses at the crack tip are used to analyse the numerical results. It is shown that the parameters (coefficients) of the expansion of the stress concentration at the crack tip can be recovered for both cases of loading from the simulated stress distribution ahead of the crack tip. Applicability of this technique is validated by recovering the parameters for a wider range of heterogeneous material properties

Dyskin, A. V.

Pant, S. R.

University of Queensland eSpace

SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
Multiscale Modelling Of The Stress Singularity At A Mode I Crack Tip 
S.R. Pant and A.V. Dyskin 
School of Civil and Resource Engineering 
University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia 
 
 
 
ABSTRACT: Direct modelling technique is used to study stresses at the tip of a mode I crack where domain 
of the material microstructure is explicitly discretised and represented by finite elements. Two cases of 
internal (at crack faces) and external uniform loading are considered. A multiscale asymptotic model 
accounting for both the external boundaries and the non-singular stresses at the crack tip are used to analyse 
the numerical results. It is shown that the parameters (coefficients) of the expansion of the stress 
concentration at the crack tip can be recovered for both cases of loading from the simulated stress 
distribution ahead of the crack tip. Applicability of this technique is validated by recovering the parameters 
for a wider range of heterogeneous material properties. 
 
1 INTRODUCTION 
Analysis of fracture phenomena controlled by complex material microstructure requires, at least as 
a first step, explicit numerical simulation of fracture propagation. Particularly attractive is a 
simulation method in which the microstructure is represented by numerical elements (e.g. finite 
element and finite difference methods), in which fractures are treated as sequences of soft elements 
with largely reduced or zero moduli. However, results of such simulations are difficult to interpret 
in terms of conventional theory of the LEFM. The applicability of this type of simulations as well as 
the LEFM interpretation of the results is considered in this paper for the case of mode I crack. The 
simulations are conducted by the finite difference code FLAC and the results of the analyses are 
compared with a multiscale asymptotic model developed to account for both external boundaries 
and the non-singular stresses at the crack tip.  
The analysis involves three scales: (1) the scale of the modelling domain, which in this particular 
case is approximated by an infinite strip; at this scale the crack is seen as a point defect. (2) the 
scale of the crack at which the process zone is seen as a point; at this scale the singular term and 
two non-singular terms are determined. (3) the scale of the process zone; at this scale the 
interpretation of the simulation results and comparison with the analytical model are conducted.  
Two cases of loading are considered, which produce different non-singular parts of the stress 
concentration – external uniform loading and uniform internal loading at crack faces. The 
macrostructural property of the material is modelled by randomly varying Young’s moduli. Stresses 
are calculated numerically at different distances ahead of the crack tip along its axis and the 
coefficients of the expansion of the stress concentration at the crack tip are recovered and compared 
with an analytical asymptotic solution. 
2 DIRECT FLAC SIMULATION OF STRESS FIELD IN MODE I CRACK 
PROPAGATION  
Figure 1 shows a finite difference mesh of 200 x 100 elements is used in this plane stress analysis. 
Crack is represented as a straight slot-like opening of one finite element width. This was achieved 
by effective removal of elements by setting it to null elements where all the stress tensors and forces 
in these elements become zero.  
Two types of loadings are considered: remote loading σyy=p, and tractions applied at the crack 
faces, ty=p. Other stress components at the sample boundary and the crack faces are zero. Owing to 
the linearity, the particular calculations are conducted with p=1 MPa. The numerical model is 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
assumed to be in equilibrium condition when the FLAC out-of-balance force is small compared to 
the total applied external force viz within an error of 0.03%.  
X
Y
σyy
σyy = 1 MPa
400 mm
150 mm
 
Figure 1. Typical finite difference mesh, loading and FLAC-element representation of crack used in the simulation. 
This type of simulation involves a natural macroscopic crack propagation criterion – a criterion 
of local failure under the stress in the element at the crack tip. The rational behind such a criterion is 
that the FLAC elements represent both structural and crack elements of the material. Therefore, 
microscopic structural details of the crack propagation cannot be simulated by such codes. 
In the numerical simulations the stress component σyy normal to the crack length was calculated 
for different types of loading at 16 points ahead of the crack tip (Figure 2) assuming that the 
material is homogeneous. In reality, rocks are heterogeneous on any scale, which inevitably 
influences fracture processes. This material heterogeneity effect was modelled by assigning 
uniformly distributed random variations of Young’s modulus in each FLAC element. The ranges of 
variation were 10, 20 and 40% (±5,  ±10 and ±20% of the mean value respectively). Figures 3 and 4 
show the stresses calculated along the crack axis (the x-axis) for the case of loading at crack faces 
for different ranges of variation of Young’s modulus. It is seen that the stress distribution is barely 
affected by these heterogeneity variations. This can be explained as follows. The lengths of the 
FLAC-elements at the crack tip represent the microstructure of the material in the numerical 
simulation. This length of the microstructural element is very small as compared to the crack length, 
i.e. d<<a, where d is the length of the microstructural element. The value of the distance from the 
crack tip, r, varies from d to 14d which is still very small in comparison to the crack length. 
Therefore, between these two scales a mesoscale associated with the dimension H of the 
representative volume element can be introduced, i.e. d << H << a. The stress field will only 
slightly differ for the one obtained from the modelling performed in an effective continuum medium 
of the scale H which corresponds to averaging of stress and strain fields over the volume elements 
of size H. This effective medium is characterised by uniform effective moduli replacing 
heterogeneous moduli of the original material. That is why the influence of the modulus 
heterogeneity is negligible everywhere except in the immediate vicinity of the crack tip where the 
effective medium modelling is no longer applicable. Correspondingly, Figures 3 and 4 show that the 
stresses near the crack tips are fluctuating as compared to the case of uniform modulus.  
3 ASYMPTOTICS OF NON-SINGULAR PART OF THE STRESS FIELD  
LEFM assumes that the process zone at the crack tip is negligibly small compared to the crack 
length or the concept of small-scale yielding [Rice, 1968]. This assumption allows reasonable 
approximation of stresses near the crack tip by considering only the singular part of stress 
Comment [AD1]: Sudeep, 
please explain 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
concentration. However, when the process zone length is not very small compared to the crack 
length, the contribution of non-singular terms becomes important. Therefore, it is convenient to 
continue usage of classical elastic crack solution, as in LEFM, but as a next step take into account 
the non-singular part of stress concentration [Dyskin, 1997]. 
The expression for normal stress at and ahead of the tip of a crack axis loaded by tractions –p(t) 
applied to the crack faces has the following form [Muskhelishvili, 1953] with the origin of the x-
axis being placed at the crack centre 
σ πy a
a
x
x a
p t a t
x t
dt( ) ( )= −
−
−−∫
1
2 2
2 2
 (1) 
where 2a is the crack length, r is the distance from the crack tip and –p(t) is the traction applied to 
the crack faces.  
Developing the full solution can be complicated for complex type of loading especially in the 
presence of finite boundaries. However, when the crack length is still large compared to the process 
zone length, equation (1) can be used to derive general asymptotic solution for stresses at the crack 
tip [Dyskin, 1997] in the following form 
( )(0) (1)( ) , ,    
2
I
y
Kr s s r o r r x a
r
σ π= + + + ∆ ∆ = = −  (2) 
where KI is the mode I stress intensity factor. 
( ) for internal uniform load
0 otherwise
O ps
−⎧= ⎨⎩  
(1)
3/ 2
3 for uniform load (external and internal)
4 2
5 for concentrated forces
2 (2 )
p
as
P
aπ
⎧⎪⎪= ⎨⎪−⎪⎩
 
0
1
2
3
4
5
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018
FLAC
Asymptotic fited
with FLAC
r (m)
σ
yy
(M
Pa
)
σ
yy
(M
Pa
)
0
1
2
3
4
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018
FLAC
Asymptotic fitted
with FLAC
r (m)
σ
yy
(M
Pa
)
σ
yy
(M
Pa
)
 
   (a)                  (b) 
Figure 2. Numerically calculated distribution of σyy along the line of crack continuation (data points) and the fitted 
asymptotics (solid line) for the case of: (a) remote loading; (b) loading at the crack faces.  
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
1
1.5
2
2.5
3
3.5
4
4.5
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
10% Variation in moduli
20% variation in moduli
40% variation in moduli
Constant Modulus
Asymptotic fitted with
FLAC (constant moduli)
σ
yy
(M
Pa
)
r (m)
σ
yy
(M
Pa
)
 
Figure 3. Numerically calculated distribution of σyy along the line of crack continuation (data points) and the fitted 
asymptotics (solid line) for the case of remote loading for heterogeneous materials. 
 
r (m)
σ
yy
(M
Pa
)
0
0.5
1
1.5
2
2.5
3
3.5
4
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016
10% Variation in moduli
20% variation in moduli
40% variation in moduli
Constant Modulus
Asymptotic fitted with
FLAC (constant moduli)
σ
yy
(M
Pa
)
 
Figure 4. Numerically calculated distribution of σyy along the line of crack continuation (data points) and the fitted 
asymptotics (solid line) for the case of loading at crack faces for heterogeneous materials. 
Now consider the interaction of this crack with the boundary of the sample. The complete 
influence of the boundary of a rectangular sample cannot be expressed in mathematically closed 
form. However, the sample geometry is such that the lateral boundaries are the closest to the crack. 
Thus, one can hope that by accounting only for these boundaries a sufficiently accurate model can 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
be obtained. This can be done by considering the crack situated in the middle of an infinite strip in x 
direction, Figure 5. 
x
2l
h p
h
y
 
Figure 5. The asymptotic model of crack in an infinite strip. 
The next simplification comes from the use of the dipole asymptotics method [e.g., Dyskin and 
Mühlhaus, 1995] to calculate the interaction between the crack and the boundary of the strip. This 
method takes into account only the main asymptotic terms characterising the interaction – the terms 
of the order of l2/h2, where 2h is the strip width. These terms represent the additional tractions on 
the crack faces reflecting the influence of the boundary.  
These additional tractions are found as follows. Firstly, the crack is considered in an infinite 
plane and stresses are calculated on the lines which correspond to the boundaries of the strip. These 
stresses are calculated with accuracy of the main asymptotic terns, i.e. only the terms of the order of 
l2/r2, where r is the distance from the crack centre. Then the corresponding components of these 
stresses, with opposite signs, are applied to the boundaries of the strip without the crack and the 
stress is calculated at the centre of the strip, i.e. at the place of the crack. Since these tractions are 
already of the order of l2/h2, their variations over the crack length are of the order l3/h3, and hence 
should be neglected within the accuracy limits of the method. This consideration allows one to 
assume the additional tractions to be uniform, which leads to an especially simple method of 
calculating these additional tractions. Then this uniform stress is applied to the crack faces and the 
additional stress intensity factor is calculated. This constitutes the correction to the stress intensity 
factor associated with the effect of the boundaries [see Dyskin et al., 2000 for details]. As a result, 
(2) becomes  
⎟⎟⎠
⎞
⎜⎜⎝
⎛ +=⎟⎟⎠
⎞
⎜⎜⎝
⎛ +=∆+++= ∗∗∗
∗
2
2
)1()1(
2
2
)1()0( 14.11,14.11,
2
)(
h
lss
h
lKKrss
r
Kr IIIy πσ  (3) 
where 2l is the crack length and 2h is the sample height. 
4 RECOVERING OF THE PARAMETERS  
Equation (3) gives asymptotic solution for the stress distribution at the crack tip. It is interesting to 
compare this result with the numerical solution to see whether the parameters K*I, s(0) and s*(1) can 
be recovered. Technically these parameters are the coefficients of multiple regressions on functions 
r-1/2, 1 and r1/2. The results of the fitting are shown in Figures 1-3 with the values of the parameters 
given in Tables 1-3. It is seen that for the case of homogeneous material the recovered parameters 
are very close to the actual values in both cases. The error of the recovery however increases with 
heterogeneity. It is also noted that average values and standard deviations of recovered parameters 
are higher for the case when cracks are loaded at its faces compared to the case where load is at 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
infinity for the same heterogeneity. This observed discrepancy can be attributed to the heterogeneity 
acting as local points of stress concentrators.  
Table 1. Comparison of recovered and actual parameters when load is applied at: (a) infinity; (b) crack faces. 
 
0.060.020s(0)
0.48273.4363.521s(1) 0.9994
0.00450.2880.292KI
Coefficient of 
determination 
Standard
errorRecoveredActual
Loaded at boundary
Parameters
      (a) 
0.0685-1.0006-1.0s(0)
0.5343.653.521s(1) 0.9994
0.004780.28960.292KI
Coefficient of 
determination 
Standard
errorRecoveredActual
Loaded at crack faces
Parameters
 
(b) 
Table 2. Comparison of recovered and actual parameters when load is applied at infinity with varied material properties.   
 
0.151-0.0570.078-0.0020.065-0.001S (0)
0.693.7680.4423.4310.393.457S (1)
0.0180.2980.0080.2920.0060.292K I
Standard 
D eviationA verage
Standard 
D eviationA verage
Standard 
D eviationA verage
40%  variation in m oduli20%  variation in m oduli10%  variation in m oduli
Loaded at boundaries
Param eter
 
Table 3. Comparison of recovered and actual parameters when load is applied at crack faces with varied material 
properties. 
0.08-1.1370.079-1.0610.045-0.983S(0)
0.493.9760.4173.9070.2673.366S(1)
0.0070.310.0060.2960.0040.29KI
Standard 
DeviationAverage
Standard 
DeviationAverage
Standard 
DeviationAverage
40% variation in moduli20% variation in moduli10% variation in moduli
Loaded at crack faces
Parameter
 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
5 CONCLUSIONS 
The problem of crack propagation in samples of heterogeneous materials is influenced by three 
scales: the scale of the sample, the scale of crack and the scale of heterogeneities. Direct modelling 
technique is used to study stresses at the tip of a mode I crack where the material microstructure is 
explicitly represented by FLAC-finite elements. Comparison of the results of the simulation with a 
multiscale asymptotic model shows that the multiscale concept offers an efficient method of 
handling problems featuring a number of distinctive scales. As a result, the parameters 
(coefficients) of the expansion of the stress concentration at the crack tip can be efficiently 
recovered. The recovery becomes more distinct for heterogeneous materials, although the error of 
recovery also increases with heterogeneity.  
ACKNOWLEDGEMENTS 
The authors acknowledge the financial support of the Australian Computational Earth Systems 
Simulator – Major National Research Facility and the 2001-2003 ARC Grant A00104937. 
REFERENCES 
Dyskin, A V and Mühlhaus, H -B: Equilibrium bifurcations in dipole asymptotic model of periodic crack 
arrays (Chapter 3), Continuum Models for Materials with Micro-Structure, eds. H -B Mühlhaus, John Wiley 
& Sons Ltd, pp. 69-104 (1995) 
Dyskin, A V, Crack growth criteria incorporating non-singular stresses: size effect in apparent fracture 
toughness, Int. J. Fracture, 83, pp. 191-206 (1997)  
Dyskin, A V, Germanovich, L N and Ustinov, K B, Asymptotic analysis of crack interaction with free 
boundary, Int. J. Solids & Structures, 37, pp. 857-886 (2000) 
Muskhekishvili, N I, Some Basic Problems of the Mathematical Theory of Elasticity, P. Noordhoff Ltd. 
Groningen-Holland (1953) 
Rice, J R, Mathematical analysis in the mechanics of fracture, In: Fracture, An Advanced Treatise, 
Mathematical Fundamentals, eds. H Liebowitz, Academic Press, 2, pp. 191-311 (1968) 
SIF2004 Structural Integrity and Fracture. http://eprint.uq.edu.au/archive/00000836 
 


Multiscale Modelling Of The Stress Singularity At A Mode I Crack Tip

http://espace.library.uq.edu.au/view/UQ:10090/pant_sif04.pdf

Multiscale Modelling Of The Stress Singularity At A Mode I Crack Tip

Abstract

Similar works

Full text

Available Versions

University of Queensland eSpace