ICF10 Honolulu (USA) 2001 Vol. A
ICF100591OR 2D AND 3D SELF-AFFINE CRACK PROPAGATION ON ALUMINUM ALLOYS M. Hinojosa, E. I. Morales, X. Guerrero, J. Aldaco and U. Ortiz Facultad de Ingeniería Mecánica y Eléctrica, Universidad Autónoma de Nuevo León, A.P. 076 Suc "F", Cd. Universitaria, San Nicolás de los Garza, N.L. MEXICO. ABSTRACT The self-affine exponents associated with the crack propagation phenomenon are evaluated on samples of aluminum alloys both on 2D and 3D experimental conditions. Fracture surfaces were generated by Charpy impact tests on samples of A319-type aluminum cast alloy. Roughness exponents and correlation length on the perpendicular and parallel directions with respect to the crack propagation direction were determined, this analysis was also performed for the arrested crack propagation front. In the two-dimensional case, cracks were propagated on notched tension specimens of aluminum foil and the resulting self-affine crack paths were recorded and analyzed, the self-affine parameters were determined for both longitudinal and perpendicular direction in order to investigate the effect of the microstructural anisotropy. Self-affine analysis was carried out using the Zmax variable bandwidth method. The combined use of different techniques (SEM, AFM, Optical microscopy and stylus profilometry) enabled the analysis over up to seven decades of length scales. The results are analyzed in terms of recent crack propagation models and the self-affine parameters are found to be correlated with microstructural characteristic lengths. KEYWORDS Roughness exponents, self-affine crack paths, fracture surfaces, aluminum alloys, crack propagation. INTRODUCTION Crack propagation and the fracture of materials are catastrophic phenomena of considerable scientific, technological an economical importance [1-4]. Despite the scale of the problem and the considerable effort that has been undertaken by engineers and scientists of different disciplines, there is at present no clear understanding of the fracture process. In recent years much interest has been devoted to the self-affine character of fracture surfaces and crack propagation [5-7]. The fractal nature of fracture surfaces was first quantitatively studied in the mid-eighties [8]. At about the same time it was suggested that the fracture of heterogeneous media has some universal properties similar to critical phase transitions [9]. Later, experimental evidence led to the conjecture of the existence of a universal roughness exponent [10] for the fracture surfaces of many different materials [11], though this is still a controversial topic [7,12]. Anyway, it is now clearly established that fracture surfaces are self-affine objects that can be quantitatively described by self-affine parameters like one or more roughness exponents and one or more characteristic lengths such as cut-off lengths separating different scaling regimes, and the correlation length.
One of the main goals of materials scientists is to find clear and useful relationships between the microstructure of materials and their macroscopic properties. In our particular field of interest this traduces to finding quantitative relationships between the microstructural features and the relevant self-affine parameters associated with the fracture surface and the crack propagation process that led to its creation. From the statistical physics point of view the question is related to how disorder affects crack propagation considering that rupture is the culmination of a self-organization of cumulative damage and cracking characterized by power-law regimes which result from the fact that disorder is present at different length scales in the form of impurities, vacancies, grain boundaries, porosity, phase boundaries and so forth. The first attempts to relate fractal parameters of the fracture surfaces [8] of maraging steels with mechanical properties were very promising and inspirational though unsuccessful, it is clear that the fractal dimension of a fracture surface is not clearly correlated with the toughness of the material. Moreover, fractal dimension is not an appropriate parameter to describe self-affine surfaces [13], the roughness or Hurst exponent should be used instead. At present [6, 7], results from experiments in a wide variety of materials tested in different kinetic conditions and analyzed with different topometric techniques over up to seven length scales [14] suggest the coexistence of two self-affine regimes, at high enough propagation speeds and/or large enough length scales the so-called universal exponent ζ ≈ 0.78 is detected, whereas at slow propagation conditions and/or small enough length scales the detected exponent has a value close to 0.5. The cut-off length separating these two regimes is apparently dependent on the propagation speed [15], it also appears to be affected by local plastic deformation at the crack tip in ductile materials. Neither of the two above mentioned exponents seems to be associated in any manner whatsoever with the microstructural features of the materials. Experiments in Al-Ti alloys suggested that the cut-off length might be linked to the size of intermetallic compounds embedded in the metallic matrix [16]. Recent results [17-19] have shown that the correlation length, i.e. the upper limit of the self-affine regime(s) is directly related with the largest heterogeneities in materials such as metals [14, 17, 18], polymers [19, 20] and certain glasses [19]. With the hope to provide more experimental evidence that can help to improve our knowledge and refine the existing theoretical models of crack propagation, in this work we report the experimental analysis of the selfaffine parameters of fracture surfaces and crack paths in aluminum alloys. A cast aluminum alloy is broken in mode I and the three associated roughness exponents are recovered along with the respective correlation length is some cases. We have also tested an aluminum foil in 2D mode I condition and have analyzed the self-affine nature of the crack paths. In both cases special attention is paid to the possible relationships between the microstructural features and the self-affine parameters. EXPERIMENTAL PROCEDURE We have performed the self-affine analysis of the fracture surfaces of aluminum samples, the same analysis was done for the crack paths generated in mode I in 2D conditions using samples of aluminum foil. The quantitative analysis was carried out using height profiles which were recorded by different techniques. The resulting topometric data sets are processed using the variable bandwith method [21] in which the following quantity was calculated: Zmax(r) = < max{z(r’)}x<r’<x+r - min{z(r’)} x<r’<x+r >x ∝ r ζ Where r is the width of the window and Zmax(r) is the difference between the maximum and the minimum height z within this window, averaged over all possible origins x of the window. A log-log plot of Zmax(r) vs. r gives a straight line for a self-affine profile.
The experimental details and results for the two cases considered in our work are presented below. 3D case, A319-type Aluminum alloy The cast aluminum alloy employed for this part of our work is an A319-type alloy, which is commonly used in the automotive industry. The typical chemical compositon is as follows (wt %): Si:7.147, Cu:3.261, Fe:0.612, Zn:0.664, Mn:0.0374, Ni:0.041, Ti:0.154, Mg:0.313, Sr:0.014, Al: balance. Fig.1 shows the microstructure of this alloy as observed by optical microscopy, the dendrites of alpha aluminum-rich phase is observed along with a numbre of different phases in the interdendritic regions. There is also a grain structure which was revealed using a special preparation. Image analysis measurements showed that the largest heterogeneities were the dendrites and the grains, with characteristic lengths identified as the primary dendrite arm length (800 µm) and the grain size (950 µm). Samples of this material were broken in Charpy impact tests according to ASTM standard E-23-93. The resulting fracture surfaces were examined by Scanning Electron Microscopy (SEM), Atomic Force Microscopy (AFM) and an stylus profilometer. These three techniques were used to obtain topometric profiles both in the perpendicular and parallel direction with respect to the crack propagation direction, see Fig. 2. Profiles of the arrested crack front were also recorded using a different procedure which is described later in this section. Using these profiles we were able to determine the perpendicular out-of-plane roughness exponent ζ⊥, the parallel out-of-plane roughness exponent ζ//, and the roughness exponent of the crack front ζf, Fig. 2. Fig.1 Optical micrography showing the microstructure of the A319-type alloy. Fig. 2 Scheme illustrating the height profiles in the parallel and perpendicular directions with respect to the propagation direction, the crack front and the three roughness exponents are also included. The SEM topometric profiles in the parallel and perpendicular directions were obtained by sectioning the surfaces previously plated with nickel, then SEM images are recorded using backscattered electrons and the profile is extracted by image analysis procedures. More details of these technique can be found in references [6, 15-17]. SEM profiles of 1024 points were obtained at magnifications ranging from 100X to 4000X. The AFM profiles are directly recorded by scanning the uncoated surfaces, we have used the contact mode in air. Profiles of 512 points with scan sizes ranging from 0.5 to a maximum of 6 µm were obtained. The stylus profilometer provided us with profiles of a maximum length of around one centimeter, a typical profile consisted of around 10,000 points with resolution of 0.25 µm.
Figures 3 and 4 show the self-affine curves obtained for the perpendicular and parallel directions, respectively. The exponents ζ ⊥, ζ// have very similar values: 0.81, 0.78, repectively. 10-3 10-2 10-1 100 101 102 103 104 10-3 10-2 10-1 100 101 102 103 ζ=0.81 Zmax(r), (µm) r, (µm) 10-1 100 101 102 103 10-1 100 101 102 ζ = 0.78 Zmax(r), (µm) r(µm) Fig. 3 Self-affine curve for the profiles in the perpendicular direction, the roughness exponent ζ ⊥ has a value of 0.81. Fig. 4 Self-affine curve for the profiles in the parallel direction, the roughness exponent ζ// has a value of 0.78. Profiles of arrested crack fronts were obtained by a very different method, we have run interrumpted torsion tests over hollow cilindrical specimens then marked the crack front using a commercial penetrating die commonly used in crack inspection and failure analysis. The specimens were then broken in the torsion machine and the arrested crack front was registered by SEM using secondary electrons, profiles were extracted by image analysis. Figure 5 shows the self-affine plot for the crack front which has a roughness exponent ζf = 0.79. 101 102 103 101 102 103 ζf = 0.79 Zmax(r), (µm) r(µm) Fig. 5 Self-affine curve for the arrested crack front, the roughness exponent ζf has a value of 0.79 These self-affine curves permit only an estimation of the correlation length. However, as it can be observed, in all cases it has a value of the order of 1 millimeter, which is very close to the size of the largest dendrites and grains. 2D case: Self-affine crack propagation in aluminum foil Tension specimens of aluminum foil (alloy 1145-O) were prepared as shown in Fig. 6, then fractured in 2D mode I condition. We have then performed the self-affinity analysis of the resulting crack paths. The purpose of these experiments was to evaluate the self-affine parameters paying special attention to the possible effect
that the anisotropic grain structure might have on the self-affine character of the crack paths. As it is shown in Fig. 7, the grains are elongated in the rolling direction, it is known that this causes anisotropic behavior of mechanical properties so one can expect an analogous effect on the self-affine parameters. Rolling direction Rolling direction 20µm Fig. 6 Scheme showing the orientation of the tension specimen with respect to the rolling direction. Fig. 7 Microstructure of the aluminum foil showing grains elongated in the rolling direction. // // // // ⊥ ⊥ ⊥ 101 102 103 100 101 102 parallel perpendicular ζ = 0.67 Zmax(µm) r(µm) Fig. 8 Samples of the recorded crack paths in the rolling direction (//) and the perpendicular direction (⊥). Fig. 9 Self-affines curves for the profiles in the parallel and perpendicular direction respectively, both curves reveal that the roughness exponent ζ has a value of 0.67 The crack paths obtained as a results of the tension test were recorded at various magnifications using SEM, optical microscopy and a conventional document scanner. Samples of the recorded self-affine paths are shown in Fig. 8 where paths belonging to cracks propagating in the rolling direction are “wider” and clearly distinguishable from those propagating in the perpendicular direction. The self-affine plot shown in Fig. 9 reveals that the roughness exponent has about the same value for both directions, ζ = 0.67, this value is in good agreement with the results predicted by the random fuse model and a 2D simulation of crack propagation reported in reference [16]. It is not possible to estimate with good precision the correlation lengths but Fig. 9 suggest that this parameter is larger for the parallel direction compared to that of the perpendicular direction, one can speculate that this can be interpreted as an effect of the elongation of the grains caused by the rolling process.
CONCLUSIONS We have determined the self-affine parameters of the fracture surface of a cast aluminum alloy. The parallel and perpendicular out-of-plane roughness exponents were determined with values of 0.78 and 0.79, respectively. The roughness exponent of the arrested crack front was also determined, with a value of 0.79. It was corroborated that the correlation length is in all the cases related to the size of the largest heterogeneities present in the complex microstructure. The analysis of the crack paths in aluminum foil as developed in 2D mode I loading allowed the determination of the respective self-affine parameters. It was found that the roughness exponent has a value of 0.67 for both the parallel and transverse direction of propagation with respect to the rolling direction. The anisotropy of the microstructure has an observable effect in the correlation length whereas the roughness exponent is apparently unaffected by this condition. ACKNOWLEDGEMENTS Author express their gratitude to E. Sánchez, E. Reyes, O. Garza, I. Suárez, L. Cruz, R. Colás and E. Velasco, the financial support of the Consejo Nacional de Ciencia y Tecnología and the UANL through the PAICYT program is also gratefully acknowledged. REFERENCES 1. M. E. Eberhart, Scientific American, p. 44-51, Oct. 1999. 2. Hellemans, Science Vol. 281, p. 943-944, 14 August 1998. 3. MRS Bulletin, Vol. 25, No. 5, May 2000. 4. J.W. Hutchinson and A. G. Evans, Acta materialia Vol. 48 p. 125-135, 2000. 5. M. Marder and J. Fineberg, Physics Today Vol. 49, p. 24, 1996. 6. E. Bouchaud, J. Phys.:Condens. Matter Vol. 9, p. 4319-4344, 1997 and references therein. 7. A. Balankin, Engineering Fracture Mechanics Vol. 57, p. 135-203, 1997 and references therein. 8. B.B. Mandelbrot, Passoja, Paullay, Nature Vol 308 19 April 1984. 9. L. de Arcangelis, S. Redner, H.J. Herrmann, Journal de Physique Lettres 46, p. 585-590, 1985. 10. E. Bouchaud, G. Lappasset, J. Planes and S. Naveos, Phys. Rev. B, Vol 48, p. 2917 ,1993. 11. K. J. Maloy, A. Hansen, E. L. Hinrichsen and S. Roux, Phys. Rev. Lett., Vol. 68, p. 213-215, 1992. 12. X. Zhang, M.A. Knackstedt, D.Y.C. Chan and L. Paterson, Europhysics Letters, Vol. 34, p. 121-126, (1996). 13. J. Feder, Fractals, Plenum Press, New York, 1988. 14. M. Hinojosa, J. Aldaco, U. Ortiz and V. González, Aluminum Transactions Vol. 3, p.53-57, 2000. 15. P. Daguier, B. Nghiem, E. Bouchaud and F. Creuzet, Phys. Rev. Lett. Vol 78, p. 1062, 1997. 16. P. Daguier, Ph. D. Thesis, Paris University, Paris, 1997. 17. M. Hinojosa, E. Bouchaud y B. Nghiem, MRS Symp. Proc. Vol. 539, p. 203-208, 1999. 18. J. Aldaco, F.J. Garza, M. Hinojosa, MRS Symp. Proc. Vol. 578, p. 351-356, 2000. 19. M. Hinojosa, J. Aldaco, U. Ortiz and J.A. González, Euromat 2000, p. 1469-1474, 2000. 20. E. Reyes, Masters Thesis, Universidad Autónoma de Nuevo León, México (in Spanish), 1999. 21. J. Schmittbuhl, J. P. Vilotte and S. Roux, Phys. Rev. E, vol. 51, p. 131, 1995.
ORAL/POSTER REFERENCE: ICF10042OR 2D SIMULATION OF TENSILE BEHAVIOR OF FIBER-REINFORCED BRITTLE MATRIX COMPOSITES WITH WEAK INTERFACE S. Ochiai, S. Kimura, M. Tanaka and M. Hojo Mesoscopic Materials Research Center, Graduate School of Engineering, Kyoto University, Sakyo-ku, Kyoto 606-8501, Japan ABSTRACT The shear lag analysis was combined with the Monte Carlo method, and applied to two-dimensional model composite to simulate the tensile behavior of unidirectional continuous fiber-reinforced brittle matrix composites with weak interface. The features particular to weakly bonded composites such as intermittent breakage of components and interfacial debonding, serrated stress-strain curve, pull-out of fibers, deleterious effects of residual stresses on strength of composite, overall unnotched strength determined by fiber bundle, longitudinal cracking arising from the tapered portion in unnotched specimen and from the notch tip in notched specimen and the notched strength given by the net stress criterion, were simulated well. KEYWORDS: tensile behavior, simulation, unidirectional composite, weak interface, damage map INTRODUCTION When the interface in brittle fiber/brittle matrix composites is strong, the crack arrest-capacity is low and high strength and toughness cannot be achieved. Thus the interface is controlled to be weak. For the design and practical use, it is needed to describe and predict the behavior of weakly bonded composites. Under tensile load, damages (breakage of fiber and matrix, and interfacial debonding) arise at many places, being distributed spatially. The damages interact mechanically to each other. Such mechanical interactions determine the species and location of the next damage, one after another. Thus the damage map and therefore the mechanical interaction among damages vary with progressing fracture. As a result of consecutive variation of them, mechanical properties such as stress-strain curve, strength and fracture morphology are determined. Thus, for description of the behavior of composites, as the damage map varies at every occurrence of new damage, the new interaction shall be calculated for the new damage map one after another One of the tools to solve this problem may be the shear lag analysis [1,2]. However, the ordinary shear lag analysis have been developed using the approximation that only fibers carry applied stress and the matrix acts only as a stress-transfer medium. Due to this approximation, it had two disadvantages; it can be applied only to polymer- and low yield stress-metal –matrix composites but not to intermetallic compound- and ceramic-matrix ones; and the residual stresses cannot be incorporated. Recently, the authors [3-5] have proposed a modified method to overcome
these disadvantages, with which the general situation (both fiber and matrix carry applied stress and also act as stress transfer media) can be described and residual stresses can be incorporated, to a first approximation. In the present study, the modified shear lag analysis mentioned above will be combined with the Monte Carlo method and be applied to 2D model, to simulate the tensile behavior of unidirectional weakly bonded brittle matrix composites. MODELING AND SIMULATION METHOD A two-dimensional model composite employed in the present work is shown in Fig.1. The components (fiber and matrix) were numbered as 1,2,…i,... to N from left to the right side. Each component was regarded to be composed of k+1 short component elements with a length Dx. The position at x=0 was numbered as 0 and then 1, 2, 3,...j... k+1 downward, in step of Dx.. The "i" component from x=(j-1) Dx to jDx was named as the (i,j)-component-element, and the interface from x=(j-1/2) Dx to (j+1/2) Dx between "i" and "i+1" components as the (i,j)-interface-element. The displacement of the (i,j)-component-element at x=jDx was defined as Ui,j. Two parameters (ai,j and gi,j) were introduced to express whether (i,j)-interface is debonded (ai,j=0) or not (ai,j=1) and whether (i,j)-component is broken (gi,j=0) or not (gi,j =1). From the spatial distribution of debonded interface-elements with ai,j=0 and broken component-elements with gi,j=0, the damage map was expressed. The values of ai,j and gi,j were determined at each occurrence of damage. The values of Ui,j were calculated by the procedure elsewhere [3,4], from which the tensile stress si,j of each component and shear stress ti,j at each interface were calculated . The simulation of the stress-strain behavior was carried out in the following procedure. (1)The strength of each component Si,j was determined by generating a random value based on the Monte Carlo procedure using the Weibull distribution. (2)Two possibilities arise for the occurrence of damage; one is the fracture of the component which occurs when the exerted tensile stress si,j exceeds the strength Si,j and another is the interfacial debonding which occurs when the exerted i j-1 j j+1 Ui,j-1 Ui,j Ui,j+1 Displacement Interface element •@•@•@(i-1,j) •@•@•@(i,j) Component element (i,j) (i,j+1) Matrix Fiber 2 j=0 •E •E •E •E j-1 j •E •E •E •E k+1 i=1 •E •E •E •E •E i•E •E •E •E •E •E•E N Dx j+1
Figure1 Modeling for simulation. shear stress ti,j exceeds the shear strength tc. To identify which occurs, si,j for all component elements were calculated and the component element having the maximum si,j/Si,j -value, say (m,n)-component , was identified. Also, the interface element with the maximum shear stress, say (m',n'), was identified. (i)If sm,n/Sm,n<1 and tm',n'/tc <1, no breakage of component and no interfacial debonding occur. Thus the applied strain was raised. (ii)If sm,n/S,m,n>=1 and tm',n'/tc <1, (m,n)-component-element is broken. If sm,n/Sm,n <1 and tm',n'/tc>=1, (m',n')-interface-element is debonded. If sm,n/Sm,n >=1 and tm',n'/tc>=1, (m,n)-component-element is broken when sm,n/S,m,n>tm',n'/tc, while (m',n')-interface- element is debonded when sm,n/Sim,n <tm',n'/tc. In this way, what kind of damage occurs is identified. Then a similar process was repeated and the next damage was identified one after another. Such a procedure was repeated until no more occurrence of damage at a given strain. (3)When no more damage occur, the applied strain ec was raised, and the procedure (2) was repeated until overall fracture of the composite. RESULTS AND DISCUSSION Figure 2 An example of the stress-strain curve and variation of damage map of the composite caused by progress of the interfacial debonding under the given geometry of the 26 ec=0.32% 13 ec=0.28% 15 14 16 17 18 19 21 22 23 24 25 20 ec=0.29% 1 2 3 4 5 6 ec=0.21% 7 9 10 8 ec=0.25% ec=0.27% Stress, 1211 Case(A) 0 100 200 300 400 500 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 s c Strain, ec ec=0% / MPa (%)
broken components. Progress of Interfacial Debonding and the Resultant Stress-strain Curve of the Composite with the Fixed Geometry of Broken Components In this work, a mini sized model composite was used to describe the fundamental process of fracture, and following values were used for calculation: N=9, k=12, df=dm=0.1mm, Dx= 0.2 mm (=2df), Ef=400GPa, Em=200GPa, Gf=160 GPa, Gm=80GPa, tc= 50- 200 MPa. and tf =0 MPa. First, in order to know the influence of pre-existent broken component element on the progress of debonding, the change of damage-map and overall stress-strain curve of composite was simulated under the given geometry of the broken components. In this case, it was assumed that no breakage of the components occurs. Figure 2 shows the stress (sc)-strain (ec) curve and damage map at various strains. The feature of the debonding progresses is read as follows. The first debonding starts at ec=0.0021, followed by the 2nd to 6th debonding at the same strain, as indicated by 1 to 6. Then the debonding stops. Due to the progress of debonding at many interface-elements, the stress-carrying capacity of composite is reduced. The reduction of stress at ec=0.0021 in the curve corresponds to such an interfacial debonding-induced loss of stress carrying capacity. After occurrence of the 1st to 6th debonding at ec =0.0021, the overall debonding stops since the shear stresses of all bonded interface-elements become lower than the critical value at this strain. Beyond ec=0.0021, no debonding occurs and the composite stress increases up to ec=0.0025, at which the 7th to 10th debonding occur one after another, resulting in loss of stress carrying capacity. After the stoppage of debonding, the composite stress again increases with increasing strain. As shown in this example, the overall debonding progresses intermittently with repetition of growth and stoppage, resulting in the serrated stress-strain curve. Stress-strain Curve of the Composite in Which Both Breakage of Components and InterfacialDebondingOccur Figure 3 shows examples of the stress-strain curve of the composite in which both breakage of components and debonding of interface occur consecutively. (a) shows the case where the coefficients of thermal expansion of fiber( af) and matrix(am) are the same (5x10-6/K) and therefore no residual stress exists and (b) the case where they are different(af =5x10-6/K and am=10x 10-6/K) and the residual stresses are introduced by cooling from 1500K to 300K. As the number of elements of broken fiber (NF), matrix (NM) and interface (NI) were quite different to each other and could not be clearly shown on the same scale, the normalized values with respect to the final values NF,f, NM,f and NI,f, respectively, are shown in this figure. Figure 4 shows the fracture process of unnotched composite with tapered grip portion. From Figs.3 and 4, following features are read. (i)The stress-strain curve shows also serration due to the intermittent breakage of the components and interfacial debonding. (ii)In case (b) in Fig.3, as af < am, the matrix and fiber have tensile and compressive residual stresses along fiber axis, respectively. In the example of Fig.3, the average failure strain of the matrix was taken to be comparable with that of fiber under no residual stresses. The authors [3,4] have shown that the tensile residual stress in the matrix enhances the breakage of matrix and also hastens the matrix breakage-induced debonding, while the compressive one in the fiber retards the fracture of fiber and also suppresses the fiber breakage-induced debonding. Comparing the variation of broken matrix-elements NM and debonded interface-elements NI with increasing applied strain in case (b) with that in case (a), the former evidently shifts to lower strain range. As known from such a difference under the existence of the present residual stresses, the matrix breakage and matrix breakage–induced debonding occur in the early stage, resulting in loss of stress carrying capacity and therefore low strength of composite. On the other hand, in the composite without residual stresses (a), the matrix is broken nearly at the same strain of breakage of fiber and the premature debonding does not arise so much until the ultimate stress, resulting in high strength. As the deleterious effect of the present residual stresses on the strength of weakly bonded composites, the strength of the composite with residual stresses (b) was 930MPa, which was far lower than 1500MPa of the composite without residual stresses (a).
(iii) In the case of (b) in Fig.3, the matrix breakage and the interfacial debonding have occurred in the early stage. Figure 3 Comparison of the stress-strain curve, strength and accumulation process of damages between the composites (a) without and (b) with residual stresses. Figure 4 (a): Schematic drawing of the longitudinal cracking in the unnotched specimen. (b): Modeling for simulation. (c) to (e): Simulated fracture process accompanied with (a) s e 0 200 400 600 800 1000 1200 1400 1600 0.0 0.2 0.4 0.6 0.8 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Stress, c /MPa Strain, c (%) F M I Without residual stress Normlized Number of Broken Elements, N F/N F,f , N M/N M,f , NI /NI,f (b) 0 200 400 600 800 1000 0.0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 Stress, s c /MPa Strain, ec (%) F M I With residual stress Normlized Number of Broken Elements, N F/N F,f, NM/N M,f, N I /NI ,f Longitudinal cracking (a) Fiber Matrix (b) (c) (d) (e)
longitudinal cracking and the fracture morphology Once such a situation has occurred, the fibers are separated to each other and behave like a fiber bundle without matrix. Thus the strength of such composite is given by the strength of the fiber bundle. Further simulation using the low failure strain-matrix composite without residual stresses showed the same feature. (iv) It has been known that, in weakly bonded composites, longitudinal cracking occurs ahead of the notch. Furthermore, even in the unnotched samples, the longitudinal cracking between the parallel gage portion and the tapered one occurs (Fig.4(a)). Such a feature is also well realized in Fig.4 (b) to (e). (v) It has been well known that the fiber is pulled-out in weakly bonded composites. The final fracture morphology of the composite, shown in Fig.4(e), describes such a feature also. (vi) In notched samples, the longitudinal cracking occurs in the whole length between the grips at lower applied stress level than that in unnotched samples. Thus, the notch does not cause mode I type fracture but enhance longitudinal cracking. As a result, the strength of notched samples could be expressed by the net strength criterion. CONCLUSIONS The variation of the damage map, stress-strain curve, strength and fracture morphology of two-dimensional model composite with weak interface were simulated by combining the shear lag analysis with the Monte Carlo method. Following features of weakly bonded composites were described. (1) Both of the breakage of components and interfacial debonding occur intermittently. (2) As a result of (1), the stress-strain curve is serrated. (3) When the fracture strain of matrix is low, the residual stresses (tensile and compressive stresses along fiber axis for matrix and fiber, respectively) enhance breakage of matrix and matrix breakage-induced debonding at low applied strain. As a result, the stress carrying capacity of the composite is reduced, resulting in low strength. (4) The strength of weakly bonded composites whose matrix has low failure strain is practically given by the strength of the fiber bundle. (5) Longitudinal cracking arises at the notch tip in notched specimens and also at the tapered corner in the unnotched specimens. (6) The notched strength is given by the net stress criterion. Acknowledgement The authors wish to express their gratitude to The Ministry of Education, Science and Culture of Japan for the grant-in-aid for Scientific Research (No.11555175). REFERENCES 1. Hedgepeth, J. M.. (1961). Stress Concentrations in Filamentary Structures. NASA TN D-882, Washington,. 2. Oh, K. P. (1979) J. Comp. Mater., 13, 311. 3. Ochiai, S., Okumura, I., Tanaka, M. and Inoue, T. (1998) Comp. Interfaces, 5, 363. 4. Ochiai, S., Tanaka, M. and Hojo, M. (1998) Comp. Interfaces, 5, 437. 5. Ochiai, S. Hojo, M. and T. Inoue. (1999) Comp. Sci. Tech., 59, 77.
ICF100844OR A BOUNDARY ELEMENT BASED MESO-ANALYSIS ON THE EVOLUTION OF MATERIAL DAMAGE H. Okada, Y. Fukui and N. Kumazawa Department of Mechanical Engineering, Kagoshima University 1-21-40 Korimoto, Kagoshima 890-0065, JAPAN ABSTRACT In this paper, we present efficient numerical formulations for the analyses of particulate composites, underging meso-structural changes such as particle fracture, stress induced phase transformation, etc. The formulations are derived based on the homogenization method and the boundary element method (BEM). Proposed formulations can efficiently deal with problems, in which particles randomly distribute and orient in the composites, by combining analytical solutions for the ellipsoidal inclusions such as Eshelby’s tensor and the boundary element method. Hence, there is no need to carry out any numerical integrations for the particles. Proposed numerical methods are computationally efficient and accurate. The formulations and numerical results for effective elastic moduli of composites and problems of particle fracture and stress induced phase transformation, are presented. INTRODUCTION In this paper, efficient boundary element formulations (see [1-6] for the boundary element method) for solids containing ellipsoidal inclusions/particles, as shown in Fig. 1, are presented (see Ashby [7] for various types of composite materials). First, a boundary element based formulation for homogenization analysis is discussed and is combined with analytical solutions for ellipsoidal inclusions in which constant initial strains are specified inside of them. A special case of the analytical solutions is given as the well known Eshelby’s tensor [8,9]. The analytical solutions for ellipsoidal inclusions such as Eshelby’s tensor, are based on the fundamental solution of linear isotropic elasticity, which is also used in the boundary element method as its kernel functions. The analytic solutions and boundary element formulation can easily be combined. Homogenization method [10-12] based on the finite element method has been applied to a various class of problems, such as identifying macroscopic elastic moduli and nonlinear behavior of meso-structure for a prescribed macroscopic deformation mode. Homogenization method assumes that the microstructure of solid is spatially periodic, and finite element analyses for a unit of periodic structure (unit cell) are carried out. However, the homogenization method can not be applied to the problems of particulate composite materials in a straight forward manner, since the orientations and distributions of particles Densely Distribute Particles Figure 1: Particulate composite material
would be somewhat random. Therefore, defining a unit cell model containing one or a few particles may be an over simplification of the problem. To accurately model such solids, a unit cell containing many particles should be analyzed. Thus, one needs to build and carry out analyses for a unit cell model having tens and hundreds of particles. Finite element model, which is required for such analyses, would be gigantic. Generating the finite element model as well as solving the problem would pose many problems. However, mechanical interactions between the particles, and between matrix material and particles are fully accounted for, when the finite element method is adopted. On the other hand, methodologies in micromechanics, such as self-consistent method [9,13,14] and Mori-Tanaka theory [9,15] have been presented. For particulate composites, Eshelby’s tensor [8,9], takes a central role (see, for example, [9,16,17]). These methodologies have advantages over the finite element method such that the effective mechanical behaviors of composites can be estimated in a closed or semi-closed form, based on the elastic moduli of matrix and particles, the distribution, orientation and volume fraction of the particles. Large scale computations are not required. However the methods in micromechanics do not account for detailed mechanical interactions between the distributed particles. The mechanical interactions may have significant roles when the particles are densely distributed or when we attempt to account for the damage evolution or meso-structural changes, such as phase transformation of the particles. Boundary element based homogenization formulation, which is developed in this paper, have advantages of both the above mentioned methodologies. Since the method is based on the boundary element method (BEM [1-6]), detailed mechanical interactions between all the material constituents can be accounted for. This nature is similar to that of the finite element method. A unit cell modeled by the boundary element method can contain many particles and the computation is simplified by using analytical solutions for ellipsoidal inclusions, such as Eshelby’s tensor [8,9]. In this regard, proposed method is similar to the methodologies in micromechanics. HOMOGENIZATION FORMULATION FOR PARTICULATE COMPOSITES In this section, equation formulations for boundary element based homogenization analysis for particulate composites are briefly discussed (see [18,19] for homogenization method based on BEM). In homogenization method, the microstructure of solid is assumed to be spatially periodic, as shown in Fig. 2. A unit of the periodic microstructure is modeled by the boundary element method and is called “unit cell”. As shown in Fig.1 in a two dimensional illustration, the size of the unit cell is represented by e . We introduce two different coordinate systems. One is ix coordinate system, which is fixed in space and the other is iy coordinate system, which is scaled by the size e of unit cell, as: i i i c y x = +e (1) where ic are the components of an arbitrary vector. Displacements iu are expressed by the two scale asymptotic expansion, by following Guedes and Kikuchi [10], as: ( ) ( ) ( ) ( ) ( ) L+ + + + = x, y x,y x,y x,y x, y 3 3 2 2 1 i i i o i i u u u u u e e e (2) For small deformation linear elasticity problem, Hooke’s law and the equation of linear momentum balance are written to be: l l x u E k ijk ij ¶ ¶ s = , 0 + = j i ij b x¶ ¶s (3) e y 1 y 2 Unit cell representing a unit of periodic microstructure Global solid composed of N unit cells (a) (b) Y M Y * Figure 2: Two dimensional illustration for a solid whose microstructure is periodic and its unit cell.
where l ijk E are the fourth order tensor, representing Hooke’s law. Based on Eqns. (2) and (3), it can be shown that ( ) u u x o i o i = [ o iu are the functions of ix only.], and one can obtain an integral equation formulation for 1 iu for the analysis of unit cell (see Okada et al. [19] for the derivation). ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] Y pq q m Y kp m m o k j Y ip m k Mijk k Y kp m Y Y jp m Y j Y pq q m y Y C t x u Y y u u Y E t Y t u C u x ¶ x ¶ ¶ ¶ x ¶ e ¶ x x ¶ x ¶ ¶ ¶ + - - - = ò ò ò ò Y * Y * * Y 1 * * 1 d d d d * y, y, y, y, l l (4) where lke are the fictitious initial strains, which are defined by, ( ) ( ) ( ) ïþ ï ý ü ïî ï í ì + - = l l l l y u x u C E E k o k M mnk mnk Mijmn ij ¶ ¶ ¶ ¶ e x, y x 1 * (5) Displacements 1 iu at the source point Y mx is evaluated by the integral equation (4). A two phase composite material is assumed and elastic constants for matrix and second phase materials are represented by Mijk E l and * l ijk E , respectively. * ip u and * ip t are the Kelvin solutions (see [4]). Y and Y¶ represent the domain and boundary of the unit cell. *Y denotes the region of second phase material in Y. We assume that the solid contains ellipsoidal particles as its second phase material and that the stresses and strains are constant values inside a particle, by following the result of Eshelby [8] and many of micromechanics analyses [9,17]. We then rewrite the volume integral term of integral equation (4), as: ( ) ( ) ( ) ( ) ( ) å ò å ò = = = L = N I I ij Y m I pij N I I Mijk j Y ip m I ij k Mijk j Y ip m I E Y y u Y E y u 1 Y 1 Y * * * * * * d d x e ¶ ¶ x e e ¶ ¶ x l l l y, y, (6) where ( ) Y m I pij x L are the analytical expressions for the integral (see Mura[9]) and I kl e are the fictitious initial strain in the Ith particle. Thus, we write: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] Y pq q m Y kp m m o k N I I ij Y m I pij k Y kp m Y Y jp m Y j Y pq q m y Y C t x u u Y t Y t u C u x ¶ x ¶ ¶ x e ¶ x x ¶ x ¶ ¶ ¶ + - - L - = ò å ò ò = Y * 1 Y 1 * * 1 d d d y, x y, y, (7) When the source point Y ix is at the interior of the Jth particle, we can derive an integral equation for the displacement gradients j iu y ¶ ¶ ˆ , by differentiating Eqn. (7) with respect to Y ix , as: ( ) ( ) ( ) ( ) ( ) J k J pqk N I j I I ij Y m I pqij Y k q Y kp m Y q Y Y jp m j J q p S u Y t Y u t y u l l e x e ¶ ¶x ¶ x ¶ ¶x x ¶ ¶ ¶ ¶ ¶ - - Q - = ÷ ÷ ø ö ç ç è æ å ò ò ¹ = * * 1 ˆ d , d , ˆ y y (8) where J pqk S l are the components of Eshelby’s tensor [8,9] for the Jth particle. Singular volume integrals, whose numerical evaluations are known to be troublesome, are replaced by the analytical formulae such as Eshelby’s tensor. Therefore, proposed integral equations are computationally efficient and accurate. Taking the last term in Eqn. (7) as the forcing term and following the standard boundary element analysis procedures by imposing the so called periodic boundary conditions on displacements 1 iu and tractions Y it , we can evaluate the displacements and tractions at the boundary of the unit cell. An initial strain iteration method is adopted to obtain the equilibrium (see [5,6] for the initial strain iteration for solving elastic-plastic problems using BEM). Thus, the solutions for given j o iu x ¶ ¶ are obtained. We repeat the analyses for six times for a three dimensional problem to determine the responses of microstructure to all the macroscopic deformation modes (six strain modes). The characteristic functions ( )y l ikF , which relate j o iu x ¶ ¶ to 1 iu , and the effective elastic moduli Hijk E l of the composite are written to be:
( ) ( ) ( ) ( ) ( ) ï ï î ï ï í ì ¶ ¶ = ÷÷ ø ö çç è æ - + - = + å å = * = * l l l l l l l l l x u u F Y y F E E Y E E Y Y E E o k ik i N I I I n mk Mijk ijk N I I Mijk ijk Mijk Hijk y y x 1 1 * 1 * 1 1 ¶ ¶ (10) ANALYSIS FOR PROGRESSIVE DAMAGE (MESO-STRACTURAL CHANGE) Here, we deal with the problems of particle fracture and of the stress induced phase transformation of particles (see [17,20] for stress induced phase transformation). We make small extensions on the integral equations (7) and (8), as follows. Particle Fracture An incremental analysis is carried out for the problems of particle fracture. The simplest scenario is assumed that when the stresses in a particle satisfy a criterion for particle fracture, the elastic modulus of the particle reduces to zero (in actual calculation, the elastic modulus is reduced to be 1/1000 of the original value). Therefore, the analysis is entirely based on the elastic analysis for the unit cell. The integral equations, which are developed in the previous section, are applied by specifying different elastic moduli for fractured and unfractured particles. Algorithms for the analysis are shown in Fig. 3. Stress Induced Phase Transformation of Particles The problems of dilatational stress induced phase transformation (see Okada et al. [17]) are considered and an incremental algorithm is adopted. It is assumed that at the instance, when hydrostatic stress ( 3 kk s ) inside a particle reaches a critical value, the particle transforms its phase and dilatational transformation strain is produced. The dilatational transformation strain is modeled as an additional initial strain. The integral equations [Eqns. (7) and (8)] are modified to include the effects of the transformation strain, as: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] Y pq q m Y kp m m o k N I I ij Y m I pij I N I I ij Y m I pij k Y kp m Y Y jp m Y j Y pq q m y Y C t x u u Y t Y t u C u x ¶ x ¶ ¶ x e b x e ¶ x x ¶ x ¶ ¶ ¶ + - + L - L - = ò å å ò ò = = Y * 1 1 Y 1 * * 1 d d d y, x y, y, (11) ( ) ( ) ( ) ( ) ( ) ( ) J ik J pqk J N I j I I ij Y m I pqij I J ik J pqk N I j I I ij Y m I pqij Y k q Y kp m Y q Y Y jp m j J q p S S u Y t Y u t y u l l l l e x e b b e x e ¶ ¶x ¶ x ¶ ¶x ¶ x ¶ ¶ ¶ ¶ + + Q - - Q - = ÷ ÷ ø ö ç ç è æ å å ò ò ¹ = ¹ = * * 1 1 ˆ d , d , ˆ y y (12) where I ij e are fictitious initial strains, which are related to the dilatational transformation strain t e , as: Start (all the particles have original properties): set load factor a=0 Homogenization analysis for given ¶ui o ¶x j Updateload factor a to satisfy the particle fracture criteron in a particle (K=1,2,3,...,N) Evaluate macroscopic stresses Reducing the elastic moduli of the fractured particle by 1/1000 Evaluate macroscopic stresses 1st step a´s ij ¶um o ¶x n ( ) a´s ij ¶um o ¶xn ( ) f a´sij K ¶u m o ¶x n ( ) ( = C Check if: (excluding fractured ones) Figure 3: Algorithms for the analysis of problem of particle fracture Start (all the particles are not transformed): set load factor a=0 Homogenization analysis for given ¶ui o ¶x j Updateload factor a to satisfy the phase transformation criteron in a particle (K=1,2,3,...,N) Evaluate macroscopic stresses Assign additional trasnformation strains to the transformed particle Evaluate macroscopic stresses 1st step a´s ij ¶um o ¶x n ( ) a´s ij ¶um o ¶x n ( ) a´ 1 3 sii K ¶um o ¶x n ( ) = sT Check if: (excluding transformed ones) Figure 4: Algorithms for the analysis of problem of stress induced phase transformation
÷ ø ö ç è æ = t k mm Mijk I ij C E e e 3 1 * l l (13) The effective stress for the solid can be given as average values of stresses in the unit cell, as: ( ) å å = = ÷ ø ö ç è æ ÷ - ø ö ç è æ - + = N I pt ijkk I N I I k Mijk ijk o M k ijk ij E u E E x Y u E 1 * 1 * 3 1 ˆ 1 e b ¶ ¶ ¶ ¶ s l l l l l (14) Algorithms for the analysis of the stress induced phase transformation problem are given in Fig. 4. RESULTS OF NUMERICAL ANALYSIS A limited number of numerical results are presented in this section due to the restriction of pages (results for the problem of particle fracture are omitted.). Unit cell models used for numerical computations are shown in Fig. 5. Randomly distributed spherical particles are assumed. The volume fraction of the particles is about 10%, and the number of the particles are 27 and 125. Linear quadrilateral boundary elements are used. Total number of boundary elements on each face of the unit cell is 36 (6 x 6). Evaluation for Effective Elastic Moduli The results are shown in Fig. 6. Isotropic elasticity is assumed for the particles and the ratio of Young’s moduli of matrix and particles are varied from 10-3 to 103. Poisson’s ratio is assumed to be 0.3 for both the material constituents. Results, which are analyzed by the 125 particle model, are presented. For a comparison purpose, the results obtained by using self-consistent method and Eshelby’s method theory (see Mura [9]) are also plotted in the figures. The results by three different methods are within an agreement. Though exact solutions are not known, three different methods are within a reasonable agreement and, therefore, the proposed technique is, at least, proven to be reliable. Stress Induced Dilatational Transformation Relationships between the effective hydrostatic stress ( 3 kk s ) and effective dilatational strain ( i o iu x ¶ ¶ ), when Young’s moduli of matrix and of particles are set to be the same, are shown in Fig. 7. Dilatational transformation strain is assumed to be 0.05 = t kk e and the phase transformation is assumed to take place when the hydrostatic stress in a particle reaches transformation stress Ts . The elastic moduli for matrix and particles are set to be the same in this case. The stress-strain relationships follow zigzag paths and have negative slopes while the phase transformation is undergoing. The path is smoother for the 125 particle model than for the 27 particle one. 7.0E-01 8.0E-01 9.0E-01 1.0E+00 1.1E+00 1.2E+00 1.3E+00 1E-03 1E-02 1E-011E+001E+011E+021E+03 7.0E-01 8.0E-01 9.0E-01 1.0E+00 1.1E+00 1.2E+00 1E-03 1E-02 1E-011E+001E+011E+021E+03 10 1 -1 10 10 -2 -3 10 10 2 10 3 1 1.1 1.2 0.9 0.8 0.7 E * / EM K / Ko 10 1 -1 10 10 -2 -3 10 10 2 10 3 1 1.1 1.2 0.9 0.8 0.7 E * / EM G / Go 1.3 : Present BEM Homogenization : Self-consistent method [9] : Eshelby's method [9] : Present BEM Homogenization : Self-consistent method [9] : Eshelby's method [9] (a) (b) Figure 6: The results for effective elastic moduli (125 particles). (a) Bulk modulus, (b) Shear modulus. (a) (b) Figure 5: Distributions of spherical particles in the unit cell. (a) 27 particles whose volume fraction is 0.1, (b) 125 particles whose volume fraction is 0.089
CONCLUDING REMARKS In this paper, a new but simple method for the analysis of particulate composite material is presented. Though the numerical results, which are presented here, are rather limited, the method is proven to be quite effective. If one carried out a three dimensional analysis with 125 randomly distributed particles in a unit cell using the finite element method, a large scale computation must be carried out and the state of art mesh generation software would be necessary for the generation of analysis model. All the numerical analyses, which are presented in this paper are carried out using a workstation within a reasonable amount of computational time. Therefore, it can be concluded that present numerical technique can deal with the meso-mechanics problems of particulate composites effectively and efficiently. REFERENCES 1. Cruse, T. A. (1969), Int. J. Solid Structures 5, 1259. 2. Cruse, T. A. (1987), Comp. Meth. Appl. Mech. Engrg 62, 227.-244, (1987). 3. Rizzo, F. J. (1967), Q. Appl. Math. 25, 83. 4. Banerjee, P. K. and Butterfield, R. (1981), Boundary Element Methods in Engineering Science, McGraw-Hill, London. 5. Okada, H. and Atluri, S.N. (1994), Solids and Structures 31, 12/13, 1735. 6. Okada, H., Rajiyah, H. and Atluri, S.N. (1988), Computers & Structures 30, 1/2, 275. 7. Ashby, M. F. (1993), Acta. Metal. Mater. 41, 5, 1313. 8. Eshelby, J. D. (1957), Proceedings of Royal Society, London A 241, 376. 9. Mura, T (1982), Micromechanics of Defects of Solid, Martinus Nijhoff. 10. Guedes, J. M. and Kikuchi, N. (1990), Comp. Meth. Appl. Mech. Engrg. 83, 143. 11. Kalamkarov, A. L. (1992), Composite and Reinforced Elements of Construction, John Wiley & Sons. 12. Hassani, B. and Hinton, E. (1998), Computers & Structures 69, 707. 13. Christensen, R. M. (1990), J. Mech. Phys. Solids 28, 3, 379. 14. Hill, R., J. Mech. Phys. Solids 13, 180. 15. Mori, T. and Tanaka, K. (1973), Acta Metal. 21, 571. 16. Nakagaki, M., Wu, Y. and Brust, F. W. (1999), Computer Modeling and Simulation in Engineering 4, 3, 186. 17. Okada, H. Tamura, T., Ramakrishnan, N., Atluri, S. N. and Epstein, J. S. (1992), Acta. Metall. Materi. 40, 6, 1421. 18. Kaminski, M. (1999), Engineering Analysis with Boundary Elements 23, 815. 19. Okada, H., Fukui, Y. and Kumazawa, N., In Advances in Computational Engineering & Sciences, pp. 1128-1133, Atluri, S. N. and Brust, F. W. (Eds.), Tech Science Press． 20. Ramakrishnan, N., Okada, H. and Atluri, S. N. (1991), Acta. Metall. Mater 39, 6, 1297. (1.0) (0.5) 0.0 0.5 1.0 1.5 2.0 0E+00 2E-03 4E-03 6E-03 0 2 4 6 （×10 ） -3 Dilatational Strain e kk Stress (s / 3) / s kk T -0.5 0.0 0.5 1.0 1.5 2.0 -1.0 (1.0) (0.5) 0.0 0.5 1.0 1.5 2.0 (4E-22) 2E-03 4E-03 6E-03 0 2 4 6 （×10 ） -3 Dilatational Strain e kk Stress (s / 3) / s kk T -0.5 0.0 0.5 1.0 1.5 2.0 -1.0 (a) (b) Figure 7: Hydrostatic pressure stress-dilatational strain curves for the problems of stress induced phase transformation. (a) 27 spherical particles, (b) 125 spherical particles. The straight lines are stress-strain curves following the results of Ramakrishnan, Okada and Atluri [20].
www.icfweb.orgRkJQdWJsaXNoZXIy MjM0NDE=