- Email: [email protected]
Micromechanical finite element analysis of metal matrix composites using nonlocal ductile failure models
Computational Materials Science 37 (2006) 29–36 www.elsevier.com/locate/commatsci
Micromechanical finite element analysis of metal matrix composites using nonlocal ductile failure models T. Drabek 1, H.J. Bo¨hm
*
Christian Doppler Laboratory for Functionally Oriented Materials Design, Institute of Lightweight Design and Structural Biomechanics, Vienna University of Technology, A-1040 Vienna, Austria
Abstract Finite element studies of ductile damage in the matrix of particle and fiber reinforced metal matrix composites are presented. Three material models capable of supporting such simulations at the constituent level are discussed within a unified framework. They comprise an element removal technique triggered by a ductile damage indicator as well as versions of the ductile rupture models of Gurson and Rousselier. Nonlocal averaging is employed for reducing the mesh dependence typically displayed by continuum damage methods upon the onset of local softening. Implementation issues specific to the use of nonlocal damage models in a continuum micromechanics framework are discussed. The efficacy of the approach in limiting the mesh sensitivity of the predicted behavior of composites subject to matrix damage is demonstrated for a simple three-dimensional matrix-particle configuration. Applications of the method to multi-particle and multi-fiber unit cells subjected to uniaxial tensile loads are presented and effects of microgeometrical parameters on the mechanical response are discussed. 2005 Elsevier B.V. All rights reserved. PACS: 46.30.Nz; 83.70.Dk; 83.10.Ff; 02.70.Dh; 11.10.Lm Keywords: Ductile damage models; Nonlocal averaging; Metal matrix composites; Micromechanics
1. Introduction In metal matrix composites (MMCs) and related materials damage and failure are caused by debonding of the matrix-reinforcement interface, by brittle fracture of the reinforcing fibers or particles, and/or by ductile failure of the matrix. The present work concentrates on the development and application of algorithms that are capable of modeling the latter failure mode. Other damage mechanisms are not accounted for.
* Corresponding author. Tel.: +43 1 58801 31712; fax: +43 1 58801 31799. E-mail address: [email protected] (H.J. Bo¨hm). URL: http://www.ilsb.tuwien.ac.at (H.J. Bo¨hm). 1 Present address: Voith–Siemens Power Generation, A-3100 St. Po¨lten, Austria.
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.12.032
Ductile damage and failure of metals are caused by the growth and coalescence of existing voids and by the nucleation of new voids by debonding from tiny inclusions such as precipitates. A number of ductile damage models have been proposed in the literature for describing the behavior of metals undergoing this damage process. The implementation of such models into a finite element code supports the investigation of the failure of ductile constituents of inhomogeneous materials under given macroscopic loading conditions at length scales where the phases are geometrically resolved, i.e., on micro- and meso-scales. Such modeling strategies are referred to as continuum micromechanical methods. They allow the overall behavior of composites undergoing ductile damage and failure of the matrix to be obtained by homogenization. Standard ductile damage models have been combined with unit cells for continuum micromechanical studies of MMCs by a number of
30
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
authors, the majority of whom used axisymmetric cells, simple three-dimensional models or planar multi-particle geometries, see, e.g., [1–3]. It is well known that in the softening regime ductile damage models (like continuum damage models in general) tend to give rise to a marked mesh dependence of the predicted results. This behavior is caused by the loss of ellipticity of the governing partial differential equations once softening sets in. In order to overcome this mesh sensitivity the present study applies a nonlocal averaging approach to three different ductile damage models. This procedure introduces an additional material parameter in the form of a characteristic length. An important precondition for the use of the above modeling strategy, viz., that the size of individual voids is much smaller than the width of the ‘‘matrix bridges’’ between neighboring particles or fibers, is usually fulfilled for micromechanical studies of MMCs and related materials. Ductile damage models that employ nonlocal averaging were introduced into continuum micromechanics by Tvergaard and Needleman [4] and have recently seen further development [5,6]. In the following, nonlocal ductile damage models are described that are specifically geared towards the requirements posed by micromechanical studies of the local and overall response of metal matrix composites that undergo ductile damage. 2. Implemented ductile damage models A number of ductile damage models can be found in the literature that are suitable for describing ductile failure processes via appropriate constitutive equations that account for void growth. In these models the yield function is typically treated in extension of J2-plasticity theory, so that expressions of the type U¼
req rf ¼ 0 D
or
b ¼ req D ¼ 0 U rf
ð1Þ
may be used to describe the yield condition. Here req stands for the von Mises equivalent stress, rf is the flow stress, and D is a parameter that indicates the progress of damage and varies between 1 and 0. A completely damage-free material (i.e., one in which no voids are present) is indicated by D ¼ 1 and follows standard J2-plasticity theory. When D reaches zero, ductile failure occurs and the von Mises stress req vanishes at the same time, so that the term req =D in Eq. (1) remains finite. Note that the above definitions differ from the ones commonly used for damage parameters in continuum damage mechanics. Formally, the ductile damage models discussed in this section differ only in the definition of D, with the mean stress rm (and, consequently, the stress triaxiality g = rm/ req) playing an important role in determining D in all three models. For this reason these models can be classified as pressure dependent plasticity models, in which the first invariant of the stress tensor must be considered in order to obtain a closed yield surface in stress space.
2.1. Ductile damage indicator triggered (DDIT) models Gunawardena et al. [7] proposed a ductile damage indicator that is capable of detecting local ductile failure of a given material point by postprocessing operations following a standard elastoplastic finite element analysis. This method, which is based on work of Rice and Tracey [8] as well as Hancock and Mackenzie [9], was discussed in depth by Fischer et al. [10]. The increment of the ductile damage indicator D can be written as exp 32 g Deeq;p ; ð2Þ DD ¼ 1:65e0 where eeq,p stands for the accumulated equivalent plastic strain and the material parameter e0 is a reference failure strain that can be calibrated by appropriate tensile tests [10]. Regions where local ductile failure has occurred are indicated by the accumulated value of D reaching or exceeding a critical value Dc, which is usually set to one. The ductile damage indicator, however, does not uniquely correlate with the degradation of the mechanical properties due to evolving damage in the range 0 < D < Dc and per se does not introduce any loss of stiffness or strength into the material model. In order to achieve a local softening response in a finite element context the ductile damage indicator can be used to trigger either the removal of the element or a marked reduction of load carrying capacity and stiffness (such that req ! 0 and rm ! 0) at the integration point(s) involved. For D < Dc the local material behavior continues to follow J2 flow theory in such a model. Accordingly, D as defined in Eq. (1) takes the form 1 for D < Dc ; D¼ ð3Þ 0 for D P Dc in descriptions of this type, which are referred to as ductile damage indicator triggered (DDIT) models in the following. 2.2. Ductile rupture models Another group of approaches to describing void growth and coalescence in ductile materials are based on continuous modifications of the von Mises flow potential in dependence on the mean stress, rm, and the volume fraction of voids, f. The best known continuum model of the above type was developed by Gurson [11] and extended by Needleman and Tvergaard [12] to include effects of void nucleation and void coalescence. For the resulting GTN model the parameter D according to Eq. (1) can be denoted as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3q2 rm ð4Þ D ¼ 2q1 f cosh þ 1 þ q3 f 2 . 2rf where q1, q2 as well as q3 are adjustable parameters. The void volume fraction at which ductile failure takes place
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
can be evaluated by rearranging Eq. (4) at D ¼ req ¼ rm ¼ 0 to give ff ¼
1 q1
for q3 ¼ q21 .
ð5Þ
An alternative ductile rupture model was proposed in the form of a thermodynamically based continuum damage mechanics method by Rousselier [13] and extended by Aboutayeb [14], with D taking the form rm e D ¼ 1 q4 f D exp ; ð6Þ q4 rf e are material parameters. An imporwhere q4 as well as D tant feature of the extended Rousselier model [14] lies in the fact that it allows to calculate the void volume fraction at failure as ff ¼
1
ð7Þ
e q4 D
in dependence of the material parameters only. The increment of the void volume fraction can be calculated in both the Gurson and Rousselier models by evolution equations of the type Df ¼ ð1 f ÞBDep ;
ð8Þ
where Dep stands for the increment of the first invariant of the tensor of plastic strains and B is defined as 1 for f 6 fc ; B¼ ð9Þ > 1 for f > fc . Eqs. (8) and (9) allow the void volume fraction to be integrated up starting from initial values f0. To avoid a discontinuity in the history of Df in Eq. (8), B can be chosen to be a linear function of the void volume fraction, Bðf Þ ¼ 1 þ
Bf 1 ðf fc Þ; ff fc
ð10Þ
for f > fc. Here Bf = B(f = ff) represents the value of B at total failure of the material point. Additional details of the models are given in [6,15]. The above three models for describing ductile damage were implemented into the general purpose Finite Element code ABAQUS/Standard V6 (Abaqus Inc., Pawtucket, RI, 2001) as user defined material subroutines (UMATs). For the DDIT model an implementation via a user defined field subroutine (USDFLD) was developed in addition. In the UMATs the consistent tangent (Jacobian) required for evaluating the stiffness matrix of the model, Ccons ¼
oDr or ¼ oDe oe
ð11Þ
was derived by the procedures proposed by Aravas [16] and extended by Zhang [17] for pressure dependent plasticity models. Predictions obtained by the basic (‘‘local’’) versions of all of the above models typically show a clear mesh depen-
31
dence once the acoustic tensor defined by the double contraction Qij ¼ C ijkl nk nl
ð12Þ
of the material tensor C and the normal vector n, which indicates possible wave propagation directions, becomes negative definite, see, e.g., the discussion by Baaser and Gross [18]. This mesh sensitivity, on the one hand, takes the form of a marked dependence of the crack path on the employed finite element mesh, the assemblage of failed elements making up the crack typically having a width of one element. On the other hand, the resulting overall force–displacement responses tend to vary in such a way that the energy required for failure decreases with increasing mesh refinement. As a consequence, the advantageous convergence properties of the Finite Element method upon mesh refinement are no longer guaranteed. These undesirable effects tend to be less marked in DDIT schemes than in ductile rupture models. 3. Nonlocal regularization technique The mesh sensitivity of ductile damage algorithms can be markedly reduced by applying regularization techniques. In such approaches an additional material parameter is introduced in the form of a characteristic length 2L. Physically this internal length scale of the material may be interpreted as the size of a process zone or as a function of the distance between individual voids or the distance between dimples in the fracture surfaces [4], but is not fully understood at present. The regularization method employed in the present work is based on smoothing the rate of an appropriate damage variable q by nonlocal averaging [19], a delocalization function proposed by Leblond [20] being employed. This nonlocal averaging procedure can be denoted in terms of the increment of the nonlocal damage variable, DqNL, as 2 32 DqNL ðxi Þ ¼
n 6 1 X 6 DqL ðb x j Þ6 4 W ðxi Þ j¼1
1þ
7 7 8 7 D Vb j ; 5 jxi b xjj 1
ð13Þ
L
the normalizing factor W being defined as 2 32 W ðxi Þ ¼
n 6 X 6 6 4 j¼1
7 7 8 7 D Vb j . 5 jxi b xjj
1þ
1
ð14Þ
L
Here qL is the local damage variable and qNL is its nonlocal counterpart, jxi b x j j stands for the distance between two integration points i and j, D Vb j is the volume associated with integration point j, and the sums run over all integration points j within a distance L from a given integration point i. The variables subjected to smoothing are the
32
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
damage indicator D in the case of the DDIT scheme and the void volume fraction f in the case of the GTN and Rousselier models. The nonlocal smoothing routines employed in the present work employ lists of neighbors for each integration point to allow efficient averaging by Eqs. (13) and (14). These lists are generated in a preprocessor step and can be updated in the course of an analysis to account for the effects of large strains. Nonlocal averaging is applied only within one constituent, viz., the metallic matrix. This feature is physically reasonable for composites consisting of a ductile and an elastic phase, but may have to be amended in cases where more than one constituent can be expected to be subject to ductile damage. Integration points where failure has occurred are excluded from the sums in Eqs. (13) and (14). Provision is also made for consistently handling nonlocal averaging in the presence of the periodicity and/or symmetry boundary conditions that are typically encountered in problems involving periodic homogenization. Fig. 1 schematically depicts the region of influence of an integration point (marked by an ‘‘x’’) that is located near the border of a two-dimensional unit cell employing periodicity boundary conditions. The user subroutine interfaces provided by ABAQUS/ Standard allow nonlocal averaging to be carried out only at the ends of increments, where the UEXTERNALDB routine can be invoked, but not at the level of the equilibrium iterations. Comparisons with published results [21] showed that the resulting modification to the ‘‘standard’’ smoothing algorithms gives rise only to minor changes in the predicted crack paths and overall stress–strain
responses [22]. Details of the averaging procedure are discussed in [6,15]. It is worth noting that the use of local damage models may be interpreted as being equivalent to implicitly specifying a distribution of characteristic lengths via the element sizes. For simple microgeometries, a judicious choice of the discretization can accordingly give rise to valid results in the softening regime. In the case of geometrically complex multi-fiber or multi-particle models, however, where the phase geometry enforces rather irregular meshes with local element sizes that are driven by the requirements of resolving the stress and strain fields, local damage models may introduce considerable inconsistencies into the material descriptions. Finally, it must be emphasized that regularization by nonlocal averaging is viable only if the mesh size in regions where ductile damage may occur is smaller than (or at most equal to) the characteristic length 2L. This condition is usually met by micromechanical studies, but may be very difficult to fulfill for macroscopic analyses of components or structures. 4. Application to micromechanical studies of MMCs A number of applications of local ductile damage models to micromechanical studies of ductile matrix composites have been reported in the literature, where the Gurson ductile rupture model, see, e.g., [1,3,23,24], and DDIT-like schemes, see, e.g., [2,25,26], were employed. Recent work of this type has also involved multi-particle unit cells [27,28]. The use of ductile damage models with nonlocal averaging in the context of continuum micromechanics was initiated by Tvergaard and Needleman [4]. Recently such approaches have been further developed for use with more complex microgeometries [5,6,22,29]. In the following the nonlocal damage routines presented in Sections 2 and 3 are employed for studying the effects of ductile damage on the mechanical behavior of particle and fiber reinforced MMCs. The material properties of the matrix correspond to Al2618-T4 and those of the fibrous or particulate reinforcements to SiC. The corresponding material parameters are listed in Table 1, where E stands for the Young’s modulus and m for the Poisson number. J2 plasticity with a modified Ludwik hardening law was used to describe the damage-free behavior of the matrix, with ry0 standing for the initial yield stress, m for the hardening modulus, and n for the hardening exponent. Results obtained with two-dimensional and threedimensional unit cell models are presented, which pertain
Table 1 Material parameters used for the elastoplastic Al2618-T4 matrix (modified Ludwik hardening law) and the elastic SiC reinforcements Fig. 1. Sketch of regions where Eqs. (13) and (14) are applied for nonlocal averaging in a finite element mesh employing periodic boundary conditions.
Al2618-T4 matrix SiC reinforcement
E (GPa)
m [1]
ry0 (MPa)
m (MPa)
n [1]
70 450
0.30 0.17
184 –
722.7 –
0.49 –
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
q4 [1]
e [1] D
f0 [1]
fc [1]
Bf [1]
1.0
2.0
0.004–0.008 (random)
0.2
3.0
to the uniaxial mechanical response of composites reinforced by continuous aligned fibers and by equiaxed particles, respectively. The investigation of a fiber reinforced MMC, Section 4.2, was based on the nonlocal DDIT scheme, the reference strain being chosen as e0 = 0.2. A nonlocal extended Rousselier model was used in the studies involving reinforcement by particles, Sections 4.1 and 4.3. The additional material parameters prescribed for this ductile damage model are listed in Table 2. For the Rousselier model the initial pore volume fractions were randomly seeded in the range 4 · 103 6 f0 6 8 · 103 at the integration points to account for the inhomogeneity of actual matrix materials. All analyses discussed in this section employed periodicity boundary conditions. It is stressed that in the damage regime the responses obtained from the following simulations are not directly comparable with experimental results on SiC/Al2618-T4. The models consider only one local damage mode, ductile failure of the matrix, whereas particle fracture, which is known to play an important role in initiating damage in this composite system [30], is not accounted for. 4.1. Mesh dependence of three-dimensional models In [22] the efficacy of the present nonlocal ductile damage scheme was discussed for simple two-dimensional phase arrangements. In order to show the applicability of the approach to three-dimensional models, a periodic microgeometry based on a unit cell containing two ‘‘vertically aligned’’ spherical particles with a volume fraction of n = 0.028 was studied. Three discretizations of this unit cell were generated, which contain some 6600, 10 400 and 22 900 10-noded quadratic tetrahedral elements, respectively, compare Fig. 2. The corresponding average element sizes are h = 0.1s, h = 0.075s and h = 0.062s, where s stands for the edge length of the cube-shaped cell. The characteristic length for the nonlocal averaging procedure
Force-ratio RF2 / RF2,max
Table 2 Material parameters prescribed for the Rousselier ductile damage model employed for the Al2618-T4 matrix
1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
33
≈ 6.600 el. ≈ 10.400 el. ≈ 22.900 el.
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 0.9 Displacement-ratio U2 / U2, max
1
Fig. 3. Normalized force–displacement responses under macroscopic tensile loading obtained with the three finite element meshes for a twoparticle unit cell shown in Fig. 2 for a fixed characteristic length of 2L = 0.175s.
was set to 2L = 0.175s in all three cases, which is slightly in excess of half the particle diameter. Macroscopic uniaxial tensile loading in the 2-direction, i.e., the vertical, was applied. Fig. 3 shows the normalized macroscopic force– displacement responses predicted for the three two-particle unit cells. Effects of ductile matrix damage become noticeable in the overall behavior for U2/U2,max ’ 0.3 and rapid softening due to ductile damage is evident for U2/ U2,max ’ 0.55. Whereas there is hardly any mesh sensitivity in the former range, a limited mesh dependence is present at elevated levels of damage. Like in two-dimensional analyses [22,29], the nonlocal averaging scheme accordingly is successful in considerably reducing the mesh sensitivity in the softening regime, so that satisfactory analyses can be carried out. However, a total suppression of the mesh dependence in the predicted damage response in micromechanical analyses of inhomogeneous materials is not achieved in general. 4.2. Progress of a ductile crack in a two-dimensional phase arrangement For demonstrating the applicability of the nonlocal ductile damage models to geometrically complex
Fig. 2. Three meshes of different resolution of a two-particle unit cell (n = 0.028). Total numbers of elements are approximately 6600 (left), 10 400 (center) and 22 900 (right).
34
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
two-dimensional micromechanical problems a modified version of a phase arrangement reported by Nakamura and Suresh [31] was studied. This unit cell describes a transverse section of a composite that is reinforced by aligned continuous fibers. The cell contains 60 fibers of spherical cross section in a pseudo-random arrangement, the fiber volume fraction being n = 0.4. The model was meshed with 36 000 6-noded triangular generalized plane strain elements and, in contrast to [31], periodic boundary conditions were prescribed. The characteristic length was set to 8.2% of the fiber diameter. Fig. 4 shows the progress of a ductile crack through the matrix in this unit cell for the load case of macroscopic uniaxial transverse loading in the vertical direction. Ductile damage initiates at a number of critical
positions. The paths of the developing cracks are continued periodically over the boundaries of the cell, a behavior that is characteristic of unit cell models. The width of the cracks can be seen to encompass at least two of elements, which is a direct effect of the nonlocal averaging scheme. The normalized force–displacement diagram obtained from the above unit cell analysis is shown in Fig. 5. The crack progresses rather rapidly through the matrix, giving rise to a limited ductility in the predicted overall response. Such behavior appears to be favored by two-dimensional (as compared to three-dimensional) model geometries and by the DDIT scheme based on the Rice/Tracey damage model (as compared to the Gurson and Rousselier ductile rupture models).
Fig. 4. Predicted progress of a ductile crack through the matrix in a unit cell describing a transversally loaded SiC/Al MMC reinforced by aligned continuous fibers (n = 0.4).
1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Displacement-ratio U2 / U2, max
0.9
1
Fig. 5. Normalized force–displacement responses under macroscopic tensile loading predicted for a transversally loaded fiber reinforced SiC/ Al MMC (n 0.4).
4.3. Multi-particle unit cells At present one of the most powerful methods for simulating the thermomechanical behavior of particle reinforced composites are multi-particle unit cell models, see, e.g., [28,33,34], in which periodic homogenization is carried out on volume elements that contain a number of particles that are positioned according to the relevant arrangement statistics. Ideally, the unit cells should be sized to approach proper reference volume elements that are statistically representative of the composite to be modeled. The choice of the size of unit cells for models involving damage, however, is an unresolved problem at present. In addition it is evident that periodic homogenization must always give rise to periodically repeating patterns of damage, compare Section 4.2, and, accordingly, is not well suited to modeling the emergence of realistic macrocracks [32]. Nevertheless, unit cell models are highly useful in situations where only microcracks are present. Fig. 6 shows three multi-particle unit cells, each of which contains 10 spherical particles of equal size in a pseudorandom arrangement. Identical positions of the particle centers, generated by a random insertion algorithm, were used in the three unit cells, the particle radii being adjusted to give total particle volume fractions of n = 0.053, n = 0.111 and n = 0.2, respectively. Each of the meshes
35
contains more than 40 000 10-noded tetrahedral elements. The number of particles in the unit cells, which was chosen in view of computational requirements, is too small for a proper RVE describing a statistically homogeneous distribution of particles. As a consequence, the actual choice of the particle positions considerably influences the predicted overall responses, especially in the inelastic regime. Due to the use of identical particle centers, however, such effects are very similar for the three microgeometries considered, so that influence of the different particle volume fractions stands out clearly in the overall responses. The same characteristic length of 2L = 0.0667s was used for the three unit cells. The normalized force–displacement curves predicted for macroscopic uniaxial tensile loading of the three phase arrangements are presented in Fig. 7. The unit cell with the highest particle volume fraction can be seen to display the highest elastic stiffness, the strongest strain hardening in the elastoplastic regime, the highest maximum stress (reached at the lowest strain), and the most rapid reduction in strength due to damage. The opposite tendencies are predicted for the lowest volume fraction in the series, and the model with n = 0.111 shows an intermediate behavior. These results are in qualitative agreement with experimental trends, viz. that increases in the particle volume fraction of MMCs lead to improved stiffness and strength but to 1 0.9
Force-ratio RF2 / RF2, max
Force-ratio RF2 / RF2, max
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
0.8 0.7 0.6 0.5 0.4 0.3 0.2
20% pvf 11.1% pvf 5.3% pvf
0.1 0
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Displacement-ratio U2 / U2,max
0.9
1
Fig. 7. Force–displacement relations under macroscopic uniaxial tensile loading predicted for the three unit cells shown in Fig. 6.
Fig. 6. Undeformed meshes of three unit cells containing 10 particles each; particle positions are identical and particle radii are chosen to obtain particle volume fractions of 20% (left), 11.1% (center), and 5.3% (right).
36
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
reduced ductility. Additional simulation runs with other microgeometries are, however, needed to bolster the statistical significance of these results. Finally, it is worth noting that a series of simulations involving three-dimensional arrangements of spherical particles in a ductile matrix consistently predicted that, in the absence of other microscopic damage modes, ductile matrix damage occurs in the immediate neighborhood of the particles at positions determined by the direction of the applied load [15]. The resulting regions of local ductile damage showed little tendency to coalesce at reasonable macroscopic loads, giving rise to a rather ductile homogenized behavior, compare Fig. 7. 5. Conclusions An implementation of nonlocal ductile damage models into a commercial Finite Element program was presented that is specially geared towards use within a continuum micromechanics framework. This modeling approach strongly reduces the mesh dependence of the predicted mechanical responses of metal matrix composites under conditions where ductile damage initiates and propagates. The macroscopic and microscopic responses of multiparticle and multi-fiber unit cell models of metal matrix composites were studied, ductile damage being incorporated as the only damage mechanism. These analyses served to verify the applicability of the nonlocal ductile damage models to micromechanical simulations of considerable geometrical complexity and reproduced typical aspects of the behavior of particle reinforced MMCs. Future work is planned, on the one hand, to aim at improvements in the numerical efficiency of the nonlocal ductile damage algorithms. On the other hand, studies that combine models for different microscopic damage modes are of obvious interest. Acknowledgement The authors wish to acknowledge the financial support of the Christian Doppler Research Society. References [1] J. Llorca, A. Needleman, S. Suresh, Acta Metall. Mater. 39 (1991) 2317–2335. [2] J. Wulf, S. Schmauder, H. Fischmeister, Comput. Mater. Sci. 3 (1994) 300–306. [3] V. Tvergaard, Acta Mater. 46 (1998) 3637–3648.
[4] V. Tvergaard, A. Needleman, Int. J. Sol. Struct. 32 (1995) 1063–1077. [5] C. Hu, S. Moorthy, S. Ghosh, in: S. Ghosh et al. (Eds.), Proc. of NUMIFORM 2004, American Institute of Physics, Melville (NY), 2004, pp. 1893–1898. [6] T. Drabek, H.J. Bo¨hm, in: Z.H. Yao et al. (Eds.), Computational Mechanics—Proc. of the Sixth World Congress on Computational Mechanics, Tsinghua University Press and Springer-Verlag, Beijing, 2004, paper #262. [7] S.R. Gunawardena, S. Jansson, F.A. Leckie, in: K. Haritos et al. (Eds.), Failure Mechanisms in High Temperature Composite Materials, ASME, New York, 1991, pp. 23–30. [8] J.R. Rice, D.M. Tracey, J. Mech. Phys. Sol. 17 (1969) 201–217. [9] J.W. Hancock, A.C. Mackenzie, J. Mech. Phys. Sol. 24 (1976) 147– 169. [10] F.D. Fischer, O. Kolednik, G.X. Shan, F.G. Rammerstorfer, Int. J. Fract. 73 (1995) 345–357. [11] A.L. Gurson, J. Eng. Mater. Technol. 99 (1977) 2–15. [12] A. Needleman, V. Tvergaard, J. Mech. Phys. Sol. 35 (1987) 151–183. [13] G. Rousselier, Nucl. Eng. Design 105 (1987) 97–111. [14] S.M. Aboutayeb, Comportement a` l’endommagement des mate´riaux me´talliques he´te´roge`nes: Simulation et expe´rience, Ph.D. Thesis, Universite´ des Sciences et Technologies de Lille, Lille, France, 2000. [15] T. Drabek, Modeling of Matrix Damage in Particle Reinforced Ductile Matrix Composites, Ph.D. Thesis, Vienna University of Technology, Vienna, Austria, 2005. [16] N. Aravas, Int. J. Num. Meth. Eng. 24 (1987) 1395–1416. [17] Z.L. Zhang, Comput. Meth. Appl. Mech. Eng. 121 (1995) 29–44. [18] H. Baaser, D. Gross, in: H.A. Mang et al. (Eds.), Proc. of the Fifth World Congress on Computational Mechanics, paper #81462, Vienna, 2002. [19] G. Pijaudier-Cabot, Z.P. Bazˇant, J. Eng. Mech. 113 (1987) 1512– 1533. [20] J.B. Leblond, G. Perrin, J. Devaux, J. Appl. Mech. 61 (1994) 236– 242. [21] H. Baaser, V. Tvergaard, Comput. Meth. Appl. Mech. Eng. 192 (2003) 107–124. [22] T. Drabek, H.J. Bo¨hm, Comput. Mater. Sci. 32 (2005) 329–336. [23] M. Ga˘ra˘jeu, P. Suquet, J. Mech. Phys. Sol. 45 (1997) 873–902. [24] M. Geni, M. Kikuchi, Acta Mater. 46 (1998) 3125–3133. [25] H. Berns, A. Melander, D. Weichert, N. Asnafi, C. Broeckmann, A. Gross-Weege, Comput. Mater. Sci. 11 (1998) 166–188. [26] L.L. Mishnaevsky, M. Dong, S. Ho¨nle, S. Schmauder, Comput. Mater. Sci. 16 (1999) 133–143. [27] J. LLorca, J. Segurado, Mater. Sci. Eng. A365 (2004) 267–274. [28] J. Segurado, Micromeca´nica computacional de materiales compuestos reforzados con partı´culas, Ph.D. Thesis, Universidad Polite´cnica de Madrid, Madrid, Spain, 2004. [29] H.J. Bo¨hm, D. Duschlbauer, T. Drabek, C. Chimani, in: Z.H. Yao et al. (Eds.), Computational Mechanics—Proc. of the Sixth World Congress on Computational Mechanics, Tsinghua University Press and Springer-Verlag, Beijing, 2004, paper #258. [30] J. LLorca, C. Gonza´lez, J. Mech. Phys. Sol. 46 (1998) 1–28. [31] T. Nakamura, S. Suresh, Acta Metall. Mater. 41 (1993) 1665–1681. [32] H.J. Bo¨hm, in: H.J. Bo¨hm (Ed.), Mechanics of Microstructured Materials, CISM Courses and Lectures, Vol. 464, Springer-Verlag, Vienna, 2004, pp. 1–28. [33] A.A. Gusev, J. Mech. Phys. Sol. 45 (1997) 1449–1459. [34] H.J. Bo¨hm, W. Han, Modell. Simul. Mater. Sci. Eng. 9 (2001) 47–65.
Micromechanical finite element analysis of metal matrix composites using nonlocal ductile failure models T. Drabek 1, H.J. Bo¨hm
*
Christian Doppler Laboratory for Functionally Oriented Materials Design, Institute of Lightweight Design and Structural Biomechanics, Vienna University of Technology, A-1040 Vienna, Austria
Abstract Finite element studies of ductile damage in the matrix of particle and fiber reinforced metal matrix composites are presented. Three material models capable of supporting such simulations at the constituent level are discussed within a unified framework. They comprise an element removal technique triggered by a ductile damage indicator as well as versions of the ductile rupture models of Gurson and Rousselier. Nonlocal averaging is employed for reducing the mesh dependence typically displayed by continuum damage methods upon the onset of local softening. Implementation issues specific to the use of nonlocal damage models in a continuum micromechanics framework are discussed. The efficacy of the approach in limiting the mesh sensitivity of the predicted behavior of composites subject to matrix damage is demonstrated for a simple three-dimensional matrix-particle configuration. Applications of the method to multi-particle and multi-fiber unit cells subjected to uniaxial tensile loads are presented and effects of microgeometrical parameters on the mechanical response are discussed. 2005 Elsevier B.V. All rights reserved. PACS: 46.30.Nz; 83.70.Dk; 83.10.Ff; 02.70.Dh; 11.10.Lm Keywords: Ductile damage models; Nonlocal averaging; Metal matrix composites; Micromechanics
1. Introduction In metal matrix composites (MMCs) and related materials damage and failure are caused by debonding of the matrix-reinforcement interface, by brittle fracture of the reinforcing fibers or particles, and/or by ductile failure of the matrix. The present work concentrates on the development and application of algorithms that are capable of modeling the latter failure mode. Other damage mechanisms are not accounted for.
* Corresponding author. Tel.: +43 1 58801 31712; fax: +43 1 58801 31799. E-mail address: [email protected] (H.J. Bo¨hm). URL: http://www.ilsb.tuwien.ac.at (H.J. Bo¨hm). 1 Present address: Voith–Siemens Power Generation, A-3100 St. Po¨lten, Austria.
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.12.032
Ductile damage and failure of metals are caused by the growth and coalescence of existing voids and by the nucleation of new voids by debonding from tiny inclusions such as precipitates. A number of ductile damage models have been proposed in the literature for describing the behavior of metals undergoing this damage process. The implementation of such models into a finite element code supports the investigation of the failure of ductile constituents of inhomogeneous materials under given macroscopic loading conditions at length scales where the phases are geometrically resolved, i.e., on micro- and meso-scales. Such modeling strategies are referred to as continuum micromechanical methods. They allow the overall behavior of composites undergoing ductile damage and failure of the matrix to be obtained by homogenization. Standard ductile damage models have been combined with unit cells for continuum micromechanical studies of MMCs by a number of
30
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
authors, the majority of whom used axisymmetric cells, simple three-dimensional models or planar multi-particle geometries, see, e.g., [1–3]. It is well known that in the softening regime ductile damage models (like continuum damage models in general) tend to give rise to a marked mesh dependence of the predicted results. This behavior is caused by the loss of ellipticity of the governing partial differential equations once softening sets in. In order to overcome this mesh sensitivity the present study applies a nonlocal averaging approach to three different ductile damage models. This procedure introduces an additional material parameter in the form of a characteristic length. An important precondition for the use of the above modeling strategy, viz., that the size of individual voids is much smaller than the width of the ‘‘matrix bridges’’ between neighboring particles or fibers, is usually fulfilled for micromechanical studies of MMCs and related materials. Ductile damage models that employ nonlocal averaging were introduced into continuum micromechanics by Tvergaard and Needleman [4] and have recently seen further development [5,6]. In the following, nonlocal ductile damage models are described that are specifically geared towards the requirements posed by micromechanical studies of the local and overall response of metal matrix composites that undergo ductile damage. 2. Implemented ductile damage models A number of ductile damage models can be found in the literature that are suitable for describing ductile failure processes via appropriate constitutive equations that account for void growth. In these models the yield function is typically treated in extension of J2-plasticity theory, so that expressions of the type U¼
req rf ¼ 0 D
or
b ¼ req D ¼ 0 U rf
ð1Þ
may be used to describe the yield condition. Here req stands for the von Mises equivalent stress, rf is the flow stress, and D is a parameter that indicates the progress of damage and varies between 1 and 0. A completely damage-free material (i.e., one in which no voids are present) is indicated by D ¼ 1 and follows standard J2-plasticity theory. When D reaches zero, ductile failure occurs and the von Mises stress req vanishes at the same time, so that the term req =D in Eq. (1) remains finite. Note that the above definitions differ from the ones commonly used for damage parameters in continuum damage mechanics. Formally, the ductile damage models discussed in this section differ only in the definition of D, with the mean stress rm (and, consequently, the stress triaxiality g = rm/ req) playing an important role in determining D in all three models. For this reason these models can be classified as pressure dependent plasticity models, in which the first invariant of the stress tensor must be considered in order to obtain a closed yield surface in stress space.
2.1. Ductile damage indicator triggered (DDIT) models Gunawardena et al. [7] proposed a ductile damage indicator that is capable of detecting local ductile failure of a given material point by postprocessing operations following a standard elastoplastic finite element analysis. This method, which is based on work of Rice and Tracey [8] as well as Hancock and Mackenzie [9], was discussed in depth by Fischer et al. [10]. The increment of the ductile damage indicator D can be written as exp 32 g Deeq;p ; ð2Þ DD ¼ 1:65e0 where eeq,p stands for the accumulated equivalent plastic strain and the material parameter e0 is a reference failure strain that can be calibrated by appropriate tensile tests [10]. Regions where local ductile failure has occurred are indicated by the accumulated value of D reaching or exceeding a critical value Dc, which is usually set to one. The ductile damage indicator, however, does not uniquely correlate with the degradation of the mechanical properties due to evolving damage in the range 0 < D < Dc and per se does not introduce any loss of stiffness or strength into the material model. In order to achieve a local softening response in a finite element context the ductile damage indicator can be used to trigger either the removal of the element or a marked reduction of load carrying capacity and stiffness (such that req ! 0 and rm ! 0) at the integration point(s) involved. For D < Dc the local material behavior continues to follow J2 flow theory in such a model. Accordingly, D as defined in Eq. (1) takes the form 1 for D < Dc ; D¼ ð3Þ 0 for D P Dc in descriptions of this type, which are referred to as ductile damage indicator triggered (DDIT) models in the following. 2.2. Ductile rupture models Another group of approaches to describing void growth and coalescence in ductile materials are based on continuous modifications of the von Mises flow potential in dependence on the mean stress, rm, and the volume fraction of voids, f. The best known continuum model of the above type was developed by Gurson [11] and extended by Needleman and Tvergaard [12] to include effects of void nucleation and void coalescence. For the resulting GTN model the parameter D according to Eq. (1) can be denoted as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3q2 rm ð4Þ D ¼ 2q1 f cosh þ 1 þ q3 f 2 . 2rf where q1, q2 as well as q3 are adjustable parameters. The void volume fraction at which ductile failure takes place
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
can be evaluated by rearranging Eq. (4) at D ¼ req ¼ rm ¼ 0 to give ff ¼
1 q1
for q3 ¼ q21 .
ð5Þ
An alternative ductile rupture model was proposed in the form of a thermodynamically based continuum damage mechanics method by Rousselier [13] and extended by Aboutayeb [14], with D taking the form rm e D ¼ 1 q4 f D exp ; ð6Þ q4 rf e are material parameters. An imporwhere q4 as well as D tant feature of the extended Rousselier model [14] lies in the fact that it allows to calculate the void volume fraction at failure as ff ¼
1
ð7Þ
e q4 D
in dependence of the material parameters only. The increment of the void volume fraction can be calculated in both the Gurson and Rousselier models by evolution equations of the type Df ¼ ð1 f ÞBDep ;
ð8Þ
where Dep stands for the increment of the first invariant of the tensor of plastic strains and B is defined as 1 for f 6 fc ; B¼ ð9Þ > 1 for f > fc . Eqs. (8) and (9) allow the void volume fraction to be integrated up starting from initial values f0. To avoid a discontinuity in the history of Df in Eq. (8), B can be chosen to be a linear function of the void volume fraction, Bðf Þ ¼ 1 þ
Bf 1 ðf fc Þ; ff fc
ð10Þ
for f > fc. Here Bf = B(f = ff) represents the value of B at total failure of the material point. Additional details of the models are given in [6,15]. The above three models for describing ductile damage were implemented into the general purpose Finite Element code ABAQUS/Standard V6 (Abaqus Inc., Pawtucket, RI, 2001) as user defined material subroutines (UMATs). For the DDIT model an implementation via a user defined field subroutine (USDFLD) was developed in addition. In the UMATs the consistent tangent (Jacobian) required for evaluating the stiffness matrix of the model, Ccons ¼
oDr or ¼ oDe oe
ð11Þ
was derived by the procedures proposed by Aravas [16] and extended by Zhang [17] for pressure dependent plasticity models. Predictions obtained by the basic (‘‘local’’) versions of all of the above models typically show a clear mesh depen-
31
dence once the acoustic tensor defined by the double contraction Qij ¼ C ijkl nk nl
ð12Þ
of the material tensor C and the normal vector n, which indicates possible wave propagation directions, becomes negative definite, see, e.g., the discussion by Baaser and Gross [18]. This mesh sensitivity, on the one hand, takes the form of a marked dependence of the crack path on the employed finite element mesh, the assemblage of failed elements making up the crack typically having a width of one element. On the other hand, the resulting overall force–displacement responses tend to vary in such a way that the energy required for failure decreases with increasing mesh refinement. As a consequence, the advantageous convergence properties of the Finite Element method upon mesh refinement are no longer guaranteed. These undesirable effects tend to be less marked in DDIT schemes than in ductile rupture models. 3. Nonlocal regularization technique The mesh sensitivity of ductile damage algorithms can be markedly reduced by applying regularization techniques. In such approaches an additional material parameter is introduced in the form of a characteristic length 2L. Physically this internal length scale of the material may be interpreted as the size of a process zone or as a function of the distance between individual voids or the distance between dimples in the fracture surfaces [4], but is not fully understood at present. The regularization method employed in the present work is based on smoothing the rate of an appropriate damage variable q by nonlocal averaging [19], a delocalization function proposed by Leblond [20] being employed. This nonlocal averaging procedure can be denoted in terms of the increment of the nonlocal damage variable, DqNL, as 2 32 DqNL ðxi Þ ¼
n 6 1 X 6 DqL ðb x j Þ6 4 W ðxi Þ j¼1
1þ
7 7 8 7 D Vb j ; 5 jxi b xjj 1
ð13Þ
L
the normalizing factor W being defined as 2 32 W ðxi Þ ¼
n 6 X 6 6 4 j¼1
7 7 8 7 D Vb j . 5 jxi b xjj
1þ
1
ð14Þ
L
Here qL is the local damage variable and qNL is its nonlocal counterpart, jxi b x j j stands for the distance between two integration points i and j, D Vb j is the volume associated with integration point j, and the sums run over all integration points j within a distance L from a given integration point i. The variables subjected to smoothing are the
32
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
damage indicator D in the case of the DDIT scheme and the void volume fraction f in the case of the GTN and Rousselier models. The nonlocal smoothing routines employed in the present work employ lists of neighbors for each integration point to allow efficient averaging by Eqs. (13) and (14). These lists are generated in a preprocessor step and can be updated in the course of an analysis to account for the effects of large strains. Nonlocal averaging is applied only within one constituent, viz., the metallic matrix. This feature is physically reasonable for composites consisting of a ductile and an elastic phase, but may have to be amended in cases where more than one constituent can be expected to be subject to ductile damage. Integration points where failure has occurred are excluded from the sums in Eqs. (13) and (14). Provision is also made for consistently handling nonlocal averaging in the presence of the periodicity and/or symmetry boundary conditions that are typically encountered in problems involving periodic homogenization. Fig. 1 schematically depicts the region of influence of an integration point (marked by an ‘‘x’’) that is located near the border of a two-dimensional unit cell employing periodicity boundary conditions. The user subroutine interfaces provided by ABAQUS/ Standard allow nonlocal averaging to be carried out only at the ends of increments, where the UEXTERNALDB routine can be invoked, but not at the level of the equilibrium iterations. Comparisons with published results [21] showed that the resulting modification to the ‘‘standard’’ smoothing algorithms gives rise only to minor changes in the predicted crack paths and overall stress–strain
responses [22]. Details of the averaging procedure are discussed in [6,15]. It is worth noting that the use of local damage models may be interpreted as being equivalent to implicitly specifying a distribution of characteristic lengths via the element sizes. For simple microgeometries, a judicious choice of the discretization can accordingly give rise to valid results in the softening regime. In the case of geometrically complex multi-fiber or multi-particle models, however, where the phase geometry enforces rather irregular meshes with local element sizes that are driven by the requirements of resolving the stress and strain fields, local damage models may introduce considerable inconsistencies into the material descriptions. Finally, it must be emphasized that regularization by nonlocal averaging is viable only if the mesh size in regions where ductile damage may occur is smaller than (or at most equal to) the characteristic length 2L. This condition is usually met by micromechanical studies, but may be very difficult to fulfill for macroscopic analyses of components or structures. 4. Application to micromechanical studies of MMCs A number of applications of local ductile damage models to micromechanical studies of ductile matrix composites have been reported in the literature, where the Gurson ductile rupture model, see, e.g., [1,3,23,24], and DDIT-like schemes, see, e.g., [2,25,26], were employed. Recent work of this type has also involved multi-particle unit cells [27,28]. The use of ductile damage models with nonlocal averaging in the context of continuum micromechanics was initiated by Tvergaard and Needleman [4]. Recently such approaches have been further developed for use with more complex microgeometries [5,6,22,29]. In the following the nonlocal damage routines presented in Sections 2 and 3 are employed for studying the effects of ductile damage on the mechanical behavior of particle and fiber reinforced MMCs. The material properties of the matrix correspond to Al2618-T4 and those of the fibrous or particulate reinforcements to SiC. The corresponding material parameters are listed in Table 1, where E stands for the Young’s modulus and m for the Poisson number. J2 plasticity with a modified Ludwik hardening law was used to describe the damage-free behavior of the matrix, with ry0 standing for the initial yield stress, m for the hardening modulus, and n for the hardening exponent. Results obtained with two-dimensional and threedimensional unit cell models are presented, which pertain
Table 1 Material parameters used for the elastoplastic Al2618-T4 matrix (modified Ludwik hardening law) and the elastic SiC reinforcements Fig. 1. Sketch of regions where Eqs. (13) and (14) are applied for nonlocal averaging in a finite element mesh employing periodic boundary conditions.
Al2618-T4 matrix SiC reinforcement
E (GPa)
m [1]
ry0 (MPa)
m (MPa)
n [1]
70 450
0.30 0.17
184 –
722.7 –
0.49 –
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
q4 [1]
e [1] D
f0 [1]
fc [1]
Bf [1]
1.0
2.0
0.004–0.008 (random)
0.2
3.0
to the uniaxial mechanical response of composites reinforced by continuous aligned fibers and by equiaxed particles, respectively. The investigation of a fiber reinforced MMC, Section 4.2, was based on the nonlocal DDIT scheme, the reference strain being chosen as e0 = 0.2. A nonlocal extended Rousselier model was used in the studies involving reinforcement by particles, Sections 4.1 and 4.3. The additional material parameters prescribed for this ductile damage model are listed in Table 2. For the Rousselier model the initial pore volume fractions were randomly seeded in the range 4 · 103 6 f0 6 8 · 103 at the integration points to account for the inhomogeneity of actual matrix materials. All analyses discussed in this section employed periodicity boundary conditions. It is stressed that in the damage regime the responses obtained from the following simulations are not directly comparable with experimental results on SiC/Al2618-T4. The models consider only one local damage mode, ductile failure of the matrix, whereas particle fracture, which is known to play an important role in initiating damage in this composite system [30], is not accounted for. 4.1. Mesh dependence of three-dimensional models In [22] the efficacy of the present nonlocal ductile damage scheme was discussed for simple two-dimensional phase arrangements. In order to show the applicability of the approach to three-dimensional models, a periodic microgeometry based on a unit cell containing two ‘‘vertically aligned’’ spherical particles with a volume fraction of n = 0.028 was studied. Three discretizations of this unit cell were generated, which contain some 6600, 10 400 and 22 900 10-noded quadratic tetrahedral elements, respectively, compare Fig. 2. The corresponding average element sizes are h = 0.1s, h = 0.075s and h = 0.062s, where s stands for the edge length of the cube-shaped cell. The characteristic length for the nonlocal averaging procedure
Force-ratio RF2 / RF2,max
Table 2 Material parameters prescribed for the Rousselier ductile damage model employed for the Al2618-T4 matrix
1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
33
≈ 6.600 el. ≈ 10.400 el. ≈ 22.900 el.
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 0.9 Displacement-ratio U2 / U2, max
1
Fig. 3. Normalized force–displacement responses under macroscopic tensile loading obtained with the three finite element meshes for a twoparticle unit cell shown in Fig. 2 for a fixed characteristic length of 2L = 0.175s.
was set to 2L = 0.175s in all three cases, which is slightly in excess of half the particle diameter. Macroscopic uniaxial tensile loading in the 2-direction, i.e., the vertical, was applied. Fig. 3 shows the normalized macroscopic force– displacement responses predicted for the three two-particle unit cells. Effects of ductile matrix damage become noticeable in the overall behavior for U2/U2,max ’ 0.3 and rapid softening due to ductile damage is evident for U2/ U2,max ’ 0.55. Whereas there is hardly any mesh sensitivity in the former range, a limited mesh dependence is present at elevated levels of damage. Like in two-dimensional analyses [22,29], the nonlocal averaging scheme accordingly is successful in considerably reducing the mesh sensitivity in the softening regime, so that satisfactory analyses can be carried out. However, a total suppression of the mesh dependence in the predicted damage response in micromechanical analyses of inhomogeneous materials is not achieved in general. 4.2. Progress of a ductile crack in a two-dimensional phase arrangement For demonstrating the applicability of the nonlocal ductile damage models to geometrically complex
Fig. 2. Three meshes of different resolution of a two-particle unit cell (n = 0.028). Total numbers of elements are approximately 6600 (left), 10 400 (center) and 22 900 (right).
34
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
two-dimensional micromechanical problems a modified version of a phase arrangement reported by Nakamura and Suresh [31] was studied. This unit cell describes a transverse section of a composite that is reinforced by aligned continuous fibers. The cell contains 60 fibers of spherical cross section in a pseudo-random arrangement, the fiber volume fraction being n = 0.4. The model was meshed with 36 000 6-noded triangular generalized plane strain elements and, in contrast to [31], periodic boundary conditions were prescribed. The characteristic length was set to 8.2% of the fiber diameter. Fig. 4 shows the progress of a ductile crack through the matrix in this unit cell for the load case of macroscopic uniaxial transverse loading in the vertical direction. Ductile damage initiates at a number of critical
positions. The paths of the developing cracks are continued periodically over the boundaries of the cell, a behavior that is characteristic of unit cell models. The width of the cracks can be seen to encompass at least two of elements, which is a direct effect of the nonlocal averaging scheme. The normalized force–displacement diagram obtained from the above unit cell analysis is shown in Fig. 5. The crack progresses rather rapidly through the matrix, giving rise to a limited ductility in the predicted overall response. Such behavior appears to be favored by two-dimensional (as compared to three-dimensional) model geometries and by the DDIT scheme based on the Rice/Tracey damage model (as compared to the Gurson and Rousselier ductile rupture models).
Fig. 4. Predicted progress of a ductile crack through the matrix in a unit cell describing a transversally loaded SiC/Al MMC reinforced by aligned continuous fibers (n = 0.4).
1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Displacement-ratio U2 / U2, max
0.9
1
Fig. 5. Normalized force–displacement responses under macroscopic tensile loading predicted for a transversally loaded fiber reinforced SiC/ Al MMC (n 0.4).
4.3. Multi-particle unit cells At present one of the most powerful methods for simulating the thermomechanical behavior of particle reinforced composites are multi-particle unit cell models, see, e.g., [28,33,34], in which periodic homogenization is carried out on volume elements that contain a number of particles that are positioned according to the relevant arrangement statistics. Ideally, the unit cells should be sized to approach proper reference volume elements that are statistically representative of the composite to be modeled. The choice of the size of unit cells for models involving damage, however, is an unresolved problem at present. In addition it is evident that periodic homogenization must always give rise to periodically repeating patterns of damage, compare Section 4.2, and, accordingly, is not well suited to modeling the emergence of realistic macrocracks [32]. Nevertheless, unit cell models are highly useful in situations where only microcracks are present. Fig. 6 shows three multi-particle unit cells, each of which contains 10 spherical particles of equal size in a pseudorandom arrangement. Identical positions of the particle centers, generated by a random insertion algorithm, were used in the three unit cells, the particle radii being adjusted to give total particle volume fractions of n = 0.053, n = 0.111 and n = 0.2, respectively. Each of the meshes
35
contains more than 40 000 10-noded tetrahedral elements. The number of particles in the unit cells, which was chosen in view of computational requirements, is too small for a proper RVE describing a statistically homogeneous distribution of particles. As a consequence, the actual choice of the particle positions considerably influences the predicted overall responses, especially in the inelastic regime. Due to the use of identical particle centers, however, such effects are very similar for the three microgeometries considered, so that influence of the different particle volume fractions stands out clearly in the overall responses. The same characteristic length of 2L = 0.0667s was used for the three unit cells. The normalized force–displacement curves predicted for macroscopic uniaxial tensile loading of the three phase arrangements are presented in Fig. 7. The unit cell with the highest particle volume fraction can be seen to display the highest elastic stiffness, the strongest strain hardening in the elastoplastic regime, the highest maximum stress (reached at the lowest strain), and the most rapid reduction in strength due to damage. The opposite tendencies are predicted for the lowest volume fraction in the series, and the model with n = 0.111 shows an intermediate behavior. These results are in qualitative agreement with experimental trends, viz. that increases in the particle volume fraction of MMCs lead to improved stiffness and strength but to 1 0.9
Force-ratio RF2 / RF2, max
Force-ratio RF2 / RF2, max
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
0.8 0.7 0.6 0.5 0.4 0.3 0.2
20% pvf 11.1% pvf 5.3% pvf
0.1 0
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Displacement-ratio U2 / U2,max
0.9
1
Fig. 7. Force–displacement relations under macroscopic uniaxial tensile loading predicted for the three unit cells shown in Fig. 6.
Fig. 6. Undeformed meshes of three unit cells containing 10 particles each; particle positions are identical and particle radii are chosen to obtain particle volume fractions of 20% (left), 11.1% (center), and 5.3% (right).
36
T. Drabek, H.J. Bo¨hm / Computational Materials Science 37 (2006) 29–36
reduced ductility. Additional simulation runs with other microgeometries are, however, needed to bolster the statistical significance of these results. Finally, it is worth noting that a series of simulations involving three-dimensional arrangements of spherical particles in a ductile matrix consistently predicted that, in the absence of other microscopic damage modes, ductile matrix damage occurs in the immediate neighborhood of the particles at positions determined by the direction of the applied load [15]. The resulting regions of local ductile damage showed little tendency to coalesce at reasonable macroscopic loads, giving rise to a rather ductile homogenized behavior, compare Fig. 7. 5. Conclusions An implementation of nonlocal ductile damage models into a commercial Finite Element program was presented that is specially geared towards use within a continuum micromechanics framework. This modeling approach strongly reduces the mesh dependence of the predicted mechanical responses of metal matrix composites under conditions where ductile damage initiates and propagates. The macroscopic and microscopic responses of multiparticle and multi-fiber unit cell models of metal matrix composites were studied, ductile damage being incorporated as the only damage mechanism. These analyses served to verify the applicability of the nonlocal ductile damage models to micromechanical simulations of considerable geometrical complexity and reproduced typical aspects of the behavior of particle reinforced MMCs. Future work is planned, on the one hand, to aim at improvements in the numerical efficiency of the nonlocal ductile damage algorithms. On the other hand, studies that combine models for different microscopic damage modes are of obvious interest. Acknowledgement The authors wish to acknowledge the financial support of the Christian Doppler Research Society. References [1] J. Llorca, A. Needleman, S. Suresh, Acta Metall. Mater. 39 (1991) 2317–2335. [2] J. Wulf, S. Schmauder, H. Fischmeister, Comput. Mater. Sci. 3 (1994) 300–306. [3] V. Tvergaard, Acta Mater. 46 (1998) 3637–3648.
[4] V. Tvergaard, A. Needleman, Int. J. Sol. Struct. 32 (1995) 1063–1077. [5] C. Hu, S. Moorthy, S. Ghosh, in: S. Ghosh et al. (Eds.), Proc. of NUMIFORM 2004, American Institute of Physics, Melville (NY), 2004, pp. 1893–1898. [6] T. Drabek, H.J. Bo¨hm, in: Z.H. Yao et al. (Eds.), Computational Mechanics—Proc. of the Sixth World Congress on Computational Mechanics, Tsinghua University Press and Springer-Verlag, Beijing, 2004, paper #262. [7] S.R. Gunawardena, S. Jansson, F.A. Leckie, in: K. Haritos et al. (Eds.), Failure Mechanisms in High Temperature Composite Materials, ASME, New York, 1991, pp. 23–30. [8] J.R. Rice, D.M. Tracey, J. Mech. Phys. Sol. 17 (1969) 201–217. [9] J.W. Hancock, A.C. Mackenzie, J. Mech. Phys. Sol. 24 (1976) 147– 169. [10] F.D. Fischer, O. Kolednik, G.X. Shan, F.G. Rammerstorfer, Int. J. Fract. 73 (1995) 345–357. [11] A.L. Gurson, J. Eng. Mater. Technol. 99 (1977) 2–15. [12] A. Needleman, V. Tvergaard, J. Mech. Phys. Sol. 35 (1987) 151–183. [13] G. Rousselier, Nucl. Eng. Design 105 (1987) 97–111. [14] S.M. Aboutayeb, Comportement a` l’endommagement des mate´riaux me´talliques he´te´roge`nes: Simulation et expe´rience, Ph.D. Thesis, Universite´ des Sciences et Technologies de Lille, Lille, France, 2000. [15] T. Drabek, Modeling of Matrix Damage in Particle Reinforced Ductile Matrix Composites, Ph.D. Thesis, Vienna University of Technology, Vienna, Austria, 2005. [16] N. Aravas, Int. J. Num. Meth. Eng. 24 (1987) 1395–1416. [17] Z.L. Zhang, Comput. Meth. Appl. Mech. Eng. 121 (1995) 29–44. [18] H. Baaser, D. Gross, in: H.A. Mang et al. (Eds.), Proc. of the Fifth World Congress on Computational Mechanics, paper #81462, Vienna, 2002. [19] G. Pijaudier-Cabot, Z.P. Bazˇant, J. Eng. Mech. 113 (1987) 1512– 1533. [20] J.B. Leblond, G. Perrin, J. Devaux, J. Appl. Mech. 61 (1994) 236– 242. [21] H. Baaser, V. Tvergaard, Comput. Meth. Appl. Mech. Eng. 192 (2003) 107–124. [22] T. Drabek, H.J. Bo¨hm, Comput. Mater. Sci. 32 (2005) 329–336. [23] M. Ga˘ra˘jeu, P. Suquet, J. Mech. Phys. Sol. 45 (1997) 873–902. [24] M. Geni, M. Kikuchi, Acta Mater. 46 (1998) 3125–3133. [25] H. Berns, A. Melander, D. Weichert, N. Asnafi, C. Broeckmann, A. Gross-Weege, Comput. Mater. Sci. 11 (1998) 166–188. [26] L.L. Mishnaevsky, M. Dong, S. Ho¨nle, S. Schmauder, Comput. Mater. Sci. 16 (1999) 133–143. [27] J. LLorca, J. Segurado, Mater. Sci. Eng. A365 (2004) 267–274. [28] J. Segurado, Micromeca´nica computacional de materiales compuestos reforzados con partı´culas, Ph.D. Thesis, Universidad Polite´cnica de Madrid, Madrid, Spain, 2004. [29] H.J. Bo¨hm, D. Duschlbauer, T. Drabek, C. Chimani, in: Z.H. Yao et al. (Eds.), Computational Mechanics—Proc. of the Sixth World Congress on Computational Mechanics, Tsinghua University Press and Springer-Verlag, Beijing, 2004, paper #258. [30] J. LLorca, C. Gonza´lez, J. Mech. Phys. Sol. 46 (1998) 1–28. [31] T. Nakamura, S. Suresh, Acta Metall. Mater. 41 (1993) 1665–1681. [32] H.J. Bo¨hm, in: H.J. Bo¨hm (Ed.), Mechanics of Microstructured Materials, CISM Courses and Lectures, Vol. 464, Springer-Verlag, Vienna, 2004, pp. 1–28. [33] A.A. Gusev, J. Mech. Phys. Sol. 45 (1997) 1449–1459. [34] H.J. Bo¨hm, W. Han, Modell. Simul. Mater. Sci. Eng. 9 (2001) 47–65.