JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

44, 1, pp. 19-30, Warsaw 2006

ALGORITHM FOR DETERMINATION OF σ̃ij(n,θ), ε̃ij(n,θ),
ũi(n,θ), dn(n), In(n) FUNCTIONS IN

HUTCHINSON-RICE-ROSENGREN SOLUTION AND ITS 3D

GENERALIZATION

Jarosław Gałkiewicz

Marcin Graba

Faculty of Mechatronics and Machine Design, Kielce University of Technology

e-mail: jgalka@eden.tu.kielce.pl; mgraba@eden.tu.kielce.pl

In the paper the algorithm to determine σ̃ij(n,θ), ε̃ij(n,θ), ũi(n,θ),
dn(n), In(n) functions that are necessary to obtain values of stresses,
strains and displacements in the crack tip neighborhood according to
the Hutchinson, Rice and Rosengren solutions is presented. The algori-
thm can also be used to determine σ̃ij(n,θ,Tz), ε̃ij(n,θ,Tz), ũi(n,θ,Tz),
dn(n,Tz), In(n,Tz) functions for a 3Dapproximate solution of the stress
field in front of the crack introduced by Guo Wanlin, where constraint
due to the thickness effect is introduced through the Tz function.

Key words: fracturemechanics, HRRfields, 3DHRRfield approximation

1. Introduction

In 1968,Hutchinson (1968) publisheda fundamentalwork,which characte-
rised stress field in front of a crack for the non-linear Ramberg-Osgood (R-O)
material. Following Williams (1952), Hutchinson proposed the Airy function
for non-linear materials in the form of a series

φ= rsφ̃1(θ)+r
tφ̃2(θ)+ . . . (1.1)

where r and θ are polar coordinates of the coordinate system located at the
crack tip. The functions φ̃i(θ) describe angular changes of components of the
stress tensor.



20 J.Gałkiewicz, M.Graba

Hutchinson limited his considerations to the first dominant element of this
series. Using the compatibility equation and the R-O relationship, Hutchison
obtained formula for the stress field in front of a crack in the form

σe = K̃r
s−2σ̃e(θ,s) σθ = K̃r

s−2σ̃θ(θ,s)

σr = K̃r
s−2σ̃r(θ,s) σrθ = K̃r

s−2σ̃rθ(θ,s)
(1.2)

where s=(2n+1)/(n+1), n is theR-O exponent, σe is the equivalent stress,
σr, σθ, σrθ are the stress tensor components in the polar coordinate system,
K̃ is the plastic stress intensity factor, which can be related to the J-integral
through the relationship (McClintock, 1971)

K̃ =
( J
ασ0ε0In

) 1
1+n

(1.3)

where: α is the R-O constant, E is Young’s modulus, σ0 is the yield stress,
ε0 is the strain related to σ0 through the relation ε0 =σ0/E.
Thus, relationships (1.2) are usually known in the form

σij =σ0
( J
ασ0ε0Inr

) 1
1+n

σ̃ij(θ,n)+ . . .

εij =αε0
( J
ασ0ε0Inr

) n
1+n

ε̃ij(θ,n)+ . . . (1.4)

ui− ûi =ασ0r
( J
ασ0ε0Inr

)
n

1+n ũi(θ,n)

Functions σ̃ij(n,θ), ε̃ij(n,θ), ũi(n,θ), In(n) must be found by solving the fo-
urth order non-linear homogenous differential equation for the plane stress
and plane strain independently (Hutchinson, 1968). In the literature, these
functions are presented for limited values of the strain hardening exponent n.
However, it is very often required to use values of these functions for other
values of n than given in the literature. The program proposed in this pa-
per allows one to obtain all functions in the HRR solution for an arbitrary
exponent n.

2. Generalization of the HRR solution to the 3D-case

Values of the functions σ̃ij, ε̃ij, ũi and In for a predetermined stain harde-
ning exponent n depend on selection of the plane stress or plane strainmodel



Algorithm for determination of... 21

of an element. In a real specimen, the plane strain or stress states can be found
in the vicinity of the symmetry axis of the specimen close to the crack edge or
near the free surface of the specimen, respectively. The remaining part of the
specimen along the crack front is dominated by three-dimensional stress and
strain fields.
Guo (1993a) defined the Tz-function as

Tz =
σ33

σ11+σ22
(2.1)

