Linear and elastic–plastic fracture mechanics revisited by use of Fourier transforms – theory and application

Computational Materials Science 16 (1999) 186±196

Linear and elastic±plastic fracture mechanics revisited by use of Fourier transforms ± theory and application K.P. Herrmann a,*, W.H. M uller b, S. Neumann a a

Laboratorium f ur Technische Mechanik, Universit at Paderborn, Pohlweg 47-49, 33098 Paderborn, Germany Department of Mechanical and Chemical Engineering, Heriot-Watt University, Edinburg EH144AS, UK

b

Abstract This paper investigates whether and how discrete Fourier transforms (DFT) can be used to compute the local stress/ strain distribution around holes in externally loaded two-dimensional representative volume elements (RVEs). To this end, the properties of DFT are ®rst reviewed and then applied to the solution of linear elastic and time-dependent elastic plastic material response. The equivalent inclusion method is used to derive a functional equation which allows for the numerical computation of stresses and strains within an RVE with heterogeneities of arbitrary shape and sti€ness. This functional equation is then specialized to the case of circular and elliptical holes of di€erent minor axes which eventually degenerate into Grith cracks. The numerically predicted stresses and strains are compared to the corresponding analytical solutions for a single circular as well as an elliptical hole in an in®nitely large plate under tension as well as to ®nite element calculations (for time-independent elastic/plastic material response). Ó 1999 Elsevier Science B.V. All rights reserved. Keywords: Discrete Fourier transforms; Composite materials; Stress concentration; Stress intensity factor

1. Introduction Inclusions and precipitates may lead to stress concentrations and to the formation of residual stresses in solid and, as a consequence, in¯uence the development of its micro-morphology. An example from brittle material is partially stabilized Zirconia (PSZ). Zirconia is stabilized at room temperature in its cubic form by addition of oxides. However, due to ¯uctuations of the stabilizer the cubic material may locally transform and be-

* Corresponding author. Tel.: +49-05251-6022832284; fax: +49-05251-603483. E-mail address: [email protected] (K.P. Herrmann).

come tetragonal. This transformation is associated with a change in shape and volume and, consequently, local stresses will arise which, in turn, in¯uence the di€usion and concentration level of the stabilizer. Eventually a point of thermo-mechanical equilibrium is reached (see Fig. 1, top left, [1]). Another example for stress-driven micromorphology changes are Ni-base superalloys where c0 -cubes of the material precipitate within a c-matrix. The lattice parameters as well as the state of order of both phases are di€erent. Due to coherency, high internal stresses arise near the interface. Over time these will lead to a growth of the cubes. In addition, externally superimposed mechanical loads lead to rafting, i.e., a directional coarsening of the c0 -phase (see Fig. 1, top right, [2]). As in the case of the PSZ, the transition from

0927-0256/99/$ - see front matter Ó 1999 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 9 9 ) 0 0 0 6 1 - 0

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

187

Fig. 1. Micro-morphology changes in PSZ-ceramics, Ni-base superalloys, tin-lead solder, and pores in a strained metal.

c to c0 results from coupling of di€usion, atomic state of order, and from the local stresses. In a similar way, the phase separation and coarsening observed in tin-lead based binary solders can be accelerated considerably by change of temperature and by application of external mechanical stresses (see Fig. 1, bottom left, [3]). Finally, micro-morphology changes under stress are also observed during the last phase of the straining of metals when pores are formed and grow (Fig. 1, bottom right, [4]). The local concentration ®eld, c, the state of order, S, and the local stresses and strains, rij and eij , are related to each other through equations of the Cahn±Hilliard±Allen type (t ± time, xi ± position): oc oJi ‡ ˆ 0; ot oxi

oS ˆ PS : ot

…1†

The di€usion ¯ux, Ji , is explicitly given by [5]   
o ow o2 c o 1ÿ ekl ÿ ekl rkl ÿ Ji  ÿ ÿ akl oxi oc oxk oxl oc 2 …2†

