Materials Science and Engineering A 425 (2006) 260–267
A Monte Carlo finite element simulation of damage and failure in SiC/Ti–Al composites P.F. Liu a,b , J.Y. Zheng b,∗ a
b
Department of Engineering Mechanics, Zhejiang University, Hangzhou 310027, China Institute of Process Equipment and Control Engineering, Zhejiang University, Hangzhou 310027, China Received 17 January 2006; received in revised form 23 March 2006; accepted 23 March 2006
Abstract A Monte Carlo 2D finite element model of the damage and failure process of a composite is established to study the stress concentration around a fiber break and the composite tensile properties. The statistical fiber strength is considered to describe the random mechanical behavior of composites. The typical damage types include fiber element failure, matrix yielding and failure, and interface debonding. A stress-based interface failure criterion combining the interfacial normal and shear effects is adopted and subsequent debonded interfaces may slip against each other obeying the Coulomb friction law. The real-time interface failure is simulated by releasing the coupled duplicate nodes at the bonded interface and the failure of fiber element is simulated by using the birth-to-death element technique. Effects of the interfacial shear strength and friction coefficient are addressed and the results are compared with those obtained by existing models and experiments. © 2006 Elsevier B.V. All rights reserved. Keywords: Mechanical properties; Finite element analysis; Monte Carlo simulation; Fiber composites
1. Introduction Continuous SiC fiber-reinforced Ti–Al composites exhibit favourable performance for their high stiffness, strength and high temperature resistance. The fracture process of composites includes a variety of complicated mechanisms such as fiber fracture, matrix yielding and failure, and interface debonding. Among the failure mechanisms above, matrix yielding and interface debonding are essentially deterministic. However, the failure of brittle fibers is essentially stochastic due to the random distributions of microscopic flaws in the fibers [1,2]. In terms of the statistical distributions of fiber strength, random fiber breaks occur during loading, leading to the stress redistributions between the intact fibers near the broken fibers. Zweben and Rosen [3] first used a statistical theory to study the stress concentration factor (SCF) around the broken fibers and the composite strength. Jansson et al. [4] first examined in detail the strength of Ti metal-based composites by relating the stress transfer around the failed fibers to the mechanical properties of the fiber and matrix as well as interface. Because fiber failure
∗
Corresponding author. Tel.: +86 571 87952110; fax: +86 571 87952110. E-mail address:
[email protected] (J.Y. Zheng).
0921-5093/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2006.03.053
due to the stress concentration may dominate the overall fracture process and the failure pattern of composite, more efforts were then focused on the composite tensile properties by exploring the loading sharing mechanism between the broken fibers and intact fibers [5–10]. Because it can effectively account for the strength distributions of fibers, the Monte Carlo method has been proved to be an effective numerical approach to simulate the failure and damage process of composites. Zhou and Wagner [11] derived the stress concentration factors for a partly debonded interface, but the random behaviors are not included. Ismar and Streicher [12] investigated the tensile stress–strain curves in detail by applying the Monte Carlo finite element method, yet neglecting the fiber break and stress concentration. Cheng et al. [13] considered the interface as an elastic-perfectly plastic body with weak bond strength to study the whole failure process of the composite. However, the friction effect is neglected and the matrix is considered only linear-elastic. Xia et al. [14–16] employed the Monte Carlo finite element method, the Green function method and the shear-lag model respectively to explore the stress concentration, and the multiscale damage and failure of composite. However, the perfectly debonded interface with zero interface bonding strength fails to describe the real-time interface debonding process.
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
This paper establishes a Monte Carlo 2D finite element model following the basic thoughts of Xia. The real-time stochastic fiber element breaks, the matrix yielding and failure, and the interface failure are simulated with increasing applied loads. The effects of the interfacial shear strength and interfacial roughness on the composite damage and failure process, the stress concentration on the surrounding intact fibers arising from the broken fiber and the tensile stress–strain relationships of composite are addressed. 2. Finite element analysis 2.1. Analytical model An array of unidirectional fibers embedded in a homogeneous matrix is adopted for the finite element analysis. The fiber radius is rf = 71 m [16], the space between two neighbouring fiber surfaces is d = 213 m and the fiber volume fraction is Vf = 40%.The finite element analysis is performed with the ANSYS (version 9.0) code. Taking into account the extreme complexity for simulating the real-time interface debonding process by using a 3D finite element model, a plane stress finite element model is established. With a single central broken fiber, the model consists of five fibers embedded in the matrix (the central fiber has two neighbours on either side). Due to the symmetry, only two intact fibers adjacent to the broken fiber are considered. The analytical model with boundary conditions and loads is shown in Fig. 1. A Cartesian coordinate system is established and the origin is located on the central axis of the broken fiber, and the x-axis represents the axial fiber direction and the y-axis denotes the direction perpendicular to the fiber axis. The symmetric constraints are exerted on the symmetric sides and the uniform displacement u is applied on the right side of model. The four-node and four-sided plane element is employed. The axial fiber length L = 8.875 mm is chosen such that the right part of model is not disturbed by the stress redistributions due to the local fiber damage at the left part of model. The length x = L/100 mm of each element in the fiber direction is selected, which is small enough to reflect the local stress concentration near the broken fiber. Fiber break is simulated by removing a single row of fiber elements where the maximum stress concentration may occur theoretically instead of using an one-quarter crack tip element. The model consists of 3497 elements and 3128 nodes, including 998 fiber elements.
Fig. 1. The analytical model with boundary conditions and loads.
261
Table 1 Mechanical properties used for numerical calculations [16]
Fiber Matrix
Young’s modulus (GPa)
Poisson’s ratio
400 115
0.17 0.30
Yielding strength (MPa)
Tensile strength (MPa)
950
1275 [17]
2.2. Material properties The SiC fiber-reinforced titanium alloy composite SiC/ IM1834 is adopted for the finite element analysis. The fiber is considered to be isotropic and linear elastic. The matrix is considered to be isotropic and elastic–plastic, and obeys the von Mises yielding criterion. According to the Ramberg–Osgood material model, the matrix stress–strain relationship is expressed as [16]: σm σm 1/b εm = + (1) Em B where b = 0.0384 is the hardening exponent, B = 1229 MPa is the hardening coefficient and Em is the matrix Young’s modulus. The material properties of the fiber and matrix are listed in Table 1. The temperature dependent material behaviors are not included in this study. According to the multi-linear hardening theory, a series of matrix stress–strain points satisfying Eq. (1) after matrix yielding are imported into the ANSYS software and the matrix material curve is effectively fitted. 2.3. Interface failure criteria Research on the stress-based interface failure criteria has been widely conducted to explore the composite mechanical properties. Cheng et al. [13] considered the interface as an interphase with independent material properties and thickness, distinctly different from the fiber and matrix. When the interface is assumed to be a spring stiffness layer with negligible thickness, Chandra and Ananth [18] presented a quadratic stress failure criterion using an isotropic debonding surface. The criterion is based on the combined action of the practical interfacial normal stress σ r and shear stress τ, and is given by 2 σr τ 2 + ≥1 (2) σrf τf where σrf and τ f are the interfacial normal and shear strength, respectively. A relatively high value σrf = 30 MPa is adopted in our simulation for exploring the effect of τ f . Fig. 2(a) shows the failure surface used for the debonding criterion. It is noted that as σrf increases to infinity in the limit case, the interface failure criterion in Eq. (2) is reduced to a constant shear strength-based criterion. Initially, duplicate nodes are created at the interface on fiber and matrix sides. All degrees of freedom of those duplicate nodes are coupled to ensure complete bonding of the fibers and matrix. When Eq. (2) is satisfied, the coupled nodes are released, leading to the debonding of inter-
262
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
Fig. 2. Schematic representations of (a) the debonding surface, and (b) the friction stress at the debonded interface.
face. After interface debonding, the debonded nodes may contact and slip against each other obeying the Coulomb friction law, yielding τ ≤ τ0 + µσr
(3)
where τ 0 is the interfacial residual shear stress before loading and can be determined once the thermal residual stress field is obtained. µ is the friction coefficient. Fig. 2(b) shows the increasing tendency of τ, which increases proportional to the normal stress up to a limit plateau. The Coulomb friction law is realized by using the augmented Lagrange contact algorithm. The fiber and matrix are defined as the “Contact” body and “Target” body, respectively. The contact normal modulus is chosen large enough to prevent the fibers and matrix from penetrating each other. The friction coefficient µ is taken as a material parameter of the interface elements. Fig. 3 shows the coupled nodes at the interface before loads. The kinks represent the node coupling and the bold lines at the interface denote the interface contact elements. 2.4. Monte Carlo simulation For the fiber elements, the maximum tensile stress criterion is introduced as the failure criterion. According to the two-parameter Weibull function and by random sampling, the strength distribution of each fiber element can be written as [16,19]: x σf m F (σf ) = 1 − exp − . (4) L0 σ 0 where F is the failure probability of fiber elements at a stress level equal to or less than the axial fiber stress σ f and is a random number within [0,1], m the shape parameter that characterizes the flaw size distributions in the fibers, and σ 0 is the scale parameter
at a gauge length L0 that describes the characteristic strength of fibers. Each fiber element is assigned a random failure strength consistent with Eq. (4), that is, a random number η within the interval [0,1] is selected and then the tensile strength σ fu of each segment is expressed as 1/m L0 σfu = σ0 − ln(1 − η) . (5) x where m = 17, L0 = 25 mm and σ 0 = 4.58 GPa for SiC fibers. When the applied tensile stress at all Gauss integration points in a fiber element or matrix element reaches the corresponding tensile strength, the element is considered to fail. Then, we multiply an extremely small factor 1 × 10−6 on the stiffness matrix of the failure element to deactivate the “birth” element to “death”. In this case, a deactivated element still remains in the finite element model, but contributes a near-zero stiffness. Also, other physical characteristics of the elements can be set to zero for deactivated elements. To model the real-time interface debonding, the failure of fiber and matrix elements, specific subroutines by using the ANSYS–APDL (ANSYS Parametric Design Language) programming language are linked to the main program. The detailed flow chart of the Monte Carlo simulation is shown in Fig. 4. The finite element analysis is applied in two steps. In the first step, by using the Newton–Raphson procedure, an equally increasing displacement is applied to determine the stress–strain curves before the critical failure point, where a restart analysis is performed after each loading step. In the second step, the critical failure point and the decreasing part in the tensile curve are obtained by activating the arc-length method. A displacementcontrolled convergence criterion is employed to control the analysis and about 20 loading steps are applied. 2.5. Tensile stress–strain relationships of composite By using the Weibull distribution theory, Phoenix and Raj [20] proposed a statistical model of the fiber bundle stress σf∗ with broken fibers and the composite tensile process including the matrix cracking and frictional sliding at the debonded interface. Also, Jansson et al. [4], Curtin [5], Ibnabdeljalil and Curtin [6], Ramamurty et al. [7,10] and Xia et al. [16] presented effective methods respectively to calculate the composite stress by using two factors: one is the fiber bundle stress σf∗ and the other is the axial matrix stress σ m . The composite tensile stress σ is
Fig. 3. Schematic representations of the coupled nodes and interface contact elements before loading.
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
263
fiber. For a fiber element of length dx, the fiber failure probability in Eq. (4) can be approximated as dx σf m dx σf m ≈ (7) dF (σf ) = 1 − exp − L 0 σ0 L0 σ 0 where σ f = k(x)σ f∞ , σ f∞ is the far-field fiber stress, and k(x) is the stress concentration factor along the x-axis. According to Weibull statistics, fiber breaks if any of these differential elements fails. Then, the survival probability of the nearestneighbour fiber F¯ is obtained m δ∗ 2 σ f∞ m F¯ = 1 − F = exp − k (x) dx . (8) L 0 σ0 0 where δ* is half a length that the axial fiber stress recovery from 0 at the break surface to σ f∞ requires. The fiber stress beyond the distance δ* is not considered due to stress relaxation. Assuming there are n nearest-neighbour fibers, the failure probability of at least one of the fibers Fn is expressed as ∗ 2n σf∞ m δ m n Fn = 1−(1 − F ) = 1− exp − k (x) dx . L 0 σ0 0 (9) Then, σf∗ is calculated as σf∗ = σf × (1 − Fn )
Fig. 4. Flow chart of the finite element analysis.
2.6. Numerical results and discussion
given by σ = Vf σf∗ + (1 − Vf )σm
(10)
(6)
According to the theories above, based on the SCF near the broken fiber, Gonz´alez et al. [9] further presented another simplified method to calculate σf∗ , which requires a calculation of the failure probability dF of a fiber element adjacent to the broken
In order that the real-time interface debonding process is included, no thermal effect is considered in this paper. Fig. 5 shows the stochastic damage and failure evolution process at τ f = 10 MPa and µ = 0.25 with increasing loads. It can be seen that the initial fiber break gives rise to the stress concentration within the neighbouring fiber elements and subsequent interface
Fig. 5. Composite damage and failure process at applied strains of (a) ε = 0.225%, (b) ε = 0.45%, (c) ε = 0.676%, and (d) ε = 0.845%, respectively.
264
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
Fig. 6. Composite damage and failure process at applied strains of (a) ε = 0.225%, (b) ε = 0.45%, (c) ε = 0.676%, and (d) ε = 0.9%, respectively (× represents a broken fiber element).
debonding. With the increase of loads, more and more nodes at the interface debond and may slip against each other. Some matrix elements begin to enter the plastic deformation stage at a applied strain of about ε = 0.507% after Fig. 5(b). At a strain slightly higher than about ε = 0.845% after Fig. 5(d), a remarkable steep drop in the composite stiffness occurs, indicating the composite failure, where the far-field matrix tensile stress reaches 918 MPa at the critical strain point and thus partial matrix elements still remain the elastic deformation stage. Due to a relatively small applied strain before the composite failure, no fiber element breaks in this simulation. Fig. 6 shows the stochastic damage and failure evolution process at τ f = 100 MPa and µ = 0.5 with increasing loads. The matrix stress–strain relationships begins to enter the plastic stage at about ε = 0.563% before Fig. 6(c) and fully transforms into the plastic deformation stage at ε = 0.9% in Fig. 6(d), where the farfield axial matrix stress reaches 952 MPa. At about ε = 0.85% before Fig. 6(d), first fiber element breaks, which randomly occurs on the right side of the model. The broken fiber element overstresses the neighbouring fiber elements and leads to a relatively large area nearby with interface debonding. At a strain level slightly higher than ε = 0.9%, the composite begins to fail. No matrix element fails in this simulation. The critical failure strain ε = 0.9% based on the plane–stress assumption approaches the SiC/Ti–Al composite ductility ε = 1% including the thermal effect, which shows the thermal residual stresses improve the composite tensile performance. The failure of the composite immediately after first break is confirmed from the theoretical simulation by Ananth et al. [21] and the experimental observations by Majumdar and Newaz [22]. By comparing Figs. 5 and 6, small τ f and µ lead to weaker ability for the composite to resist the deformation and damage than large τ f and µ at the same applied strain. The introduction of the central fiber break induces significant changes in the local stresses around the break. With the
increase of applied strains, Fig. 7 shows the stress concentration factor (SCF, the local axial fiber stress normalized by the far-field fiber stress, and averaged over the fiber cross-section) of two neighbouring fibers around the broken fiber. The SCF increases with increasing applied strains and is slightly larger for a strong interface than for a weak interface, which is in good agreement with the result obtained by Fiedler et al. [23] by using the linear-elastic matrix. However, the interfacial bonding properties have smaller effects on the SCF at large applied strain, where relatively large damage in the composite occurs. By comparison, several statistical theoretical models which describe the composite failure evolution behavior have been developed such as the globle load sharing (GLS) model [5,10,24], the local load sharing (LLS) model [6,7], the HVD shear-lag model [8], the HVDP shear-lag model [25] and Xia’s
Fig. 7. Axial SCF at different τ f and µ in terms of two neighbouring fibers (‘First’ represents the nearest-neighbour fiber and ‘Second’ denotes the second nearest-neighbour fiber).
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
265
Table 2 SCF at an applied strain of 0.788%
1 2
τ f = 10 MPa
τ f = 100 MPa
µ = 0.25
µ = 0.5
1.2445 1.1017
1.2547 1.1057
GLS model [5,6]
HVDP shear-lag model [25]
1.0011 1.0011
Xia’s model [16]
1.1046 1.0100
µ = 0.25
µ = 0.5
1.038 1.017
1.057 1.105
Table 3 Composite tensile strength, as predicted and as measured τ f = 10 MPa, µ = 0.25
UTS (MPa)
1854
τ f = 100 MPa, µ = 0.5
1973
Curtin [5]
Xia’s model [16]
µ = 0.25
µ = 0.5
µ = 0.25
µ = 0.5
2225
2302
2268
2317.5
model [16]. The GLS model assumed all intact fibers share the applied loads equally and thus equally carry the loads lost from broken fibers, where sufficiently low interfacial shear stress is assumed. In the LLS model, only the immediate unbroken neighbours carry these lost loads, thereby causing more severe overloads than those that occur by the GLS model. The HVD shear-lag model considered a fiber composite consisting of infinite hexagonally aligned fibers embedded in the elastic matrix and is more valid than the GLS model because HVD accounted for the effect of elastic deformation of the fibers and matrix on the redistributions of tensile stress based on the shear-lag theory. The HVDP shear-lag model considered the geometry of finite periodic patches based on the HVD shear-lag model and applied a so-called general influence function technique to solve a 3D stress field around fiber breaks. Xia et al. [16] established a 3D finite element model of a hexagonal array of fiber composite and considered high thermal processing temperature enough to debond the interface initially compared with the models above neglecting the thermal effects. Table 2 compares the SCF values obtained by us with those by existing models at an applied strain of about 0.788%. The SCFs obtained by the GLS model are the smallest among all mod-
Experiment conducted by Leucht [26] 2200–2300
els. The GLS model decreases the stress concentration around the broken fiber and underestimates the sharing portions of the loads borne by the neighbouring intact fibers. The SCFs obtained by the HVDP model and Xia’s model are both slightly smaller than those obtained by us. This can be interpreted as: the HVDP model neglects the matrix tensile stress and friction stress at the debonded interface, and assumes the matrix to be in a state of only pure elastic shear. Therefore, the stress concentration on the surrounding intact fibers is relaxed. Xia et al. [16] considered zero interfacial shear and normal strength as well as high interfacial thermal stresses so as to debond the interface sufficiently before the applied loads. Then, the debonded interface can contact due to high residual stresess and slip against each other obeying the Coulomb friction law. Because the axial fiber stress and SCF are both functions of the applied strain, the stress recovery length δ* always changes in the loading process. Fig. 8 shows the stress profile along the central broken fiber as well as those of the first and second neighbouring fibers. The stress recovery process for the broken fiber is non-linear compared with the linear recovery process in the shear-lag model, in which the stress recovery length around the broken fiber can be estimated using a constant friction stress
Fig. 8. Stress profile along the central broken fiber as well as those of the first and second neighbouring fibers normalized by the far-field fiber stress at (a) τ f = 10 MPa, µ = 0.25 and an applied strain of 0.845%, and (b) τ f = 100 MPa, µ = 0.5 and an applied strain of 0.9%.
266
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
Fig. 9. Average axial fiber and matrix stresses with increasing applied strains.
as δ* = rf σ f∞ /2τ. This simulation predicts higher values for the stress recovery length δ* than those obtained by the shear-lag model, where the matrix is considered to be elastic and the friction stress is neglected after interface debonding. The stress recovery length obtained by us in Fig. 8(a) due to a weak interface reaches about 7.5 mm, larger than about 5.5 mm in Fig. 8(b) with a strong interface. For the nearest-neibour fiber around the broken fiber, by using the ORIGIN data software, the stress concentration factor k(x) as a function of the position x within the stress recovery region in Fig. 8(a) and (b) can be fitted as, respectively k(x) = 0.99182 + 0.10458 exp(−0.4821x), 0 < x < 7.5 mm (11) k(x) = 0.99666 + 0.20013 exp(−1.137x), 0 < x < 5.5 mm (12) Fig. 9 shows the average axial fiber and matrix stresses with increasing applied strains for the simulation above. It can be clearly seen that the curves are approximately the same at the elastic stage for different τ f and µ in terms of the fiber and matrix, respectively. The matrix stress–strain curves takes on the distinct elastic–plastic properties. Fig. 10 shows the tensile stress–strain curves of the composite. Two curves take on the same tendency at the elastic stage. The plastic stages in these curves occur before the composite failure. Table 3 shows the UTS of the composite obtained by using the theoretical and experimental methods, respectively. Increasing τ f and µ result in larger values for the UTS. The UTS values obtained by us are smaller than the theoretical and experimental values obtained by Curtin [5], Xia et al. [16] and Leucht et al. [26], respectively. The conclusion can be interpreted as: Xia and Leucht considered 750 ◦ C thermal residual temperature due to thermal processing and thus the interfacial radial and shear stresses increase, leading to a stronger ability for the composite to resist the failure. If the composite is considered as a homogeneous material with perfectly bonded interface, the Young’s modulus of composite is 229 GPa at the fiber volume fraction
Fig. 10. Composite tensile stress–strain curves with increasing applied strains.
40% and then the UTS will reach 1935, 2061 and 2290 MPa at applied strains of 0.845%, 0.9% and 1%, respectively. However, when the fiber break and interface failure are commonly considered, the practical UTS of the composite should be smaller, which indicates the validity of our results in a qualitative sense. Also, by using the classical GLS model, Curtin [5] derived an expression of the UTS σ u , given by 1/(m+1) 2 m+1 (13) σu = Vf σc + (1 − Vf )σm m+2 m+2
m 1/m m/(m+1) σ0 τL0 1/(m+1) σ 0 rf L 0 (14) , δc = σc = rf τ where σ c is the characteristic fiber strength at a gauge length δc . For the SiC/IM1834 composite, τ = 208µ MPa [16] is considered as the friction stress at the debonded interface when the thermal residual stresses are included. The GLS model predicts an accurate value for the UTS of composite at a low friction coefficient µ, but overestimates the UTS at higher µ. The conclusion can be interpreted as: the efficiency of stress transfer at the debonded interface is low at low µ when the thermal effect is considered and thus the SCF around the broken fiber is relatively low, which approaches the case that the GLS model assumes. However, at a higher µ larger than 0.5, the local loading sharing models [6,7] are strongly required to emphasize the effect of relatively high SCF on the UTS of composite. 3. Conclusions A Monte Carlo finite element model of the failure and damage evolution process of a fiber composite is presented. The random behavior in the tensile process is represented by the statistical characteristic of fiber strength distributions. A stressbased interface failure criterion combining the interfacial normal and shear effect is adopted and subsequent debonded interfaces slip against each other obeying the Coulomb friction law. The real-time interface failure is simulated by releasing the coupled
P.F. Liu, J.Y. Zheng / Materials Science and Engineering A 425 (2006) 260–267
duplicate nodes at the bonded interface and the real-time break of fiber elements is simulated by using the birth-to-death element technique. The primary work is placed on the stress concentration around the broken fiber and the tensile stress–strain curves of composite. The following conclusions are obtained: (1) At low applied strain, a strong interface results in a large SCF than a weak interface. However, the interfacial bonding properties have smaller effects on the SCF at large applied strain. (2) To some extent, the GLS model and HVDP shear-lag model underestimate the stress concentration on the neighbouring intact fibers around the broken fiber. (3) For the prediction of the UTS of composite, the GLS model and Xia’s model are applicable only to the case of a relatively low friction coefficient at the debonded interface. (4) The thermal residual stresses due to thermal processing increase the UTS of composite. Theoretically, our model is more general than Xia’s model [16] because it can account for the real-time interface debonding process in the absence of thermal stresses.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
Acknowledgement The work at Zhejiang University was supported by the National Natural Science Foundation of P.R. China (Project No. 10372091).
267
[22] [23] [24] [25] [26]
E.M. Wu, C.S. Robinson, Comp. Sci. Tech. 58 (1997) 1421–1432. B.A. Bednarcyk, S.M. Arnold, Comp. Sci. Tech. 61 (2001) 705–729. C. Zweben, B.W. Rosen, J. Mech. Phys. Solids 18 (1970) 189–206. S. Jansson, H. Deve, A.V. Evans, Metall. Trans. 22A (1991) 2975–2983. W.A. Curtin, J. Am. Ceram. Soc. 74 (1991) 2837–2845. M. Ibnabdeljalil, W.A. Curtin, Int. J. Solids Struct. 34 (1997) 2649–2668. U. Ramamurty, F.W. Zok, F.A. Leckie, H. Deve, Acta Mater. 45 (1997) 4603–4613. S. Mahesh, I.J. Beyerlein, S.L. Phoenix, Phys. D 133 (1999) 371–389. C. Gonz´alez, J. Llorca, Acta Mater. 49 (2001) 3505–3519. U. Ramamurty, Comp. Sci. Tech. 65 (2005) 1815–1825. X.F. Zhou, H.D. Wagner, Comp. Sci. Tech. 59 (1999) 1063–1071. H. Ismar, F. Streicher, Comput. Mater. Sci. 16 (1999) 17–24. T.L. Cheng, R. Qiao, Y.M. Xia, Comp. Sci. Tech. 64 (2004) 2251–2260. Z. Xia, W.A. Curtin, T. Okabe, Comp. Sci. Tech. 62 (2002) 1279–1288. Z. Xia, T. Okabe, W.A. Curtin, Comp. Sci. Tech. 62 (2002) 1141–1149. Z. Xia, W.A. Curtin, P.W.M. Peters, Acta Mater. 49 (2001) 273–287. X.N. Mao, L. Zhou, Y.G. Zhou, Rare Met. Mater., Eng. 33 (2004) 620–623 (in Chinese). N. Chandra, C.R. Ananth, Comp. Sci. Tech. 54 (1995) 87–100. T. Okabe, N. Takeda, Y. Kamoshida, M. Shimizu, W.A. Curtin, Comp. Sci. Tech. 61 (2001) 1773–1787. S.L. Phoenix, R. Raj, Acta Metall. 40 (1992) 2813–2828. C.R. Ananth, S.R. Voleti, N. Chandra, Comp. Part A 29 (1998) 1203–1211. B.S. Majumdar, G.M. Newaz, Philos. Mag. A 66 (1992) 187–212. B. Fiedler, A. Klisch, K. Schulte, Comp. Part A 29 (1998) 1013–1019. K.M. Gaffney, J. Botsis, Comp. Sci. Tech. 59 (1999) 1847–1859. J.M. Hedgepeth, P. Van Dyke, J. Comp. Mater. 1 (1967) 294–309. R. Leucht, H.J. Dudek, Mater. Sci. Eng. A 188 (1994) 201–210.