For the non-linear plastic materials, Tz is equal to 0 for plane stress and 0.5
for plane strain. Thus, Tz changes from 0.5 to 0 along the crack front from
the specimen axis to the specimen surface.
Using function (2.1), Guo postulated the Maxwell stress function in the

form
φi = K̃r

s(Tz)φ̃i(θ,Tz) (2.2)

where the functions φ̃i(θ,Tz) describe the angular changes of stress tensor
components. In the functions φi both the s exponent and φ̃i functions were
assumed to be dependent on the function Tz. The φi functions were used to
obtain a solution analogous to the HRR field. However, Guo’s solution is not
limited to the plane strain or plane stress cases. Since the Tz function changes
along the crack front, Guo’s solution covers also these layers of the material
along the crack front which are in the 3D state of the stress and strain field.
The only requirement is to know the Tz(x1,x3) function.
Guo showed that the s(Tz) function is equal to s in the HRR solution

for the plane strain (Tz = 0.5) or plane stress (Tz = 0) cases only. However,
between the specimen axis and the specimen surface s(0<Tz < 0.5) 6= sHRR.
Moreover, for all these points along the crack front, the Jintegral is not path
independent.
Because at the above problems, Guo postulated a general approximate

formula for a quasi-three-dimensional case in the form

σij =σ0
( Jfar
ασ0ε0In(n,Tz)r

) 1
1+n

σ̃ij(θ,n,Tz) (2.3)

where: Jfar is the J-integral computed for the contour of integration drawn
over the domain dominated by the plane stress (r >B where B is specimen
thickness). In Eq. (2.3), the out of plane effect was taken into account by the
value of Tz. Guo demonstrated, by comparison with numerical results, that
using Eq. (2.3) the error in stress values along the crack edge was not greater
than 7% for n = 10 and decreased with n. Equation (2.3) can be used for



22 J.Gałkiewicz, M.Graba

each point along the crack front as well as for the mean value of Tz through
the specimen thickness. In that case, one obtains an intermediate case with
respect to the plane strain and plane stress models. In the polar coordinate
system, the Guo solution is

σr =Kr
s−2σ̃r(n,θ,Tz) σθ =Kr

s−2σ̃θ(n,θ,Tz)
(2.4)

σrθ =Kr
s−2σ̃rθ(n,θ,Tz)

where

σ̃r(n,θ,Tz)= sφ̃+
∂2φ̃

∂θ2
σ̃θ(n,θ,Tz)= s(s−1)φ̃

(2.5)

σ̃rθ(n,θ,Tz)= (1−s)
∂φ̃

∂θ

and

σe =Kr
s−2σ̃e(n,θ,Tz)

(2.6)

σz =TzKr
s−2[σ̃r(n,θ,Tz)+ σ̃θ(n,θ,Tz)]

In the following section the functions in Eqs (2.5) will be determined.

3. Calculation of values of σ̃ij(n,θ,Tz), ε̃ij(n,θ,Tz), ũi(n,θ,Tz),
dn(n,Tz), In(n,Tz)

The constitutive relation of a homogeneous isotropic elastoplastic continu-
um can be expressed by

εij =(1+ν)Sij +
1−2ν

3
σkkδij +

3

2
ασn−1e Sij (3.1)

where Sij is the stress deviator, δij is the Kronecker delta, ν is the Poisson
ratio and σe is the von Misses equivalent stress.
The strain components in the cylindrical coordinate system can bewritten

in the form

εrr =(1+ν)σrr− (1+Tz)ν(σrr+σθθ)+
3

2
ασn−1eff

[
σrr−

1+Tz
2
(σrr +σθθ)

]

εθθ =(1+ν)σθθ− (1+Tz)ν(σrr +σθθ)+
3

2
ασn−1eff

[
σθθ−

1+Tz
2
(σrr+σθθ)

]

εrθ =(1+ν)σrθ+
3

2
ασn−1eff σrθ (3.2)



Algorithm for determination of... 23

Making use of Eq. (2.4) and Eq. (3.2) in the compatibility equation (3.3)

1

r

∂2

∂r2
(rεθθ)+

1

r2
∂2εrr
∂θ2
−
1

r

∂εrr
∂r
−
2

r2
∂

∂r

(
r
∂εrθ
∂θ

)
=0 (3.3)

one can obtain forth order non-linear homogeneous differential equation (3.4)
with twounknowns: the functions φ̃(θ,Tz) andthe singularityparameter s(Tz)

[n(s−2)+1][n(s−2)]σ̃n−1e