and for the production of the order parameter, PS [2],   
ow o2 S o 1ÿ  ÿ ÿ bkl ekl ÿ ekl rkl : PS  ÿ oS oxk oxl oS 2 …3† In these relations, w stands for the local part of the free energy density excluding mechanical contributions, eij are the eigenstrains (due to thermal stresses and/or phase transitions), akl is related to the matrix of surface tensions, and bkl is a coecient matrix complementary to akl , but for the order parameter, S. The precise knowledge of the local stress/strain distributions is a necessary prerequisite for a simulation of the evolving microstructures shown in Fig. 1. Note that the local stresses will change while the micro-structure evolves. This necessitates use and development of algorithms which quickly react to changes in geometry. In order to solve equations of such complexity, discrete Fourier transforms (DFTs) have recently been used in a

188

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

very promising manner [2,5,6]. They allow to obtain, ®rst, the local stresses and strains in heterogeneous (non-linear) materials which are, second, used to predict the micro-morphology, i.e., for the solution of Eqs. (1)±(3) [2,5]. In this paper, the focus is on the prediction of local strains and stresses, and on an investigation as to whether and how DFT can also be used for their accurate computation under the following constraints: · large di€erences in elastic constants between the heterogeneity and the matrix (in particular for the extreme case of a hole in a matrix); · a stress-concentrating geometry: elliptic holes degenerate into slit-like cracks; · non-linear material behavior: a time-independent elastic/plastic matrix interspersed with holes. We proceed as follows: · Section 2 presents a summary of the equations used for linear as well as non-linear stress/strain analysis and how DFT can be used for their solution. · Section 3 is dedicated to the presentation of ®rst results, such as a comparison of analytically and FE-based stress-strain data to data obtained by application of DFT to KirschÕs problem, as well

as the calculation of stress intensity factors (SIFs).

2. Discrete Fourier transforms and solid-state problems 2.1. Fourier's theorem and di€erentiation rules in discrete form Following the procedures outlined in [7] consider a periodically arranged array of representative volume elements (RVEs), as indicated in Fig. 2, which is discretized by N
x  Ny  Nz points ÿ at positions x ˆ ax hx ; ay hy ; az hz in physical space. De®ne discrete ®eld variables, f …a†, in each point, a. Their DFT, f^…s†, follows from f^…s† ˆ Y … f …a†† ˆ ÿ

1

N x ÿ1 X

N y ÿ1 X

N z ÿ1 X

f … a†
1=2 Nx Ny Nz az ˆ0 ax ˆ0 ay ˆ0    ax s x ay s y az s z ‡ : …4† ‡  exp i 2p Nx Ny Nz

Inversely, if f^…s† is known, and due to the periodicity conditions which are assumed to hold

Fig. 2. 2D-Heterogeneity Problem: RVE before and after discretization.

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

across the RVE, the corresponding quantities in physical space, f …a†, are given by   f …a† ˆ Y ÿ1 f^…s† 1

y ÿ1 N N x ÿ1 N z ÿ1 X X X

f^…s† Nx Ny Nz sx ˆ0 sy ˆ0 sz ˆ0    ax s x ay s y az s z  exp ÿ i2p ‡ : ‡ Nx Ny Nz

ˆ ÿ


1=2

…5†

It should be pointed out that these are ®nite sums which can be determined exactly, e.g., by FFT, and that this is not an approximation of the continuous Fourier theorem. Di€erentiation rules hold:     ÿ
of ˆ ÿnj sj Y ‰ f Š‡O h2j ; j 2 f x;y;zg; Y oxj   …6† ÿ
i sj ; nj sj ˆ ÿ sin 2p Nj hj  Y

 ÿ
o2 f ˆ nij Y ‰ f Š ‡ O hi hj ; oxi oxj     ÿ
2 sj ÿ1 ; njj sj ˆ 2 cos 2p Nj hj     ÿ
1 si sj nij sj ˆ ‡ cos 2p Ni Nj 2hi hj    si sj ÿ cos 2p ÿ ; Ni Nj

…7†

which follow by approximations of the di€erentials by central di€erence quotients (cf. [7]).

189

2.2. Solution of continuum mechanics problems by DFT Starting point are the equations for static equilibrium of forces, HookeÕs law, and linearized kinematic conditions (ui denotes the displacements): ÿ
orij ˆ 0; rij ˆ Cijkl ekl ÿ ekl ; oxj   1 ouk oul ‡ : ekl ˆ 2 oxl oxk

…8†

