. . .,.ERlAIS SCIEMCEa E_IMEERIMe ELSEVIER
A
Materials Science and Engineering A209 (1996) 82-90
Modelling of interfacial fracture V.L. Rabinovich, V.K. Sarin Manufacturing Engineering, Boston University, Boston, MA 02215, USA
Abstract This paper adresses the problem of modelling the complex non-linear phenomena of indentation and fracture of brittle materials. These include frictional contact interaction, crack propagation and debonding of laminated brittle three-dimensional finite bodies. An original unified methodology for the solution of the related boundary value problems, based on combination of the variational and boundary integral multi-region approaches, has been developed. The discretization of the involved surfaces and unknown functions leads to a finite dimensional quadratic programming problem, solved by a modification of the gradient-projection method. The discretized version of the algorithm has been implemented in the form of computer codes for both conventional serial and parallel Connection Machines CM2/5, multi-processor SGI Challenge Series. The model components were tested against analytical and numerical results. Crack propagation and debonding of the interfacial boundaries in a quasi-static scenario were investigated. The model has been used for numerical simulation of the Micro-Scratch Test to aid in the development of procedures to evaluate adhesion for coatings. The model demonstrated that for the Micro-Scratch Test (using the Knoop indenter) the dominant mechanism of the interfacial debonding is different than in indentation fracture of the homogeneous materials. Numerous experiments with the microscratch technique show that this test, with the aid of numerical model, can quantify a coating/substrate interfacial strength and characterize the factors influencing it. Keywords: Modelling; Interfacial fracture; Coatings
1. Introduction
A model for studying the micromechanics of the fracture initiation and propagation in and near interfacial regions in multi-layered bodies under arbitrary loading, including frictional contact loading has been developed. The solution method for the contact problem for homogeneous finite bodies [1], and its modification for composite bodies [2] was presented previously. Additionally a fracture submodel and contact fracture in homogeneous finite bodies with application to the indentation fracture of brittle materials has also been reported [3]. Therefore the focus of this paper deals with model adjustments for application to interfacial fracture in composite bodies due to a contact load with sliding friction. One of the most important applications of the model under consideration is the evaluation of adhesion of hard coatings in the framework of the Micro-Scratch Test [4] (see Fig. 1).
A planar surface of the coated articles is metallographically polished, normal to the interface, exposing a portion of the interface between coating and substrate. A static load, via a diamond indenter of predetermined mass and geometry, is applied to a starting point on the substrate in a direction generally normal to the planar surface. The indenter is moved across the surface from the starting point in a straight line towards and across the substrate/coating interface. Loading and scratching steps are repeated with increasing loads. The value of the applied load at which debonding (or interfacial crack formation) is observed, and the length of the crack, are used to characterize adherence. Simulations with the model demonstrated that for the Micro-Scratch Test (using a Knoop indenter) the dominant mechanism of the interfacial debonding was different than in indentation fracture of the homogeneous materials [5] (though the latter can be of much greater importance for the initiation of the same cracks). Numerous experiments with the micro-scratch technique 0921-5093/96/$15.00 © 1996 -
Elsevier Science SA All rights reserved SSDl 0921-5093(95)10141-1
v.L.
Rabinovich, VK. Sarin
Materials Science and Engineering A209 (1996) 82-90
show that this test, with the aid of a numerical model, can be used to quantify coating/substrate interfacial strength and to characterize the factors influencing it.
2. Problem statement
Two three dimensional linearly elastic bodies, one of which is a homogeneous isotropic body ~y = QI, and another one Q2 is composed of two homogeneous isotropic bodies Qi and Qk bonded through a interfacial boundary Pk, are interacting, as shown in Fig, 2. We assume that the surface displacements to be fixed on f~) and that on f} surface traction to be specified (x=I,2):
UX(O=O, t'(O =
x = i, k;
justifiable for the case of the scratching by diamond indenter. The solid advantage is that it circumvents the necessity to work with nonconvex and nondifferentiable (after Gateaux) functional in the full version of the Coulomb friction formulation, or convex and nondifferentiable functional in the version with sticking allowed, which creates significant technical difficulties [6]. In the equilibrium configuration (Navier-Cauchy equations): in
~ E f},
(I)
where t~ is an external traction distribution. On the f:~ -uncracked part of pk, the perfect joint conditions are imposed:
~Ef:~
(2)
The interfacial crack region tion free:
f(k is assumed to be trac-
ri(0=0,
tk(O=O,
~Eflk
(3)
The contact area is assumed to be small so that only displacements perpendicular to the common tangential plane at a representative contact point have to be taken into account. Furthermore, the two surfaces q, (x = 1, 2), can be identified with their projection r c to the tangential plane (contact plane). These assumptions are crucial for simplification and linearization procedures in derivation of the general contact conditions [6], which can be rewritten for our particular case [3] as
p(O=o
for
p(O~O
for w,,=g(o,
w,,>g(O ~Efc
(4)
includes the real contact area [3], p(x) is the contact pressure, w" is relative normal displacements of foundation Q2 with respect to the indenter QI, and g(x) = () Mx), () is the depth of indentation and Mx) is the gap between foundation and indenter before the deformation. For the tangential component, we have a sliding contact with Coulomb's friction law, or more precisely, restricted Coulomb friction with no sticking [6]: =
(6)
3.1. Boundary integral equations
According to Ref. [7]. displacement u and traction t fields on the regular surface of each body under consideration admit the representation:
±u(l7) +
(5)
k lric It,,1
We believe that this additional assumption
IS
quite
i
T(x,
to .u(x) dr, =
i
U(x, 17) . t(x) dr, (7)
where
16
I (I
1 )-[(3-4v)6",I/+r",r,,]
U/""
=
T"III
= - 8n\~ -21') ~ [~ren (6 I'
nil
(8)
- v r
r~
11 /fl
.1'",1',,)
+-31 2v
- n"r", + n",rl/ ]
(9)
I ~ nI, n ~ 3. In addition, ;. and fl are Lame coefficients, and)! = ;./(2(1c - fl)) is the Poisson ratio, of a corresponding body, 611 ,,/ is Kroneker's symbol. In Eqs. (8) and (9) r is the distance between the field point x and the load point t]: r =
f c is the candidate contact zone, which a pnon
Itrl
(J = i, k, I
the problem becomes one of finding a displacement field, ui ; stress field, au; the contact pressure, p; and stress-intensity factors K or energy release rate G.
ui(O = uk(o, ti(O= -tk(O,
Q/I,
3. Modified model
~Ef1
t~(O,
83
rL",~3
I
(x", - 17",)2
]L2
(10)
and all differentiation is done with respect to the field point x. The interior stresses could be obtained by differentiation of the Somigliana's Identity [7]:
a(x)+
r uT(O'S(x,Odf= Jrr tr(O'D(x,Odf(ll)
JI
with
Ski.j= 4n( I
i1
- v.) r 1
{ar . . . 3- [(I-2v)o,/rk +v(Oik r /+ 6 r i) an
.
' ./k .
- 5r,r,r k ]+ 3v(n ir/. k +n/I.k) +(1- 2v) x [3n kr ir j
+ nAk + ni()/k]- (I
- 4v )nk6 u
~
(12)
84
V.L. Rabinovich, v.K. Sarin / Materials Science and Engineering A209 (/996) 82-90
I
D kIJ·· = 8n(1 _ v)r2 {(I - 2v)[J Ik,rJ + Jkr .- J·r k] J.I II •
+ 3r,ir/,d
Indenter Candidate contact zone
(13)
3.2. Multiregion boundary element technique Discretization of the boundaries r /l , fJ = i, k, I into Mil planar triangular elements, and approximation of the surface displacements ull and tractions til
Coating
t-------,
!
Substrate r.:
Mil
ufl(x) =
I
(14)
(15)
'" Crack
i~l
Mil
tll(x) =
I
Fig. 1. A schematic diagram of the Micro-Scratch Test.
j=1
where
Mil
1ull (1J.i) + I
In=
Tfnj' ull(xm ) = I
I
In=
Ufnj ' tll(xm )
(16)
1
which can be rewritten, using given boundary conditions, into canonical form (a = 1,2): A~'
X"= B"
3.4. Approximation of the variational inequality in the contact problem Approximating the unknown contact pressure p E H1/2(r c) similar to the surface displacements and tractions (Eqs. (14) and (15»:
(17)
The latter system was solved by means of Crout LV decomposition with subsequent backsubstitution [8]. Multi-Region Boundary Element submodel was used whenever a need for solution for the mixed fundamental boundary value problem (BVP) occurs.
Me
p(X)
=
I
where Pm ~ 0 E R,
3.3. Reciprocal variational formulation for the contact pressure Rather complex non-linear contact boundary conditions (Eq. (4)) dictated for an adequate method for its treatment. For this purpose a variational approach has been adopted. Introducing a Green's operator for corresponding BVP, and decomposing all fields into two components, which accounts for boundary tractions and contact pressure, one can derive variational inequality for the contact pressure p [3]: pEKR:
VqEK R
(20)
m=l
min
M .(p) =
pEKMC
(
Me
Me
1m~lj~1 I I AmjPmPi - i~l I BmPm
where (22)
After solving for the contact pressure distribution, multi-region boundary element technique was used to
(18)
where K R = {q E (H 1/ 2 (rc )'Iq ~ a}, C -)ry is the duality pairing on (Hl/2(r~)' x Hl/2(r~), yt -normal trace operator, and g = g - y~ (w) modified gap function. Equivalent constrained ~inimization problem on K R can be derived from: min (p) = 1
(21 )
.
Contact Zone
n
(19)
In general the Green's operators are not available in explicit closed form, except for a few simple cases such as semi-infinite domains. Thus the relative normal trace of the surface displacements wand the modified gap function g must generally be approximated by the boundary element technique.
Crack Fig. 2. Problem formulation.
v.L. RahillOl'ich, V.X Sarin
85
Materials Science and Enl{ineerinl{ A209 (1996) 82 -90
recover remain surface and internal fields. The numerical procedure for solution of the contact problem is very computationally intensive. To overcome this difficulty, a parallel algorithm was developed and implemented on the massively parallel computer CM-2/S using eM-FORTRAN (FORTRAN 90 dialect). For the case of 476 boundary elements the parallel version was about 10 times faster than the serial on the 64 nodes CM-5 and Sun 4/480 platforms respectively [9]. However parallel code had severe memory limitations about 500 boundary elements, which was crucial for the complicated cases of indentation scratching fracture modelling. Much better results have been achieved by using parallelized BLAS (Basic Linear Algebra Subprograms) routines in our code on SGI Challenge Series multiprocessor system (8 MIPS processors). In the case with 964 boundary elements the use of BLAS routines speeded up the serial code on the same machine by factor 16. 3.5. Sliding contact with Coulomb's Fiction
y
o
x
Fig, .\. Fragmcnt of the crack triangulation. Shadcd area corresponds to thc crack front perturbation,
geometries. For sake of efficiency, ultimate linear system (Eq. (17)) was solved by means of bordered system solver [2, 10]. Crack growth phenomena is not addressed in this paper. In all examples below the pre-existing cracks had fixed prescribed geometry, and the energy release rates due to the indentation stress/strain fields were calculated as a fracture criterion. Von Mises yield criterion was used to assess a location and a shape of the plastic zone:
A more accurate treatment was given to the modelling of friction. A coupling between the tangential component of the contact interaction (friction traction) and the normal component (contact pressure) was taken into consideration in the following iterative manner. First distribution of the contact pressure was obtained for non-frictional case, as discussed above. Next for current pressure friction tractions were found based on Eq. (5). This distribution of the tangential tractions effects the normal displacements in the contact zone, which in turn alters the pressure distribution. To find a new pressure, we need to estimate a correction H'/Iil to the modified gap function, which describes a local geometry for the contacting bodies. This correction was found as a surface displacement for the tangential field (Eq. (5)), and new modified gap function
where (Ji are the principal stresses, Y is a yield stress for simple tension. Inside the plastic zone, the stress distribution was evaluated through the revised expanding cavity model for the elastic/plastic indentation of the brittle materials [II]:
g = g - y t (,t' + w/lie )
(j=---,
(23)
Then we iterated for modified g, repeating the same scheme, until contact pressure stabilized (it usually happened after 3-5 iterations). 3.6. Fracture mechanics and effect oj'plasticit)'
For evaluation of the interfacial fracture, a general criterion of critical energy release rate was adopted: (5U
G=(5A
(24)
where (5 U is a change in elastic strain energy owing to the perturbation (5A of the crack at some point at its front (see Fig. 3). (5 U was obtained as a result of two model simulations for original (unperturbed) and perturbed crack
I'p
(26)
(1 - 1')
where p, (j are the vertical and lateral pressure respectively, I' is Poisson ratio for the substrate and coating. Stress distribution assessment for the interfacial plane was performed on the same principle as in Eq. (26). Namely, elementary cube which in this case consist of two materials subjected to different vertical pressure pi, p2, experience horizontal tractions (j\, (j,. The values of (j\ and (j, for which the strain energy density is a minimum are the one that we are looking for [11]. It turned out that these values (j" = (j, are quite close to the average of (j obtained for substrate and coating from Eq. (26). Plastic contact pressure p was assessed based on the Chiang-Marshall-Evans model [12], using the values of the ratio E / H. This part of the model was implemented with a Mathematica script. I<
86
VL Rabinovich, V.K. Sarin / Materials Science and Engineering A209 (/996) 82-90
3.7. General modelling schemes
4.2. Micro-Scratch test modelling
The contact and crack problems were solved in a quasi-static state. The position and the shape of the interfacial crack were considered based on experimental data, and were fixed in each instance of modelling. Two different approaches were used for modelling of the contact stress field acting upon the crack. In the first case these stresses were obtained for the uncracked specimen in the location of the assumed crack as a combination of the plastic near-field and elastic far-field owing to the sliding contact of this uncracked specimen and an indenter. This method of construction of the stress field for fracture problems is traditional for indentation fracture studies [13,14]. We used this approach successfully for modelling indentation fracture in ceramics [2]. In the case of interfacial fracture this method did not produce satisfactory results, even after changing to the more detailed calculations of the near and far fields. In the second approach the stresses on crack edges were obtained from the solution of the contact problem using a preexisting crack of given geometry with sliding friction for the indenter and a specimen. We also used literature results [12] to evaluate the plastic contact pressure based on the solution for the elastic contact problem. Based on stress distribution on the crack edge releated to the position of the indenter, the energy release rates G were calculated. Finally, the energy release rate for two corner surface points of the crack front which positions determine the surface crack length, which in turn is available from experiment, was maximized with respect to the location of the indenter. These maximized values of G represent the quantitative measure of the interfacial fracture toughness.
All calculations were performed for a titanium carbide coated cemented carbide TiC/WC-Co specimen. The elastic constants were chosen to be Ewc = 695 GPa, Vwc = 0.24, and ETic = 460 GPa, VTiC = 0.188. [16]; Ywc = 18.4 GPa, YTiC = 24.2 GPa. Estimation of the yield stress Y was obtailled by using the relationship Y = H /3 [17]. Though this approximation is somewhat coarse for this type of materials, the influence of Y on the shape of the plastic zone is minimal due to the steep changes in stress field under the indenter. The thickness of the coating was measured to be 6 jim, and the geometry of the modelled system is shown in Figs. I and 4. Our first modelling strategy was based on a well-established set of hypotheses for indentation fracture of brittle materials [13], namely: (a) cracking occurs as a Mode I fracture; (b) there are two major components of the crack edge loading: an elastic field outside the plastic zone, and a residual field inside the plastic zone. Since fracture (debonding) was assumed to take place at the interfacial plane, we focused on the normal component of the tractions at the interface. Unlike the majority of the models which use idealized elastic fields (such as the Boussinesq concentrated point load field), our approach uses detailed stress distributions, based on accurate formulation and rigorous solution of the 3-D contact problem. Fig. 5 demonstrates the contact stress distribution under the Knoop indenter, positioned at the center of the 6 jim thick coating. The effects of material "softening" near the specimen edge free boundary (small values of the y coordinate) and material hardening in proximity to the interface with the stiff substrate (large ys) can be observed through the relative decrease and increase in contact pressure in corresponding areas. The same type of behavior for the pressure was reported for different shapes of indenter in [3,4,15].
4. New results and discussion
4.1. Model verification
Numerous verification examples for the various submodels have been presented in previous publications. Articles [3,4,15] contain comparisons of the model simulations to the few existing analytical results for 3-D contact problems for semi-infinite bodies. In Ref. [2], test examples for internal stress evaluation, and solutions for surface and buried crack problems were demonstrated. In the same publication application of the model to the fracture of Si 3 N 4 for Knoop scheme identation was considered. Better agreement with experimental data was shown in comparison with some other numerical and analytical models for the case of low-load indentation fracture.
Fig. 4. Wireframe model for Micro-Scratch Test simulation.
v.L. Rabinovich, v.K. Sarin
Materials Science and Enr;ineerinr; A209 (1996) 82-90
87
b)
d)
c)
Fig. 5. Contact pressure distribution.
The contact pressure variation for the changing indenter position from inside the substrate to the final one in the center of the coating, is presented in Fig. 6. In accordance with intuition, the contact pressure gradually decreases, while the size of the contact zone increase correspondently (substrate has higher Young's modulus than coating, and all simulations were performed for the fixed indentation load). As soon as the contact pressure distribution was found, the location and shape of the plastic zone was evaluated, based on the internal stress field and von Mises criterion (Eq. (25)). To maximize the impact of the contact field at the interface, all the simulations in this part were done for the indenter tip positioned at the substrate/coating interface. The cross-section of the plastic zone at the interfacial plane is shown in Fig. 7. The inner shape is related to the coating yield stress YTiC , and the external shape corresponds to the substrate yield stress Ywe. Making a reasonable assumption that owing to diffusion, the interfacial boundary is a thin zone of continuously changing properties from substrate to coating, rather than a geometrical plane with a sharp jump in properties, we took an average value for the interfacial yield stress, and obtained a cross-sectional shape of the plastic zone as intermediate silhouette in Fig. 6(a). In Fig. 6(b) this shape is shown
again with its elJiptical approximation, which was used in the model. Inside the plastic zone the crack opening tractions were calculated as in Section 3.6. Regarding the elastic component (interfacial tractions outside of the plastic zone), we used two different approaches. In the first one, the far-field approximation was used. However, instead of idealized Boussinesq field [13], we used a stress distribution under the Knoop indenter, obtained from our contact submodel simulation. This strategy worked well for the indentation fracture modelling in homogeneous materials [2]. In the case of interfacial fracture of the inhomogeneous samples this method yielded very low values for energy release rate, less than 0.5 J m -1. In order to provide a fracture mechanics submodel with more accurate distribution of the interfacial tractions inside the crack, we attempted to switch from the far-field approximation to the elastic stress distribution according to the approach [II]. This method effectively promotes an idea, that the stresses in the elastic matrix outside of the yielded zone can be found by distributing the stress field (p, q, q) from Eq. (26) over the surface of the yielded zone. Thus, to find the elastic component of the interfacial tractions we removed the material comprising the plastic zone from the specimen, distributing the tractions, corresponding to the field (p, q, q), over the cavity
V.L. Rabinovich, v.K. Sarin / Materials Science and Engineering A209 (/996) 82-90
88
-8 10
l...•....•...,....•.....,
".. _
~.~
..
--')
Fig. 6. Contact pressure distributions for moving indenter positioned (a) inside the substrate; (b) indenter tip still in the substrate, but close to the interfacial boundary; (c) indenter tip just crossed the interface; (d) in the center of coating.
surface, and then solved emerging BVP. Fig. 4 demonstrates the wireframe model for the evaluation of the elastic component of the crack opening tractions. Owing to symmetry, only half of the specimen was computed. The plastic zone exclusion can be seen. The resulting combination of the residual and elastic components of the crack opening tractions for 0.245 N indentation load is shown in Fig. 8.
Fig. 7. Plastic zone shapes.
Fig. 8. Crack opening tractions.
As can be seen, inside and in close proximity to the plastic zone, tractions tend to be crack opening. This corresponds to the negative values in our setup. Close to the periphery, they become crack closing (positive). For the final phase of the simulations, the interfacial tractions were superimposed on the interfacial crack geometry obtained from experiment. Scratch Test experiments were performed on the uncracked specimens (Micro-Scratch Test experiments were performed by Johanna Lindkvist at Boston University as a part of her thesis at the Stockholm Royal Institute of Technology) (see for details of experiment Ref. [I]). Currently only the surface lengths of the interfacial cracks are available. Therefore semi-circular and semi-elliptical cracks have been considered to date, and it was found that the influence of crack shape on the energy release rate at the surface point of the crack front is marginal. With a revised approach to the evaluation of the elastic component of the crack opening tractions, we managed to obtain surface energy release rates up to 1.0 J m - 2 for the 0.245 N indentation load and an average semi-elliptic crack of 13.3 lim
v.L. Rahinorieli, v.K. Sarin
lvlilleria/.\ Science and Engineering A209 (/996) 82 90
89
Table I Experimental data and results of numerical simulation or the Micro-Scratch Test ----
-----------------
Load
Crack length, experiment
(g) 10
•....
_.
(jim)
Conlidencc interval 1'01' experimental data (1 1m )
Energy release rate range (Jm C)
13, 12, 14, 13,
:121726,IU274:
26.4 54.3
16.4-43.4
__
.~
.. __. _ - -
12, 14, 14, 12 25
34, 26, 20, 28 36. 28
:22.0561, "<'2773 :
22.9 56.1
11.2-4X.4
50
50. 55, 65, 62, 60
: 50.152. 66.648:
22.3 48.2
6.1-46.3
surface length and a depth of 10 jI m. This level of energy release rate was considered to be very low. This can be attributed to the subsidiary feature of the Mode I fracture driven by the combination of the residual and elastic stresses under the assumptions listed at the beginning of this section. However, some inaccuracies in the input data for the model, and!or inadequacy of the simplified plasticity modelling for this rather complicated object of simulation, could also be factors. To investigate this hypothesis, we pursued a different strategy: namely, to solve the contact and crack problems simultaneously, assuming more general Mixed Mode crack propagation as opposed to Mode I fracture. In other words, the contact stress distribution was determined for the cracked specimen, and then the energy release rate for this stress/strain state was assessed for the given crack geometry. As a final step, the value of the energy release rate at the surface point of the crack front was maximized with respect to the position of the indenter on the sample. This approach has an additional advantage that the crack surface stays free of tractions (which means that crack edges are not overlapping), thus there is no need to address the very complicated issue of oscillations in the stress field in proximity to the interfacial crack front. All simulations in this approach were performed for one coated specimen (see description above) for three levels of indentation load: 0.098, 0.245 and 0.49 N (10, 25 and 50 g). The results are summarized in Table I. In column 2, experimental values of the surface length for the interfacial cracks are given. Column 3 contains the confidence intervals for the experimental data from the previous column at a 95'/';, confidence level. The next two columns present the results of numerical simulations for the crack sizes from column 3: first value corresponds to the lower limit of the confidence interval, and the second number represents the maximized G for the upper limit. The elastic contact pressure, and plastic pressure distributions were used to obtain the results listed in columns 4 and 5 respectively. The latter was obtained under following assumptions: the boundary of the plastic zone coincides with the elastic
contact zone boundary, and the pressure distribution in each subregion (substrate and coating) is uniform and proportional to the pressure from model [12] for the corresponding materia\. As the indentation load increases, the position of the indenter at which energy release rate G is maximized moves further away from the interface toward the free edge of the coating. For a load of 109 the distance from this point to the interface was 2.7 lim, for 25 g it was 4.2 jim, and for 50 g, 4.5 jim. It can be seen that maximum G's from elastic and plastic distributions are quite close. Moreover, this mechanism interfacial fracture, in the framework of the Micro-Scratch Test, looks more appropriate, than the one that was previously used to study the indentation fracture in homogeneous materials. However the Mode I fracture driven by the residual field can be dominant mechanism at the initial stages of the fracture. We think that implementation of the full-fledged elastic! plastic 3-D model can improve to some extent the outcome by taking into consideration plasticity at the crack tip and other effects associated with plastic flow, but can give rise to severe technical problems in implementation.
5. Conclusions The model for interfacial fracture in composite materials was applied to the complicated simulation of the Micro-Scratch Test for adhesion evaluation. Based on both experimental data and the results of the numerical simulation, values of energy release rate for interfacial cracks were obtained. Different mechanisms for interfacial fracture were considered. Possible applications for the model include: • Estimation of the bond strength or adherence of multi-layered material systems: • Design and development of new/improved laminated materials and coatings: • Estimation of the interfacial strength between fihers and the matrix: • Development of composite materials.
v.L.
90
Rabinovich,
v.K.
Sarin / Materials Science alld
References [J] V.L. Rabinovich, S.R. Sipcic and V.K. Sarin, Three dimensional unilateral contact problem for finite bodies, J. Appl. Mech., 6! (1994) 54 59. [2] S.R. Sipcic, V.L. Rabinovich and R.A. Rugg, Boundary variational formulation of a frictionless contact problem for composite finite bodies, AIAA J., 32 (1994) 365-370. [3] V.L. Rabinovich and V.K. Sarin, Modelling of the indentation fracture in brittle materials. 1. Mater. Sci. En~., A206 (1996) 208 214. [4] V.K. Sarin, Micro-scratch test for adhesion evaluation of thin films, J. Adhesion Sci. Techno!., 7 (1993) 1265-1278. [5] D.B. Marshall, Controlled flaws in ceramics: A comparison of knoop and vickers indentation, 1. Am. Ceram. Soc., 66 (1983) 127-131. [6] N. Kikuchi and J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities alld Finite Element Methods,
Studies in Applied Mathematics, SIAM, Philadelphia, 1988. [7] T.A. Cruse, Numerical solutions in three dimensional elastostatics. IIlI. J. Solids and Structures, 5 (1969) 1259-1274. [8] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes ill C. Cambridge University Press, Cambridge,
En~ineering
A209 (/996) 82-90
1988. [9] V.L. Rabinovich and S.R. Sipcic, Unilateral contact problem for finite bodies ~. parallel implementation, Computational Mechanics, 13 (1994) 414-426. [10] W. Govaerts, Stable solvers and block elimination for bordered systems, SIAM 1. matrix anal. appl., !2 (3) (1991) 469-483. [II] E.H. Yolfe, Elastic stress fields caused by indenting brittle materials, Philos. Mag., 46 (1982) 617 628. [12] S.S. Chiang, D.B. Marshall and A.G. Evans, The response of solids to elastic/plastic indentation. i. stresses and residual stresses, 1. Appl. Phys., 53 (1982) 298-311. [13] B.R. Lawn, A.G. Evans and D.B. Marshall, Elastic/plastic indentation cracking in glasses and ceramics, 1. Am. Caam. Soc., 63 (1980) 574-581. [14] L.M. Keer, T.N. Farris and J.-c. Lee, Knoop and vickers indentation in ceramics analyzed as a three-dimensional fracture, 1. Am. Ceram. Soc., 69 (1986) 392 - 396. [15] V.L. Rabinovich and V.K. Sarin, Three dimensional modelling of elastic stresses due to indentation of coated samples, J. Hard Mater., 5 (1994) 229-239. [16] J.F. Shackelford, (ed.), The CRC Materials Science alld EIl~i neerin~ Handbook, CRC Press, Boca Raton, FL, 1992. [17] R. Hill, The Mathematical Theory of Plasticity, Clarendon, Oxford, 1950.