[
φ̃
(
s
(3(s−1)−s(1+Tz)

2

))
−
∂2φ̃

∂θ2

(1+Tz
2

)]
+

+
∂2

∂θ2

{
σ̃n−1e

[
φ̃
(
s
(3−s(1+Tz)

2

))
+
∂2φ̃

∂θ2

(2−Tz
2

)]}
+

(3.4)

−n(s−2)σ̃n−1e

[
φ̃
(
s
(3−s(1+Tz)

2

))
+
∂2φ̃

∂θ2

(2−Tz
2

)]
+

−2[n(s−2)+1]
∂

∂θ

{3
2
(1−s)σ̃n−1e

∂φ̃

∂θ

}
=0

In order to find values of σ̃ij(n,θ,Tz), ε̃ij(n,θ,Tz), ũi(n,θ,Tz), dn(n,Tz),
In(n,Tz), Eq. (3.4) should be solved first. After transformation of Eq. (3.4)
one can obtain

∂4φ̃

∂θ4
=
NUMERATOR

DENOMINATOR
(3.5)

where

NUMERATOR=

−
[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+Z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−1
2
[
Z9φ̃−Z10

∂2φ̃

∂θ2

]
+

−
n−1

2

[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+Z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−3
2

·

·
[∂φ̃
∂θ

(
2Z16φ̃+

∂2φ̃

∂θ2
(2Z15+Z17)

)
+
∂3φ̃

∂θ3

(
Z17φ̃+2Z12

∂2φ̃

∂θ2

)][
Z11
∂φ̃

∂θ
+2Z5

∂3φ̃

∂3θ

]
+

−
n−1

2

n−3

2

[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−5
2

· (3.6)

·
[∂φ̃
∂θ

(
2Z16φ̃+

∂2φ̃

∂θ2
(2Z15+Z17)

)
+
∂3φ̃

∂θ3

(
Z17φ̃+2Z12

∂2φ̃

∂θ2

)]2(
Z4φ̃+Z5

∂2φ̃

∂θ2

)
+

−
n−1

2

[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+Z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−3
2

·



24 J.Gałkiewicz, M.Graba

·
[∂φ̃
∂θ

(
2Z16
∂φ̃

∂θ
+
∂3φ̃

∂θ3
(2Z15+2Z17)

)
+
∂2φ̃

∂θ2

(
2Z16φ̃+

∂2φ̃

∂θ2
(2Z15+Z17)

)
+

+2Z12
(∂3φ̃
∂θ3

)2](
Z4φ̃+Z5

∂2φ̃

∂θ2

)

DENOMINATOR=Z5
[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−1
2

+

+
n−1

2

[
Z16φ̃

2+Z15
(∂φ̃
∂θ

)2
+
(
Z17φ̃+Z12

∂2φ̃

∂θ2

)∂2φ̃
∂θ2

]n−3
2

· (3.7)

·
(
Z17φ̃+2Z12

∂2φ̃

∂θ2

)(
Z4φ̃+Z5

∂2φ̃

∂θ2

)

and

Z1 = [n(s−2)+1][n(s−2)] Z10 =Z1Z3+Z8Z5−Z4+Z6Z7

Z2 = s
(3(s−1)−s(1+Tz)

2

)
Z16 = s

2Z12+Z12Z13−sZ14

Z3 =
1+Tz
2

Z12 =1−Tz +T
2
z

Z4 = s
(3−s(1+Tz)

2

)
Z13 = s

2(s−1)2

Z5 =
2−Tz
2

Z14 =(1+2Tz −2T
2
z )s(s−1)

Z6 =2[n(s−2)+1] Z15 =3(s−1)
2

Z7 =
3

2
(1−s) Z11 =2Z4−Z6Z7

Z8 =n(s−2) Z17 =2sZ12−Z14

Z9 =Z1Z2−Z8Z4

To solve Eq. (3.5), a combination of numerical methods was used (see
Fig.1). The function φ̃(θ) was found by the forth order Runge-Kutta method
(Burden andFaires, 1985), and in order to find the initial value of ∂2φ̃(0)/∂θ2

and the value of the power exponent s(Tz), the shootingmethod (Burden and
Faires, 1985; Marciniak et al., 2000) was used.

It turns out that during theprocess of numerical computations,more itera-
tions are requiredwhen the strain hardening exponent n increases. A solution
to the plane strainproblem ismore easily obtained.A satisfactory convergence



Algorithm for determination of... 25

Fig. 1. The algorithm for numerical solution to Eq. (3.2)



26 J.Gałkiewicz, M.Graba

was found for the strain hardening exponent n¬ 20 in the case of the plane
stress, and for n ¬ 30 in the plane strain. The convergence depends on Tz
and theworst situationwas observed for Tz ∈ [0.23−0.27]. Exemplary results
for the strain hardening exponent n= 20 and the plane stress are presented
in Fig.2.

Fig. 2. Exemplary results for n=3, Tz =0; (a) stress functions σ̃(θ), (b) strain
functions ε̃(θ), (c) displacement functions ũ(θ) and their first derivatives,

(d) function φ(θ) and its derivatives

4. Calculations of dn and In functions

In order to compute the crack opening displacement (Neimitz, 1998), one
must know the value of dn(α,ε,n,Tz)

δT = dn(α,ε,n,Tz)
J

σ0
(4.1)



Algorithm for determination of... 27

Using Eq. (2.1) and themethod proposed by Shih (1981), (see Fig.3), one can
compute dn for anymaterial from the following formula

dn =
2

In
ũ2(π,n,Tz)

(ασ0
E
[ũ1(π,n,Tz)+ ũ2(π,n,Tz)]

)1
n

(4.2)

Fig. 3. Definition of CTOD

The value of the function In(n,Tz) follows form the path-independency of
the J-integral. The J-integral is path-independent when

I(Tz,n)=

π∫

−π

{ n
n+1

σ̃n+1e cosθ−
3

2

(
sinθ[σ̃r(ũθ− ũ

′

r)− σ̃rθ(ũr + ũ
′

θ)]
)
]+

(4.3)

+
3

2
cosθ[n(s−2)+1](σ̃rũr + σ̃rθũθ)

}
dθ

5. Comparison of results

Thevalues obtainedwith thehelp of theprogramhrr par.exe for In(n) and
(π/In)

1/(n+1) which depend on the R-O power exponent n, may be compared
with the Hutchinson results obtained for Tz = 0.5 or 0 (Hutchinson, 1968).
Differences are presented in Table 1. The differences are smaller than 0.3%
both for the plane stress and plane strain.

The values of the singularity exponent s(n,Tz) are close to Guo’s results
(Guo, 1993a,b). The differences are less than 0.15% almost for all cases (Ta-
ble 2). Only for n=8 and Tz =0.45 this difference is about 1.35%. In Fig.4,



28 J.Gałkiewicz, M.Graba

the function dn (obtained using hrr par.exe) is presented. Results are close to
those presented by Guo (1995).

Table 1.Values of In(n) function, (π/In)
1/(n+1)

n In(n) In(n)HRR Difference (π/In)
1

n+1 (π/In)
1

n+1

∣∣
HRR

Difference

Plane stress (Tz =0)

3 3.85 3.86 0.26% 0.950 0.949 0.15%

5 3.41 3.41 0.00% 0.986 0.987 0.06%

9 3.03 3.03 0.00% 1.004 1.004 0.04%

13 2.87 2.87 0.00% 1.006 1.006 0.05%

Plane strain (Tz =0.5)

3 5.51 5.51 0.00% 0.869 0.869 0.00%

5 5.02 5.01 0.20% 0.925 0.925 0.02%

9 4.60 4.60 0.00% 0.963 0.963 0.04%

13 4.40 4.40 0.00% 0.976 0.976 0.02%

Table 2.Values of singularity exponent s(n,Tz)

n
Tz

0 0.3 0.4 0.45 0.5

Our results

3 −0.2500000 −0.2380839 −0.2380157 −0.2426839 −0.2500000

3.6364 −0.2156929 −0.2031756 −0.2034906 −0.2084886 −0.2156854

8 −0.1111266 −0.0997175 −0.1019225 −0.1059254 −0.1111112

10 −0.0909101 −0.0803697 −0.0828057 −0.0863690 −0.0909083

Guo results

3 −0.250000 −0.237825 −0.237730 −0.242500 −0.250000

3.6364 −0.215686 −0.203186 −0.203186 −0.208336 −0.215686

8 −0.111111 −0.099617 −0.101835 −0.104511 −0.111111

10 −0.090909 −0.080280 −0.082909 −0.086409 −0.090909

Differences between our results andGuo ones (Shih, 1981)

3 0.00% 0.11% 0.12% 0.08% 0.00%

3.6364 0.00% 0.01% 0.15% 0.07% 0.00%

8 0.01% 0.10% 0.09% 1.35% 0.00%

10 0.00% 0.11% 0.12% 0.05% 0.00%

Theresultspresented in this section concernan infiniteplatewith the crack
loaded at infinity. The results for a finite body can be found inGałkiewicz and
Graba (2004).



Algorithm for determination of... 29

Fig. 4. Results of normalized crack tip opening displacement for different n vs. Tz

6. Conclusions

Using the program hrr par.exe, one can obtain results for a wide range of
materials which are characterised by n, σ0, ε0, E in an easy and fast way.
The obtained results are accurate when compared to those published in the
literature.

The programhrr par.exe – the source version (for Delphi 6) and the compiled
source are available on:

http://www.tu.kielce.pl/∼pgfm/HRR.htm

http://www.tu.kielce.pl/∼mgraba → Fracture

http://www.tu.kielce.pl/∼mgraba/Fracture 03.htm

Acknowledgements

We are pleased to acknowledge helpful discussions with prof. A. Neimitz from
Kielce University of Technology.

Thework presented in this paperwas carried outwith the support of Polish State

Committee for Scientific Research; grant No. 5T07C00425.

References

1. Burden R.L, Faires J.D., 1985, Numerical Analysis, third edition, PWS-
KENTPublishing Company, Boston

2. Gałkiewicz J., GrabaM., 2004,Aproksymowanie rozwiązań3Dprzed fron-
tem szczeliny poprzez wprowadzenie więzóww kierunku grubości,XX Sympo-



30 J.Gałkiewicz, M.Graba

zjum Zmęczenie i Mechanika Pękania, BydgoszczPieczyska, 65-71 (in Polish)

3. GuoW., 1993a,Elastoplastic three dimensional crackborder field – I. Singular
structure of the field, Engineering Fracture Mechanics, 46, 1, 93-104

4. GuoW., 1993b,Elastoplastic threedimensional crackborderfield – II.Asymp-
totic solution for the field, Engineering Fracture Mechanics, 46, 1, 105-113

5. GuoW., 1995,Elastoplastic three dimensional crackborder field - III. Asymp-
totic solution for the field, Engineering Fracture Mechanics, 51, 1, 51-71

6. Hutchinson J.W., 1968, Singular behaviour at the end of a tensile crack in a
hardeningmaterial, Journal of the Mechanics and Physics of Solids, 16, 13-31

7. Marciniak A., Gregulec D., Kaczmarek J., 2000,Podstawowe procedury
numeryczne w języku Turbo Pascal, WydawnictwoNakom, Poznań (in Polish).

8. McClintock F.A., 1971, Plasticity aspects of fracture, In:Fracture – an Ad-
vanced Treatise, H. Liebowitz (edit.)

9. NeimitzA., 1998,Mechanika pękania,WydawnictwoNaukowePWN,Warsza-
wa (in Polish)

10. ShihC.F., 1981,Relationbetween the J-integral and the crackopeningdispla-
cement for stationary and extending cracks, Journal of Mechanics and Physics
of Solids, 29, 305-329

11. ShihC.F., 1983,TablesofHutchinson-Rise-Rosengrensingularfieldquantities,
BrownUniversity Report, MRLE-147

12. WilliamsM.L., 1952,Stress singularities resulting fromvariousboundarycon-
ditions in angular corners of plates in extension, Journal of Applied Mechanics,
19, 526-528, 3, 47-225

Algorytm wyznaczania funkcji σ̃ij(n,θ), ε̃ij(n,θ), ũi(n,θ), dn(n), In(n)
w rozwiązaniu HRR i jego trójwymiarowym uogólnieniu

Streszczenie

Wartykule zaprezentowanoalgorytmpozwalającyna określenie funkcji σ̃ij(n,θ),
ε̃ij(n,θ), ũi(n,θ), dn(n), In(n) niezbędnych do opisu pola naprężeń, odkształceń
i przemieszczeńwmateriałach nieliniowychwedług prawaRamberga-Osgooda.Algo-
rytm pozwala uzyskać również wartości funkcji σ̃ij(n,θ,Tz), ε̃ij(n,θ,Tz), ũi(n,θ,Tz)
dla próbek trójwymiarowychprzywykorzystaniu parametru Tz, który zależy od gru-
bości próbki i może być wyznaczony numerycznie.

Manuscript received May 18, 2005; accepted for print December 17, 2005