Moreover, Cijkl is the sti€ness matrix which, due to the heterogeneities, will depend on space:   ‡ ‡ ÿ ÿ h…x† Cijkl ÿ Cijkl ; Cijkl …x† ˆ Cijkl  …9† 1 if x 2 ÿ; h…x† ˆ 0 if x 2 ‡; where ‡ refers to the matrix and ÿ to the heterogeneities. An average ¯ow of stress is assumed to cross the periphery of the RVE: Z

Z Ly

rxx …x† dy ˆ Ly r0x ;

Lx

ryy …x† dx ˆ Lx r0y :

…10†

We now follow Eshelby [8] and apply the equivalent inclusion method to rewrite HookeÕs law as follows:

Fig. 3. KirschÕs problem for a circular hole and the normal stress in horizontal direction.

190

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196



 rij …x† ˆ Cijkl …x† ekl …x† ÿ ekl …x†   ! H ˆ Cijkl ekl …x† ÿ ekl …x† ÿ eH kl …x† ;

H ‡ ˆ Cijkl ; Cijkl

…11†

H denotes an auxiliary sti€ness which is where Cijkl constant in space, and an auxiliary strain ®eld, eH kl … x†, to account for the x-dependence of the actual sti€ness matrix. This paper is exclusively concerned with the case of holes and then it suf®ces to put

ÿ Cijkl ˆ 0;

ekl …x† ˆ 0:

…12†

If these equations are inserted into (8)1 a set of equations of the Navier-type results: ‡ Cijkl

o2 o H ‡ uk …x† ˆ Cijkl e …x†: oxj oxl oxj kl

…13†

This, by virtue of the di€erentiation rules shown in Eqs. (6) and (7) can be mapped onto a system of

Fig. 4. Left: vertical cuts through normal stress ®eld, numerical and analytical, for 40, 64 and 80 pixel size holes; right: magni®cation of the transition zone.

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

algebraic equations in Fourier space, so that one particular solution of (13) is given by:

…n‡1† eH ij …x†

  0 eH eij …x† ˆ Y ÿ1 A^‡ mn …x† ‡ eij ; ijmn ^
‡ 1 ÿ ni Njk ‡ nj Nik Cklmn ˆ nl ; 2D Nij ‡ : Mik ˆ Cijkl njl ; Mijÿ1 ˆ D

A^‡ ijmn

…14†

The strain ®eld, e0ij , represents the strain in a homogeneous RVE subjected to external (normal) loads (see Fig. 2). This also accounts for the boundary conditions   ÿ ÿ1
‡ 0 r0x 0 0 0 : …15† rkl ˆ eij ˆ C ijkl rkl ; 0 r0y Unfortunately, this is only a formal solution since the auxiliary ®eld, eH kl , is still unknown. However, it can be determined iteratively, if Eq. (14) is inserted into Eq. (11), n   o ‡ 0 eH ‰1 ÿ h…x†Š Y ÿ1 A^‡ Cijkl rs …x† ‡ ekl klrs ^ n o   ‡ H 0 H ^ … † … † e Y ÿ1 A^‡ x ‡ e ÿ e x : …16† ˆ Cijkl kl rs klrs rs This functional equation is solved by a Neumann iteration:

ÿ
‡ ˆ h…x† C ÿ1 ijmn 8 9 0 1 < = …n† ‡ 0 A eH  Cmnop Y ÿ1 @ A^‡ rs …x† ‡ eop oprs ^ : ;

…0†

eH rs ˆ 0:

191

; …17†

2.3. The solution of the time-independent elastic/ plastic heterogeneity problem It will now be assumed that the matrix material can be described by time-independent bilinear elastic±plastic material theory on the basis of incremental J2 -von Mises ¯ow theory [9]. Consequently, instead of Eq. (11) a di€erential form of HookeÕs law applies   ‡ H … † … † dekl …x† ÿ depl x ÿ de x ; …18† drij …x† ˆ Cijkl kl kl where depl kl …x† denotes the increment in plastic strain. Standard algorithms, such as the radial return method (see [9]), can now be applied for a numerical solution. The ®eld variables are determined for the grid points of the DFT mesh shown in Fig. 2 instead of the GaussÕ points in a ®nite element mesh.

Fig. 5. KirschÕs problem for an elliptical hole and the normal stress in vertical direction.

192

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

Fig. 6. Left: horizontal cuts through normal stress ®eld, numerical and analytical, for a ˆ 20; b ˆ 10; a ˆ 40; b ˆ 20; a ˆ 60; b ˆ 20, and a ˆ 100; b ˆ 20 elliptic holes; right: magni®cation of the transition zone.

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

3. Results 3.1. Kirsch's problem for circular holes in a linearelastic matrix As a ®rst spot check for the accuracy of the stress/strain ®elds obtained from DFT KirschÕs problem of a circular hole in a large (in®nite), linear-elastic matrix subjected to tension at in®nity is considered: Fig. 3, left inset. The inset on the right shows the normal stress ®eld in horizontal direction as it results from a DFT analysis for an RVE of 1024  1024 pixels. The sequence of pictures in Fig. 4 presents central vertical cuts through that stress ®eld and compares them to data obtained from the corresponding analytical solution [10]. It turns out that the numerical solution tends to underestimate the peripheral stresses, especially if larger holes are considered. This should be attributed to a hole-to-hole interaction between the neighboring RVEs and to the numerical representation of the shape function, h…x†. 3.2. Kirsch's problem for elliptical holes in a linearelastic matrix Fig. 5 presents KirschÕs problem for an elliptical hole (left inset) and the normal stress ®eld in ver-

193

tical direction obtained from DFT for an RVE of 1024  1024 pixels. As before a central (horizontal) cut is applied and the corresponding stress data along this line are compared with data from the analytical solution [11]. The ®rst two rows investigate the case of a constant ratio a/b, whereas the last two focus on large, slender ellipses with more pixels in the interior. This was done as an attempt to increase the degree of accuracy. As for the circular hole the numerical solution underestimates the stresses at the periphery of the elliptical hole, and it appears that some interaction between the RVEs is present. 3.3. Grith cracks and stress intensity factors It is tempting to investigate whether the stress information from slender elliptical holes, such as the ones shown in Fig. 6, can be used to obtain stress intensity factors (SIFs). It has been pointed out that DFT tends to underestimate the stresses in the immediate vicinity of the ``crack tip''. Hence it is advisable to use all the stress information within a certain region rather than a single point. A similar approach has been suggested in photoelasticity [12], the digital multiple-point method: A Williams expansion is evaluated for each pixel within a ring around the crack tip (Fig. 7, left inset), and the resulting over-determined linear

Fig. 7. Half circle segment used for SIF evaluation and SIF results.

194

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

system of equations for the unknown coecients, ai , is solved through Householder transforms: rxx …r; u† ‡ ryy …r; u† r0 ÿ2 cos …u=2† a1 2 sin …u=2† a2 p p p ÿ p ˆ r0 a r0 a a a  u p paa 4 ‡ 4a3 ‡ 2 cos a r0 2  u p paa aa6 5 a ÿ 24 sin ‡ 4 cos …u†a r0 r0 2 aa7 ÿ …23 sin …u† ‡ 3 sin …3u††a : r0

…19†

The inset on the right of Fig. 7 shows the resulting (normalized) SIFs as they result from DFT in combination with Eq. (19) by evaluation within a half circle segment at di€erent distances, Rmin , from the crack tip, for an equal width of 10 pixels. Obviously, as it is known from photoelasticity experiments, there is a certain critical evaluation region in front of the crack tip for which the SIF come out best (in the present case of a single Grith crack the SIF should be identical to one). Further results for various other crack problems were recently obtained and can be found in [13].

Fig. 8. A quarter of the stress distribution for KirschÕs problem (tensile load beyond yield limits, 300 Mpa, in 1 and 2-direction, FE-simulation (left) vs. DFT (right).

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

3.4. Kirsch's problem for circular holes in an elastic± plastic matrix Fig. 8 presents an overview of the stress distribution for a circular hole in an elastic±plastic matrix subjected to tension in vertical direction. The plot on the left stems from an FE-calculation (ABAQUS), the plot on the right is the corresponding result from DFT. In Fig. 9 horizontal and vertical cuts have been applied to the lines of symmetry for that problem which allows for a direct comparison of stress data. 4. Conclusions and outlook · For KirschÕs problem with a circular hole, the agreement between the numerical and the analytical results is excellent, even if evaluated directly at the hole periphery.

195

· For elliptic holes, the overall quality of the numerically obtained stress data is good and satisfactory at the periphery of the ellipse. · The accuracy of the numerical results for the ellipse depends on its curvature. · The number of pixels within the hole has an in¯uence on the accuracy of the numerical solution. · A certain interaction between neighboring RVEs could be observed for circular as well as for elliptical holes. · The transition from an elliptical hole to a crack is challenging, ®rst results for stress-intensity factors exist. · The multi-parameter pixel techniques used during a photoelastic analysis of crack tip ®elds can directly be used for DFT. As in the experiment, best results are obtained when the evaluation is performed within a critical region in front of the crack tip.

Fig. 9. Stresses along axes of symmetry of KirschÕs problem (tensile load beyond yield limit, ry ˆ 300 MPa plus linear hardening), in 1 and 2-directions.

196

K.P. Herrmann et al. / Computational Materials Science 16 (1999) 186±196

· Di€erent discretization widths for the di€erence quotients in fx; y; zg-direction might lead to a further improvement of the quality of the results. · DFT and time-independent von-Mises ¯ow theory has been used to study KirschÕs problem of a circular hole in an elastic±plastic matrix to obtain ®rst results when compared to FE results. Further improvement of the code is necessary to make DFT competitive to commercial FEprograms as far as the speed of convergence is concerned. Acknowledgements The authors wish to thank the Deutsche Forschungsgemeinschaft (DFG) for ®nancial support through a grant for the research project Zur Analyse des elasto-plastischen Verhaltens von Verbundwerksto€en mit Hilfe der diskreten Fouriertransformation. One of the authors (WHM) also wants to express his gratitude for ®nancial assistance for participation in IWCMM8 through the organizers. References [1] L.H. Schoenlein, M. R uhle, A.H. Heuer, In situ straining experiments of Mg-PSZ single crystals, in: N. Claussen, M. R uhle, A.H. Heuer (Eds.), Advances in Ceramics, vol. 12, Science and Technology of Zirconia II. The American Ceramic Society, Columbus, Ohio, 1984, pp. 275±282.

[2] W. Dreyer, Development of microstructure based viscoplastic models for an advanced design of single crystal hot section components, in: J. Olschewski (Ed.), Development of Microstructural Based Viscoplastic Models for an Advanced Design of Single Crystal Hot Section Components, 1994. [3] P.G. Harris, K.S. Chaggar, M.A. Whitmore, The e€ect of ageing on the microstructure of 60:40 tin-lead solders, Soldering and Surface Mount Technol. 7 (1991) 20±23. [4] A. Zavaliangos, L. Anand, Thermo-elasto-viscoplasticity of porous isotropic metals, J. Mech. Phys. Sol. 41 (6) (1993) 1087±1118. [5] W. Dreyer, W.H. M uller, A study of the coarsening in tin/ lead solders, IJSS, accepted, 1999. [6] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, C.R. Acad. Sci. Paris 318 (II) (1994) 1417± 1423. [7] W.H. M uller, Fourier transforms and their application to the formation of textures and changes of morphology in solids, in: Y.A. Bahei-El-Din, G.J. Dvorak (Eds.), Proceedings of the IUTAM Symposium on Transformation Problems in Composite and Active Materials, Cairo, Kluwer Academic Publishers, The Netherlands, 1998, pp. 61±72. [8] T. Mura, Micromechanics of Defects in Solids, second revised edition, Martinus Nijho€ Publishers, Dordrecht, The Netherlands, 1987. [9] K.-J. Bathe, Finite-Elemente-Methoden, Springer, Berlin, 1987. [10] H.G. Hahn, Elastizit atstheorie, Teubner, Stuttgart, 1985. [11] H.G. Hahn, Bruchmechanik, Teubner, Stuttgart, 1976. [12] F. Ferber, K.P. Herrmann, Experimental and numerical simulations of failure mechanisms in thermomechanically loaded composite material models, in: D. Hui (Ed.), ICCE/ 2, New Orleans, 21±24 August 1995, pp. 225±226. [13] S. Neumann, K.P. Herrmann, W.H. M uller, Computation of stresses and strains in heterogeneous bodies by use of the discrete Fourier transform, GAMM 99, Metz, 12±16 April 1999.