A numerical method for computing the overall response of nonlinear composites with complex microstructure
AA numerical method for computing the overallresponse of nonlinear composites with complexmicrostructure
H. Moulinec, P. Suquet
L.M.A./ C.N.R.S.31 Chemin Joseph Aiguier13402. Marseille. Cedex 20. France.
Abstract
The local and overall responses of nonlinear composites are classically investigated by the Finite Ele-ment Method. We propose an alternate method based on Fourier series which avoids meshing and whichmakes direct use of microstructure images. It is based on the exact expression of the Green function ofa linear elastic and homogeneous comparison material. First the case of elastic nonhomogeneous con-stituents is considered and an iterative procedure is proposed to solve the Lippman-Schwinger equationwhich naturally arises in the problem. Then, the method is extended to nonlinear constituents by a step-by-step integration in time. The accuracy of the method is assessed by varying the spatial resolution of themicrostructures. The flexibility of the method allows it to serve for a large variety of microstructures.
This study is devoted to a numerical method introduced by Moulinec and Suquet [1], [2] to determine thelocal and overall responses of nonlinear composites. Numerous studies dealt with nonlinear cell calculationsby the Finite Element Method (FEM) (see for example Adams and Donner [3], Christman et al [4],Tvergaard [5], Michel and Suquet [6]). Most of them are limited to “simple” microstructures, one ortwo inclusions embedded in a volume of matrix. The need to incorporate more detailed information on themicrostructure is clearly recognized. Recently, several studies have considered “complex” microstructuresinvolving a significant number of inclusions with irregular shape. Brockenborough et al [7], B¨ohm et al [8],Nakamura and Suresh [9], Dietrich et al [10], Becker and Richmond [11] are some of the contributions tothis recently developed subject. All were based on the FEM. The difficulties due to meshing and to the largenumber of degrees of freedom required by the analysis limit the complexity of the microstructures which canbe investigated by this method.A typical example of a complex microstructure which is difficult to mesh and therefore to handle by meansof the FEM is shown in Figure 10 taken from the work of Bornert [12]. The digital image of this Iron/Silverblend was obtained by Scanning Electron Microscopy (SEM). The initial idea of the method proposed in[1] was to make direct use of these digital images of the real microstructure in the numerical simulation. Asimilar idea can be found in Garboczi and Day [13] who used a spring network technique.The proposed method avoids the difficulty due to meshing. It makes use of Fast Fourier Transforms (FFT)to solve the unit cell problem , even when the constituents have a nonlinear behavior. FFT algorithmsrequire data sampled in a grid of regular spacing, allowing the direct use of digital images of the microstruc-ture. The second difficulty (size of the problem) is partially overcome by an iterative method not requiringthe formation of a stiffness matrix.The interest in numerical simulations of the nonlinear response of composites has recently been strength-ened by the development of theoretical methods which analytically predict the nonlinear overall behavior ofcomposites ( Willis [15], Ponte Casta˜neda [16], Suquet [17]). Part of the present study provides precisenumerical results for uniaxial loadings which could serve as guidelines for theoretical predictions.The body of the method and the resulting algorithms are presented in section 2. In section 3, the accuracy ofthe method and several numerical points are discussed (choice of the reference medium, spatial resolution ....).In section 4 the method is applied to determine the local and overall responses of composites with ”random”microstructures. In all the cases considered in this study the models have been limited to two dimensional During the revision of this paper, the attention of the authors was called on a similar work by M¨uller [14] concerning phasetransformation. a r X i v : . [ c s . C E ] D ec pproximations. The first reason for this approximation is the limitation on current computational capability.The second reason is that many microstructural observations are two dimensional. The overall behavior of a composite is governed by the individual behavior of its constituents and by its mi-crostructure. Its effective response to a prescribed path of macroscopic strains or stresses may be determinednumerically via the resolution of the so-called ”local problem” on a representative volume element (r.v.e.) V .In this study, the ”representative” information on the microstructure is provided by an image (micrograph)of the microstructure with arbitrary complexity. The image contains N pixels, and independent mechanicalproperties are assigned individually to each pixel. Most applications involve only a limited number of phases,although in principle each pixel could be considered as an individual constituent.The local problem consists of equilibrium equations, constitutive equations, and boundary and interfaceconditions. All different phases are assumed to be perfectly bonded (displacements and tractions are contin-uous across interfaces). Displacements and tractions along the boundary of the r.v.e. are left undeterminedand the local problem is ill-posed. We choose to close the problem with periodic boundary conditions whichcan be expressed as follows. The local strain field ε ( u ( x )) is split into its average E and a fluctuation term ε ( u ∗ ( x )): ε ( u ( x )) = ε ( u ∗ ( x )) + E or equivalently u ( x ) = u ∗ ( x ) + E . x . By assuming periodic boundary conditions it is assumed that the fluctuating term u ∗ is periodic (notation: u ∗ σ . n is anti-periodic in order to meet the equilibrium equations on the boundarybetween two neighboring cells (notation: σ . n − First we consider the preliminary problem of a homogeneous linear elastic body with stiffness c subjectedto a polarization field τ ( x ). σ ( x ) = c : ε ( u ∗ ( x )) + τ ( x ) ∀ x ∈ V div σ ( x ) = ∀ x ∈ V, u ∗ , σ . n − (1)The solution of (1) can be expressed in real and Fourier spaces, respectively, by means of the periodic Greenoperator Γ associated with c : ε ( u ∗ ( x )) = − Γ ∗ τ ( x ) ∀ x ∈ V, (2)or ˆ ε ( ξ ) = − ˆ Γ ( ξ ) : ˆ τ ( ξ ) ∀ ξ (cid:54) = , ˆ ε ( ) = (3)The operator Γ is explicitly known in Fourier space (see appendix A). When the reference material isisotropic (with Lam´e coefficients λ et µ ) it takes the form :ˆΓ ijkh ( ξ ) = 14 µ | ξ | ( δ ki ξ h ξ j + δ hi ξ k ξ j + δ kj ξ h ξ i + δ hj ξ k ξ i ) − λ + µ µ ( λ + 2 µ ) ξ i ξ j ξ k ξ h | ξ | . (4) The auxiliary problem can be used to solve the problem of an inhomogeneous elastic composite materialwith stiffness c ( x ) at point x under prescribed strain E : σ ( x ) = c ( x ) : (cid:0) ε ( u ∗ ( x )) + E (cid:1) ∀ x ∈ V div σ ( x ) = ∀ x ∈ V, u ∗ , σ . n − (5)For simplicity E is assumed to be prescribed, although other average conditions could be considered aswell (see appendix B for prescribed stresses). A homogeneous reference material with elastic stiffness c isintroduced and a polarization tensor τ ( x ), which is unknown a priori , is defined as : τ ( x ) = δ c ( x ) : ε ( u ( x )) , δ c ( x ) = c ( x ) − c . (6)Thus, the problem reduces to the periodic Lippmann-Schwinger equation ( Kr¨oner [20]), which reads, in realspace and Fourier space respectively: ( u ( x )) = − Γ ( x ) ∗ τ ( x ) + E , (cid:98) ε ( ξ ) = − (cid:98) Γ ( ξ ) : (cid:98) τ ( ξ ) ∀ ξ (cid:54) = , (cid:98) ε ( ) = E (7)where τ is given by (6). The Lippman-Schwinger equation is an integral equation for ε ( u ∗ ). The principle of the algorithm is to use alternately (6) and (7), in real space and Fourier space, respectively,in an iterative scheme, to solve (5):
Initialization : ε ( x ) = E , ∀ x ∈ V, σ ( x ) = c ( x ) : ε ( x ) , ∀ x ∈ V,Iterate i + 1 : ε i and σ i being known a ) τ i ( x ) = σ i ( x ) − c : ε i ( x ) ,b ) (cid:98) τ i = F ( τ i ) ,c ) Convergence test ,d ) (cid:98) ε i+1 ( ξ ) = − (cid:98) Γ ( ξ ) : (cid:98) τ i ( ξ ) ∀ ξ (cid:54) = and (cid:98) ε i+1 ( ) = E ,e ) ε i+1 = F − ( (cid:98) ε i ) ,f ) σ i+1 ( x ) = c ( x ) : ε i+1 ( x ) . (8) F and F − denote the Fourier transform and the inverse Fourier transform. This algorithm can be furthersimplified by noting that Γ ∗ ( c : ε ) = ε . The modified algorithm reads :
Initialization : ε ( x ) = E , ∀ x ∈ V, σ ( x ) = c ( x ) : ε ( x ) , ∀ x ∈ V,Iterate i + 1 : ε i and σ i being known a ) ˆ σ i = F ( σ i ) ,b ) Convergence test ,c ) ˆ ε i+1 ( ξ ) = ˆ ε i ( ξ ) − ˆ Γ ( ξ ) : ˆ σ i ( ξ ) ∀ ξ (cid:54) = and ˆ ε i+1 ( ) = E ,d ) ε i+1 = F − (ˆ ε i+1 ) e ) σ i+1 ( x ) = c ( x ) : ε i+1 ( x ) , ∀ x ∈ V, (9)Convergence is reached when σ i+1 is in equilibrium. The error serving to check convergence is : e i = (cid:0) < || div( σ i ) || > (cid:1) / || < σ i > || = (cid:0) < || ξ . ˆ σ i ( ξ ) || > (cid:1) / || ˆ σ i ( ) || . The iterative procedure is stopped when the error e is smaller than a prescribed value (typically 10 − in ourcalculations). The unit cell is discretized into a regular grid consisting of N × N pixels (two-dimensional problem),or N × N × N ”voxels” (tri-dimensional problem). The data and the unknowns used in the numericalcalculations are images sampled on this grid ( N × N or N × N × N arrays). In two dimensions, thecoordinates of the pixel labeled by i , i are x d ( i , i ) = (cid:18) ( i − · T N , ( i − · T N (cid:19) , i = 1 , ...N , i = 1 , ...N , where T j is the period of the unit cell in j th direction ( j = 1 , Initialization : ε ( x d ) = E , ∀ x d ∈ V, σ ( x d ) = c ( x d ) : ε ( x d ) , ∀ x d ∈ V,Iterate i + 1 : ε i and σ i known at every x d a ) ˆ σ i = FFT ( σ i ) ,b ) Convergence test ,c ) ˆ ε i+1 ( ξ d ) = ˆ ε i ( ξ d ) − ˆ Γ ( ξ d ) : ˆ σ i ( ξ d ) ∀ ξ d (cid:54) = , ˆ ε i+1 ( ) = E ,d ) ε i+1 = FFT − (ˆ ε i ) e ) σ i+1 ( x d ) = c ( x d ) : ε i+1 ( x d ) , ∀ x d ∈ V (10)where x d denote the coordinates of pixels in real space, and ξ d denote the N × N corresponding frequenciesin Fourier space. To be more specific, the discrete frequencies are (in dimension 2) when N j is even : ξ j = ( − N j T j , ( − N j T j , ..., − T j , , T j , ..., ( N j −
1) 1 T j , N j T j , and when N j is odd : ξ j = − N j −
12 1 T j , ..., − T j , , T j , ..., N j −
12 1 T j . The discrete error serving to check convergence is : e i = (cid:32) N (cid:88) d || ξ d · ˆ σ i ( ξ d ) || (cid:33) / || ˆ σ i ( ) || (where N = N × N is the total number of pixels).When the spatial resolution is low and when the number N j of discretization point is even, a specialattention must be paid to the highest frequency ξ j = ± (cid:16) N j (cid:17) T j , j = 1 or 2. In most FFT packages,the Fourier expansion at these frequencies consists of either cos( ξ j x j ) or exp( − i ξ j x j ), instead of the correctexpression consisting of the two terms exp( − i ξ j x j ) and exp(i ξ j x j ). Therefore, even when the stress σ iscorrectly approached by its Fourier expansion in step a) of the algorithm (10), the result of step d) may notapproach accurately the Fourier expansion of the strain ε at these particular frequencies. This is because ˆ Γ is neither even nor odd with respect to each individual component ξ j . Oscillations were observed when (4)was used with relatively small values of N j (lower than 128). This problem was fixed by using a differentexpression of ˆ Γ in algorithm (10) at these frequenciesˆ Γ = (cid:0) c (cid:1) − . In other terms, the stress σ is forced to by the algorithm at these frequencies when convergence is reached. The algorithm can be extended to the case in which the individual constituents obey a nonlinear law, writteneither in terms of stresses and strains (nonlinear elasticity at infinitesimal strain) or in incremental formrelating strain-rates and stress-rates (flow theory). The nonlinearity requires an appropriate modificationof step e ) in algorithm (10). In the present study, special attention will be paid to phases exhibiting anincremental elastic-plastic behavior at small strains governed by a J -flow theory with isotropic hardening(although more general constitutive laws can be considered) :˙ σ = c : ( ˙ ε − ˙ ε p ) , ˙ ε p = ˙ p s σ eq , σ eq − σ ( p ) ≤ , ˙ p ≥ . (11) ε p denotes the plastic strain, s denotes the stress deviator and p denotes the hardening parameter, whichcoincides with the cumulated plastic strain˙ p ( t ) = (cid:18)
23 ˙ ε pij ( t ) ˙ ε pij ( t ) (cid:19) , p ( t ) = (cid:90) t ˙ p ( s ) ds σ eq = (cid:18) s ij s ij (cid:19) . An error had crept into the expression of the highest frequency in the original paper published in Comput Methods ApplMech Eng. The authors thank Anthony Rollett for pointing it out. HM 11/12/2020 he integration in time of the constitutive law (11) is achieved by means of an implicit scheme which isclassical in the analysis of elastic-plastic structures by the FEM method. The time interval (or, alternatively,the loading path) is discretized into subintervals [ t n , t n +1 ]. The field equations are solved for ( ε n , σ n , p n ),which denote strain, stress and hardening parameter at time t n . Assuming that these fields are known atstep n (time t n ), the principal unknown at step n + 1 is ε n +1 . The incremental equations (11) are discretizedby an implicit scheme. The unknown ε n +1 is a compatible strain field such that the associated stress field (bythe constitutive law) is in equilibrium. The resulting system of equations to be solved for ε n +1 is nonlinear.The algorithm for the determination of ε n +1 reads (for simplicity the lowerscript ( n + 1) is omitted below;superscripts i and i+1 refer to the iterative loop within the step) : Initialization : ε ( x d ) given by (13) , Compute σ and p from ( ε , σ n , ε n , p n ) ,Iterate i + 1 : ε i and σ i are known a ) ˆ σ i = FFT ( σ i ) ,b ) Convergence test ,c ) ˆ ε i+1 ( ξ d ) = ˆ ε i ( ξ d ) − ˆ Γ ( ξ d ) : ˆ σ i ( ξ d ) ∀ ξ d (cid:54) = , ˆ ε i+1 ( ) = E n +1 ,d ) ε i+1 = FFT − (ˆ ε i+1 ) e ) Compute σ i+1 and p i+1 from ( ε i+1 , σ n , ε n , p n ) (12)More specificallya) The initial strain ε at time t n +1 is extrapolated (linearly) from ε n and ε n − at the two previous timesteps t n and t n − : ε ( x d ) = ε n ( x d ) + t n +1 − t n t n − t n − ( ε n ( x d ) − ε n − ( x d )) , ∀ x d ∈ V. (13)This choice significantly improves the convergence of the iterative process within the time step.b) σ i and p i are computed from ( ε i , σ n , ε n , p n ) (step e ) in algorithm (12)) by a radial return method (seeappendix C). The rate of convergence of the algorithm depends drastically on the Lam´e coefficients λ and µ of thereference material. After several tests, the best rate of convergence was observed with λ = (cid:18) inf x ∈ V λ ( x ) + sup x ∈ V λ ( x ) (cid:19) µ = (cid:18) inf x ∈ V µ ( x ) + sup x ∈ V µ ( x ) (cid:19) (14)The number of iterations at convergence is significantly influenced by several other parameters. First,as shown in Figure 1, it increases with the contrast between the phases (typically the ratio between theelastic moduli of the phases). When the contrast is infinite (rigid inclusions or voids in an elastic matrix),the algorithm no longer converges. Second, the number of iterations at convergence also depends on thecomplexity of the solution itself. In the example of an elastic ideally plastic matrix reinforced by stiffinclusions, the computing time increases with the tortuosity of the bands where the strain tends to localize(see below). The constitutive law acts locally in real space ( i.e. applies separately to each individual point x ). Similarly,Green’s function Γ acts locally in Fourier space, ( i.e. applies separately to each individual frequency ξ ). From a computational standpoint, the corresponding steps (c and e in the algorithms (10) or (12) )are performed by independent loops on each individual pixel in real or Fourier space. These steps canconsequently be vectorized or parallelized. In addition, optimized FFT packages are available on most vectoror parallel computers. The whole algorithm can therefore be efficiently implemented on these machines.It follows from the same argument that the time spent in the steps corresponding to the constitutive lawand to the Lippman Schwinger equation varies linearly with the number N of pixels. The CPU time for aFT varies as N · log N . The time required by the other steps of the algorithm are comparable to the timerequired by the FFTs. The CPU time t for one iteration can be estimated by k × N ≤ t ≤ k × N log N, where k and k are expected to be independent of the size N of the problem. The dependence of the CPUtime on the size of the problem is shown in Figure 2. The square unit cell shown in Figure 5 is subjectedto uniaxial transverse tension at 0 . The volume fraction of fibers is 47.5%. Both the fibers and the matrixare assumed to be elastic with elastic constants given by (17) and (18). The dependence of the CPU timeon the size of the problem is approximately linear. Optimizing the memory occupancy . The Fourier transform of a real valued function has the symmetryproperty ˆ f ( − ξ ) = ˆ f ( ξ ) . Since all quantities under consideration in our computation are real, this symmetry property allows us torestrict our attention to positive frequencies (the values of the fields for negative frequencies being immedi-ately deduced). The size of the arrays can therefore be divided by 2, provided the FFT package allows forthe storage of real numbers as complex numbers with the same memory occupancy.
Performances . Most computations were run on a Cray YMP with peak performance of 333
M F lops . Theperformance observed with our algorithm was (cid:39)
M F lops on the elastic-plastic problem described insection 4 with unit cells discretized into 1024 × pixels . The typical CPU time on one processor of thiscomputer is less than 30 seconds for an elastic problem (with a spatial resolution of 1024 × pixels ,the ratio between the Young moduli being approximately 6). When the matrix is elastic plastic, the typicalCPU time for a run as described in section 4 is 4000 seconds. To assess the accuracy and the stability of the method we examined two cases for which analytical solutionsare available.
Laminates . The first example concerns layered materials. As is well-known, the strain field is then uniformwithin each individual layer and takes different values from one layer to another. The example shownin Figure 3 corresponds to a two-phase material, both phases having equal volume fraction. The layersare parallel to the plane ( x , x ). The constitutive materials of the layers were linear elastic with elasticcharacteristic given by (17) and (18). The applied loading was pure shear parallel to the layersΣ arbitrary , Σ = Σ = Σ = Σ = Σ = 0 . The image was discretized into 32 ×
32 pixels (good results were obtained with an even cruder resolution).The computed local strain field ε is plotted in Figure 3 and shows no oscillation. In addition the numericalsolution coincides with the exact solution. Circular fiber at dilute concentration . The second example concerns the elastic strain field generatedby stiff circular fibers placed at the nodes of a square lattice in a more compliant matrix. The exact solutionto this problem (with periodic boundary conditions) is not known in closed form (to the authors’ knowledge).However when the volume fraction of fibers is small this solution can be accurately approximated by thesolution of a simpler problem, where a circular fiber (with radius a ) is surrounded by a circular shell ofmatrix (with radius b ) and subject to the boundary condition u ( x ) = E . x when r = b, where the overall strain E is the same as in the original periodic problem. When the imposed loading is anin-plane shear E (cid:54) = 0, other E ij = 0, the displacement field has the form u r ( r, θ ) = (cid:0) Ar + Br + C ø r + D ø r (cid:1) sin(2 θ ) ,u θ ( r, θ ) = (cid:0) l + 3 µ ø lAr + Br + µ ø l + 2 µC ø r − D ø r (cid:1) cos(2 θ ) , where r and θ are the polar coordinates in the plane. A, B, C, D, take different values in the matrix and inthe fiber. They solve a system of linear equations expressing the boundary condition at r = b , the absenceof singularity at r = 0, the continuity of tractions and displacements at r = a .According to Saint Venant’s principle, the local strain fields in the two problems coincide far from theboundary of the cell. Therefore at low volume fraction of fibers ( a /b (cid:28) a/b = 1 /
16. The spatial discretization used in the numerical calculation was1024 × ε of the strain field in a square window of width c = 4 a is shown in Figure4 (note that the unit cell itself with width 2 b is much larger than the window shown). There is almost nodifference between the analytical and the numerical solutions shown in (a) and (b) respectively. A moreexplicit comparison is made in Figure 4 (c) which shows an horizontal cut through the field ε at x = 0.Except from little undulations inside the inclusion, there is no significant oscillations at the fiber boundarywhere the field ε is discontinuous. In addition the accuracy of the numerical solution is observed to increasewith the spatial resolution. The discrepancy between the numerical and the analytical solutions dependson the spatial resolution and should not be attributed to a Gibbs phenomenon, i.e. to an oscillation ofthe Fourier series of a function in the vicinity of a discontinuity point. This oscillation is attached to thesummation of the Fourier series which is not what the discrete inverse Fourier transform performs. Discrete Fourier transform . The discrete Fourier transform, when applied to an image discretized into N × N pixels, is the exact Fourier transform of the image when two requirements are met : ( Brault andWhite [22]) C1 the image is periodic with the same period ( T , T ) as the unit cell, C2 the image cut-off frequency f c ( i.e. the frequency above which the Fourier transform of the imagevanishes identically) is less than half of the sampling frequency (Shannon’s theorem): f cj < N j T j j = 1 , As already stated the influence of the spatial resolution depends on the stress and strain gradients withinthe phases and therefore on the strength of the phases nonlinearities. The following examples illustrate thesegeneral considerations. The method has been applied to simulate the local and overall response of compositesreinforced by unidirectional long fibers aligned along the e direction. The geometry of these composites isdescribed by a two-dimensional image of their cross section. Generalized plane strains were assumed : u ( x ) = u ( x , x ) , u ( x ) = u ( x , x ) , u ( x ) = E x . (15)The overall strain E has four independent components E , E , E , E (the other two are equal to 0).The overall stress Σ also has four independent components. It is possible to prescribe either a path inthe space of strains, or a path in the space of stresses, or alternatively some components of the strain andthe other components of the stress. Classical plane strains are a particular case of the more general settingconsidered in (15). It corresponds to a path in the space of strains along which E is identically 0. The needto introduce generalized plane strain is illustrated by uniaxial tension in the 0 direction, which correspondsto a path in the space of stresses along whichΣ arbitrary , Σ = Σ = Σ = 0 . (16)The axial component E of the strain is unknown and determined a posteriori by the condition Σ = 0.The assumption of generalized plane strains reduces (5) to a two-dimensional problem for the two unknowns( u ∗ , u ∗ ).Two classical configurations were investigated in which the fibers were placed at the nodes of a square orhexagonal lattice. The fibers were assumed to be elastic, isotropic, and characterized by a Young modulusand a Poisson ratio : E f = 400 GPa , ν f = 0 . . (17)The fiber volume fraction was 47.5 % (for comparison, we chose the same volume fraction as in [8]). Thebehavior of the matrix was varied from linear elasticity to elasto-plasticity with hardening so as to study theeffect of the nonlinearity on the accuracy of the method. All the constitutive laws of the matrix which wereonsidered can be put in the incremental form (11). Its isotropic elastic properties were characterized by aYoung’s modulus and Poisson coefficient E m = 68 . , ν m = 0 . . (18)The plastic properties of the matrix were governed by the Von Mises criterion σ eq ≤ σ + Hp. (19)The initial yield stress σ was either infinite (pure linear elasticity) or given by σ = 68 . H was either 0 (perfectly plastic behavior) or H = 1 171 MPa (isotropic linear hardening).The influence of spatial resolution on the accuracy of the results was studied. The spatial resolution of theimage is determined here through the square root of the total number of pixels contained in the image dividedby the number of fibers in the image. For the square array, with N × N pixels and a single fiber in the unitcell, the spatial resolution is exactly N . The hexagonal array can be viewed as a rectangular array, thusallowing the use of the Fourier technique in orthogonal coordinates, instead of the natural nonorthogonalcoordinates defined by the two unit vectors of the hexagonal lattice (see Figure 5). The rectangular unit cellcontains 1 + 4 × = 2 fibers. The number of pixels along the first direction x is 2 times larger than thenumber of pixels in the second direction x . The spatial step in x is 2 √ / x .Therefore in the hexagonal array, the spatial definition as defined above is again N for an image containing2 N × N pixels.Both unit cells were submitted to uniaxial tension at 0 and 45 in the sense of (16). The results ofthe overall response of the composite are shown in Tables 1 to 6. The initial response of the composite islinear and its slope defines the overall Young’s modulus of the composite. When the matrix is elastic ideallyplastic the overall stress applied to the composite in the direction of tension reaches (asymptotically) a limitwhich defines the overall flow stress of the composite. When the matrix is governed by a linear hardening,the stress-strain curve of the composite exhibits a nonlinear transition to an asymptotically linear (affine)response. The slope of this limit response is the overall hardening modulus of the composite.Each table gives an overall material constant as a function of the spatial resolution of the image. The”error” was estimated as the relative difference between the result at a given resolution and the result at thefinest resolution.These results suggest the following remarks.1. When both constituents are linearly elastic, the overall stiffness is not very sensitive to spatial reso-lution. Even at the lowest resolution (32 ×
32 pixels/fiber), the estimated error was under 1% in allcases.2. When the matrix is elastic plastic, the local and overall responses are sensitive to spatial resolution.The strain fields exhibit a strong tendency to concentrate in thin bands. The higher the nonlinearity,the thinner the bands. These stiff gradients in strain require high spatial resolution to be correctlycaptured.3. The solutions may even be discontinuous when the matrix is elastic-perfectly plastic. This explains therelatively high errors at low resolution: about 15% for the square array of fibers in an elastic-perfectlyplastic matrix under tension at 0 ◦ , with a resolution of 32 ×
32 pixels/fiber. Shear bands can form inthe matrix under tension at 45 ◦ . These shear bands correspond to a mode of deformation of the r.v.e.in plane strains. Therefore, for this particular loading, the effective behavior of the composite dependsonly on the behavior of the matrix. The overall flow stress of the composite coincides with the flowstress of the matrix under plane strain conditions, i.e. σ √ . The formation of a slip plane through thematrix is well captured by the numerical method and explains the precision of the numerical result forthis particular loading.4. When the matrix has linear hardening, the strain fields are more regular than in the perfectly plasticcase. The local and overall responses of the composite are less sensitive to spatial resolution. The erroron the hardening modulus is about 7.5% with a resolution of 32 ×
32 pixels/fiber.This study of the influence of spatial resolution led us to use a resolution of 128 ×
128 pixels/fiber in mostof the examples presented in the next section.
Fiber arrangement
In this section we investigate the influence of the geometrical arrangement of the fibers on the local andoverall responses of nonlinear composites. Attention is again restricted to two-dimensional problems, i.e. tocomposites reinforced by aligned fibers. The fiber arrangement is determined by a two-dimensional image ofthe composite cross section.
Two classes of fiber arrangement, regular and random, were considered. The fibers were identical circulardisks and they were not allowed to overlap (impenetrability condition) except in section 4.3. In mostsimulations the fiber volume fraction was prescribed to 47 . Standard fiber distribution . The “standard” configurations consist of a single fiber placed at the nodesof a square or an hexagonal lattice (see preceeding section). Most F.E.M. cell calculations reported in theliterature are based on these standard configurations with the exceptions of Brockenborough et al (1991) andB¨ohm et al (1993) who investigated the effect of disorder in the fiber arrangement on the overall transverseproperties of composites.
Random fiber distribution
In the ”random” configurations, the centers of the fibers were placed atrandom in the unit cell, subject only to the constraints of impenetrability and periodicity. The latterconstraint implies that, when a fiber overlaps the boundary of the unit cell, it is split into two parts I andII (see Figure 6 ) to fit in the unit cell. The size of the images was the largest one allowed by the memoryon our computer and compatible with a resolution of 128 ×
128 pixels per fiber. These two constraints ledto unit cells discretized into 1024 × Twenty three different configurations of 64 impenetrable fibers were generated randomly in the unit cell.The fibers were assumed to be elastic with material properties given by (17). The matrix was an elasticplastic material governed by a J flow theory (11) with material properties given by (18) (19). The localand global responses of each configuration to a transverse uniaxial tension in the 0 direction (according to(16)) were computed with the above described method. The square array and hexagonal array were alsosubjected to transverse tension in the 0 and 45 directions. The stress-strain curves predicted by the simulation are shown in Figure 7. The solid line corresponds tothe mean response (average of the stress-strain curves over the 23 configurations).These results call for the following comments :1. The fibers were stiff and perfectly bonded to the matrix. Therefore, although the strain E in theaxial direction was not imposed a priori (Σ was prescribed to 0), it was relatively small along thewhole loading path. The strain state was consequently close to the plane strain state, explaining thestrain concentrations observed in the perfectly plastic matrices. As is well known, plane strain is morefavorable to these strain concentrations than is pure uniaxial tension.2. The square lattice has a marked transverse anisotropy which is strengthened by the nonlinear behavior,which gives raise to different responses when the direction of tension makes an angle of 0 or 45 withone of the axes of the square lattice. The low value of the flow stress in the diagonal direction (45 ) isdue to a shear plane passing through the matrix. Indeed, when a plane of shear can be passed throughthe weakest phase of a composite, the shear strength of the composite is exactly the strength of theweakest phase ( Drucker (1959)). In tension (under plane strains) in a direction inclined at 45 on thisplane, the transverse flow stress of the composite is 2 σ m / √
3. This is the flow stress observed in Figure7 and Table 4 (2 σ m / √ (cid:39) .
56 MPa). In conclusion, except at low volume fractions, the squarearray should not be used to investigate the transverse properties of transversely isotropic nonlinearcomposites.3. The hexagonal lattice approaches transverse isotropy. When the matrix is a hardening material, thepredictions obtained with the hexagonal lattice underestimate the stiffness of the composite, or atleast are located below the average of the predictions for the random configurations in the range ofoverall deformations considered. Another computation, not reported here, was performed up to 30%of transverse strain, with no modification in the conclusions. A similar observation was made byBrockenborough et al (1991) for another system. When the matrix is ideally plastic, the low value ofhe flow stress in the diagonal direction (45 ) is again due to a shear plane passing through the matrix.In conclusion, the hexagonal lattice should be used with care to predict the transverse properties ofnonlinear composite systems, even for hardening matrices.4. The deviation from the average of the transverse Young’s moduli computed on the different configu-rations is small. By contrast, the deviations in the other properties (flow stress, hardening modulus)are higher and may be attributed to the combined effects of nonlinearity and incompressibility.5. The local plastic strains showed significant differences between the ideally plastic case and the hardeningcase. For the former, the strain concentrates in thin bands in the matrix. In most configurations, onlya small percentage of the matrix contributes to the plastic dissipation. The overall flow stress ofthe composite is observed to be directly related to the ”tortuosity” of these bands. Two differentconfigurations with the corresponding zones of strain concentration are shown in Figure 8. In thefirst configuration slip bands inclined at approximately 45 on the direction of traction can be passedthrough the matrix, resulting in a low flow stress. Conversely, the fiber arrangement in the seconddirection inhibits long-range slip bands and causes these bands to deviate or the plastic deformationto spread into wider zones. The plastic dissipation and the flow stress are higher in the secondconfiguration than in the first one. Adding more fibers in the undeformed zones would not change theplastic dissipation, or in other terms, would not affect the flow stress of the composite. These resultslead us to think that, when the matrix is perfectly plastic, the geometrical parameter which governs(at first order) the flow stress of the composite is not the volume fraction of the fibers but, instead, thelength of the shortest path passing through the matrix at an angle of approximately 45 in tension, or0 in shear.6. When the matrix is a hardening material, the plastic strain spreads all over the matrix (see Figure8). The whole matrix contributes (although non homogeneously) to the plastic dissipation and, con-sequently, to the overall strengthening of the composite. In this case, the volume fraction of the fibersseems to be the relevant geometrical information (at least to first order) to predict the overall hardeningof the composite.7. In spite of the differences in the maps of plastic strains in the ideally plastic material and in thehardening matrix, the ”stiffest” (respectively the ”weakest”) configurations in the ideally plastic caseremain the stiffest (respectively the weakest) configurations in the hardening case. The present section deals with the ”representativity” of a unit cell in two aspects. First, does the unit cellcontain enough heterogeneities so that the computed effective properties no longer depend on the cell size?Second, how much do different unit cells randomly generated with the same volume fraction and number ofheterogeneties differ from each other?Several series of microstructures containing 4, 9, 16, 36, 64 or 256 impenetrable fibers randomly placedin the unit cell were generated. The volume fraction of fibers was identical in all simulations (47 . ×
128 pixels/fiber). The total number of pixels in each imagewas therefore the number of fibers multiplied by × × (see (16)). Statistical data on the computed Young’s modulias a function of the number of fibers in the unit cell are reported in Table 7. The mean Young’s modulusand its standard deviation are defined as¯ E = 1 N s (cid:88) i =1 ,N s E i , σ ( E ) = (cid:115) N s − (cid:88) i =1 ,N s ( E i − ¯ E ) where E i is the Young’s modulus of the i th microstructure and N s is the number of different microstructures.The error on the mean is classically estimated by the ratio σ ( E )¯ E √ N s Similar data on the overall flow stress of the composite are given in Table 8. The number of fibers in theunit cell does not significantly influence the mean overall properties, provided a lower number of fibers iscompensated by a higher number of configurations. The mean Young’s modulus and the mean flow stressof configurations with four fibers differ from those of configurations with 256 fibers by 0 .
56% and 0 . .
13% and 0 .
23% for theYoung’s modulus and the flow stress for configurations with 256 fibers). This is an illustration of the ergodicroperty : spatial averaging on one large sample is equivalent to ensemble averaging on many small samples.A related observation is that the standard deviations of the overall properties decrease as the number offibers increases.
In the above analyses, the fibers were placed randomly in the unit cell with impenetrability as the onlyrestriction. The effects of imposing a minimal space between fibers are of interest for at least two reasons.First, when the minimal spacing between the centers of the fibers increases, the ordering of the microstructureincreases. As a limit case, when this minimal spacing reaches (cid:113) √ . SN ( S is the surface of the unit cell,N is the number of fibers), the microstructure is completely determined and coincides with the centeredhexagonal arrangement. Second, numerical difficulties could be expected when two neighboring fibers arenearly touching. Indeed, when the spatial resolution is not fine enough, the method cannot capture the highstrain gradients in the necks between the two fibers.Ten configurations with 64 fibers were generated, and a minimal space of 4 pixels between two neighboringfibers was imposed. This distance seemed sufficient to correctly describe strain concentration. The resultsof this study suggest the following comments :1. When the matrix is elastic ideally plastic, the mean overall flow stress is Σ = 86 . .
44 MPa). This value is 2 .
0% smaller than the value obtained with no restrictionon the space between fibers. It lies slightly below the flow stress of the hexagonal array subjected totension at 0 (Σ = 87 . or 45 (Σ = 79 . .2. When the matrix is elastic-plastic with linear hardening, the effective hardening modulus drops signif-icantly : H = 9382 MPa (estimation error = 123 . H = 7100 MPa at 0 , H = 7420 MPa at 45 ).In conclusion, it seems that the “safety coating” around the fibers leads to a decrease in the overall mechanicalproperties of the composite, at least at the volume fraction which has been investigated. The above analyses show that the overall flow stress of the composite and, to a lesser extent, its overallhardening depend primarily on the tortuosity of shear bands passing through the matrix. Obviously, thevolume fraction of the reinforcing phase plays a role in the possibility that such bands are formed, but fora fixed volume fraction, significant differences arise from the differences in the patterning of bands. Theseshear bands are locked or deviated by the fibers. The overall flow stress of the composite can (empirically) berelated to the length of the shortest path passing through the matrix and making an angle of approximately45 with the tensile direction. It can be expected that the shape of the fibers, which act as ”shear bandsbarriers”, is important in their capacity to inhibit shear bands. The shape of fibers is important at two levels.First it affects the arrangement of fibers in the unit cell. For instance, it can be favorable to clustering ofparticles, leaving large areas of inclusions-free matrix where plastic strain is likely to localize. At a smallerscale an elongated particle perpendicular to a shear band will form an effective barrier.Random microstructures were generated with three shapes of fibers : circular, elliptical (aspect ratio=3.333), equilateral triangles. The volume fraction was 47 . × .
9% higher). However the flow stress is significantly higher for the composite with triangular inclusions(5 .
2% higher). This ”hardening” effect can be attributed to the fact that at a given volume fraction trian-gles form more efficient barriers to shear band formation. This efficiency can be related to the length of theprojection of the fiber orthogonally to the shear bands. The minimal length, the maximal length, and theaverage length over all possible orientations are reported in Table 10 for each shape of fibers at a given area s . For circular fibers these three quantities are equal to the radius of the fiber (2 (cid:112) s/π ). .3 Penetrable fibers When the matrix is elastic ideally plastic, the overall response of the composite is strongly influenced bythe existence of continuous paths in the matrix, connected from one cell to the other. The contiguity of thematrix obviously plays a crucial role in the formation of these paths, which are ruled out when the matrixis not contiguous.In order to study this effect, different configurations at different volume fractions were generated with penetrable fibers. The centers of the fibers were first chosen at random. Then the volume fraction of thereinforcing phase was controlled by increasing the radius of the fibers (all fibers at a given volume fractionhad identical radius). The matrix was assumed to be elastic perfectly plastic. The results of the simulationscan be analyzed as follows :1. When the fiber volume fraction is small, shear bands can be passed through the matrix. Accordingto Drucker’s remark, the resulting overall flow stress of the composite coincides with the flow stressof the matrix under plane strains, 2 σ / √
3. However, when the fiber volume fraction is very small, anearly homogeneous deformation of the matrix is more favorable (less energy is dissipated in the plasticdeformation) and no strain concentration is observed. Then the overall flow stress of the compositestands between the flow stress of the matrix σ and the flow stress of the matrix under plane strains(2 σ / √ ± o on the tensile direction. Periodiccontinuous paths can again be passed through the matrix but they are tortuous. The bands wherethe plastic strain concentrates have a nonvanishing width. The stress-strain response of the compositeagain reaches a limit value, one higher than the flow stress of the matrix under plane strains. Theoverall flow stress increases with the volume fraction of fibers, and the increase is closely related to thetortuosity of the ”shear” bands.3. When the fibers percolate and form a contiguous phase, the matrix loses contiguity. No periodiccontinuous path can be passed through the matrix. This leads to a drastic modification of the stress-strain curve of the composite, which is no longer limited. The composite behaves asymptotically as anelastic plastic material with linear hardening . To illustrate the capability of the method to deal with complex microstructures, we have considered a realmicrostructure taken from the work of Bornert [12] (see also [21]). The materials studied in [21] were two-phase iron/silver blends, manufactured with powder metallurgy techniques. The digital image was obtainedby Scanning Electron Microscopy. The microstructure is shown in Figure 10 (a). Clearly meshing thismicrostructure for application of the FEM would be a considerable task. The present numerical method canhandle such a microstructure as easily as the simpler ones shown in previous examples. In the numericalsimulation each phase is considered elastic-plastic following a J -flow theory with isotropic hardening of theVon-Mises type. The stress/strain curves for each constituent under uniaxial tension are shown in Figure 10(c). The applied loading is uniaxial tension in the horizontal direction. The map of equivalent strain is shownin Figure 10 (b) at an overall strain E = 3 . A new numerical technique has been developed to investigate the local and overall response of nonlinearcomposites. The advantages of the method are the following :1. Images of microstructures can be directly used in the analysis, which avoids meshing the microstructure.Complex microstructures can be investigated. Part of the efficiency of the method is due to the use ofFFT packages.. The iterative procedure does not require the formation or inversion of a stiffness matrix.3. Convergence is fast.However the method has some limitations.1. Convergence is not ensured for materials containing voids or rigid inclusions.2. The number of degrees of freedom is high by comparison with the FEM (typically an image with1024 × Acknowledgements . Most computations were carried out at the Institut M´editerran´een de Technologie inMarseille; the funds being provided by the PACA region. The other computations were carried out at theInstitut du D´eveloppement et des Ressources en Informatique Scientifique funded by CNRS. The authorsare indebted to Michel Bornert for fruitful discussions and for providing the image of the microstructureshown in 10 (a).
Green’s operator of a linear elastic material
The auxiliary problem of a homogeneous material with stiffness c subject to a periodic polarization field τ plays an important role in the method which has been proposed. Its solution, which can be found in severaltextbooks ( e.g. Mura [23]), can be expressed in terms of the Fourier transform of the polarization field bymeans of the Fourier transform of the Green’s operator of the following systems of equations σ ( x ) = c : ε ( u ∗ ( x )) + τ ( x ) ∀ x ∈ V div σ ( x ) = ∀ x ∈ V, σ . n − , u ∗ (20)In Fourier space, these equations take the formˆ σ ij ( ξ ) = i c ijkh ξ h ˆ u ∗ k ( ξ ) + ˆ τ ij ( ξ ) , i ˆ σ ij ( ξ ) ξ j = 0 . (21)(It is hoped that the index i will not be confused with the complex number i = √− σ ij between the two equations in (21) yields K ik ( ξ ) .u ∗ k = ˆ τ ij ( ξ ) ξ j , where K ( ξ ) denotes the acoustic tensor of the homogeneous material, K ik ( ξ ) = c ijkh ξ h ξ j . Thenˆ u ∗ k ( ξ ) = i N ki ( ξ ) ˆ τ ij ( ξ ) ξ j = i2 ( N ki ( ξ ) ξ j + N kj ( ξ ) ξ i ) ˆ τ ij ( ξ ) , where the symmetry of τ has been used and where N ( ξ ) denotes the inverse of K ( ξ ). Thereforeˆ ε kh ( u ∗ ) = i2 ( ξ h ˆ u ∗ k ( ξ ) + ξ k ˆ u ∗ h ( ξ )) = ˆΓ khij ( ξ ) ˆ τ ij ( ξ ) , (22)with ˆΓ khij = 14 (cid:0) N hi ( ξ ) ξ j ξ k + N ki ( ξ ) ξ j ξ h + N hj ( ξ ) ξ i ξ k + N kj ( ξ ) ξ i ξ h (cid:1) , (23)and ˆ τ ij ( ξ ) = < τ ij ( x ) e − i ξ. x > . (24)The strain field induced at each point x of the unit cell V by an initial stress τ can be determined from (22),(23) and (24). These formulas give the explicit form of the operator Γ and of the operation ∗ considered insection 2: (cid:15) ( u ∗ ) = − Γ ∗ τ . Deatailed expressions of Γ can be found in Mura [23] for different types of anisotropy for the referencemedium. Its expression is particularly simple when the material is isotropic with Lam´e coefficients λ and µ ; the above expression becomes : c ijkh = λ δ ij δ kh + µ ( δ ik δ jh + δ ih δ jk ) .K ij ( ξ ) = (cid:0) λ + µ (cid:1) ξ i ξ j + µ | ξ | δ ij N ij ( ξ ) = 1 µ | ξ | (cid:18) δ ij − ξ i ξ j | ξ | λ + µ λ + 2 µ (cid:19) . Therefore: ˆΓ khij ( ξ ) = 14 µ | ξ | ( δ ki ξ h ξ j + δ hi ξ k ξ j + δ kj ξ h ξ i + δ hj ξ k ξ i ) − λ + µ µ ( λ + 2 µ ) ξ i ξ j ξ k ξ h | ξ | Imposing a macroscopic stress direction.
In the above described algorithm the overall strain is prescribed by assessing the value of the Fourier trans-form of the strain field at the zero frequency : ˆ ε ( ) = E . It is often convenient (or necessary) to impose the overall stress Σ , rather than the overall strain E . Atypical example is provided by uniaxial tension in the transverse direction as described by (16). In stronglynonlinear problems it is even necessary to impose only the direction of the overall stress and to drive theloading by means of an auxiliary parameter (arc length method). The algorithm can be modified to accountfor loadings in the form Σ = k S and E : S = t, (25)where S is the prescribed direction of overall stress (by direction of stress we refer to a direction in the6-dimensional space of stresses), k is the unknown level of overall stress and t , which serves as a loadingparameter, is the component of the overall strain in this direction. Then, the overall strain and stress E i et Σ i have to be determined by means of (25). For this purpose, at iterate i , σ i − and ε i − being known, theloading level t i being known but k i being unknown, E i and Σ i are subject to : Σ i − c : E i = < σ i − > − c : < ε i − > Σ i = k i S , E i : S = t i (26)Elimination of Σ i yields E i = k i c − : S − c − : < σ i − > + < ε i − > (27)and k i = t i + ( c − : < σ i − > − < ε i − > ) : S c − : S : S Therefore the modification brought into the algorithm (12) is an additional step to determine E i accordingto (27), which is then prescribed as the overall strain through :ˆ ε i ( ) = E i . It is worth noting that the condition < ε i > = E i is met at each step of the iterative procedure, whereasthe equality < σ i > = Σ i is met only at convergence. The difference arises from the fact that σ i is deducedfrom the constitutive law, whereas Σ i is deduced from (26). Indeed, once convergence is reached, one has E i = E i − = < ε i − >, and, according to (26), Σ i = < σ i − > = < σ i > . Radial return algorithm
The equations governing a plastic material obeying a J flow theory with isotropic hardening read :˙ σ = c : ( ˙ ε − ˙ ε p ) , ˙ ε p = 32 ˙ p s σ eq , (28)˙ p = 0 when σ eq − σ ( p ) < , ˙ p > σ eq − σ ( p ) = 0 . [ B (29) ε p is the plastic strain, p is the equivalent plastic strain ˙ p = (cid:0) ˙ ε p : ˙ ε p (cid:1) / . c is the stiffness tensor, assumedto be isotropic and characterized by a bulk modulus k and a shear modulus µ .Time is discretized into intervals [ t n , t n +1 ]. F n denotes the value of a function F at time t n . ε n , σ n and p n denote the strain, stress and equivalent plastic strain at time t n . Given the mechanical fields at step n , andgiven the strain field ε n +1 at step n + 1, the constitutive law amounts to finding the stress field σ n +1 andthe equivalent plastic strain field p n +1 . Replacing time differentiation by a finite difference in (28) provides σ n +1 − σ n = c : (cid:0) ε n +1 − ε n − ˙ ε pn +1 × ( t n +1 − t n ) (cid:1) . The elastic prediction is σ n +1 T = σ n + c : (cid:0) ε n +1 − ε n (cid:1) . (30)After due account of plastic incompressibility, (30) gives σ n +1 = σ n +1 T − t n +1 − t n ) µ ˙ ε pn +1 . Alternatively, making use of the flow rule (28) and of the decomposition of σ n +1 into a spherical stress anddeviator stress tr( σ n +1 ) = tr( σ n +1 T ) = tr( σ n ) + 3 k tr( ε n +1 − ε n ) (31) s n +1 = s n +1 T − t n +1 − t n ) µ ˙ p n +1 σ n +1 eq s n +1 (32)(31) can be re-written, assuming that there are no initial stresses or strains at time t ,tr( σ n +1 ) = 3 k tr( ε n +1 )The radial return method is based on the observation that, according to (32), the deviators s n +1 and s n +1 T are proportional. The Von Mises stresses associated with σ n +1 and σ n +1 T are therefore related through : σ n +1 eq = ( σ n +1 T ) eq − µ ( p n +1 − p n ) (33)- If ( σ n +1 T ) eq < σ ( p n ), the step is purely elastic, σ n +1 = σ n +1 T , p n +1 = p n . - If ( σ n +1 T ) eq ≥ σ ( p n ), the material plastifies at step n + 1, σ n +1 eq = σ ( p n +1 ) and (33) reduces to : σ ( p n +1 ) + 3 µ p n +1 = ( σ n +1 T ) eq + 3 µ p n Assuming that hardening is positive (no softening), the function h ( p ) = σ ( p ) + 3 µp can be invertedto give p n +1 = h − (cid:0) ( σ n +1 T ) eq + 3 µp n (cid:1) (34)The case of linear hardening leads to simple inversion. Indeed, in this case, σ ( p ) = σ + Hp and (34)reduces to p n +1 = 3 µH + 3 µ p n + ( σ n +1 T ) eq − σ H + 3 µ . The case of a perfectly plastic material, corresponding to H = 0, is covered by the above relation.When h − is not available in a closed form, it can be approximated by linear interpolation. When k ∈ [ h ( p l ) , h ( p l +1 )], p = h − ( k ) is approximated by p l + ( k − h ( p l )) p l +1 − p l h ( p l +1 ) − h ( p l ) .inally, the algorithm used in our computations reads : ε n , σ n , p n , ε n +1 being known , Compute s n +1 T = s n + 2 µ ( ε n +1 − ε n ) , ( σ n +1 T ) eq = (cid:0) s n +1 T : s n +1 T (cid:1) / Test If ( σ n +1 T ) eq < σ ( p n ) p n +1 = p n s n +1 = s n +1 T Else p n +1 = h − (cid:0) ( σ n +1 T ) eq + 3 µp n (cid:1) s n +1 = σ ( p n +1 )( σ n +1 T ) eq s n +1 T End of testUpdate tr( σ n +1 ) = tr( σ n ) + 3 k tr( ε n +1 − ε n ) σ n +1 = tr( σ n +1 ) IId + s n +1 (35) number of iterations E fi b e r / E m a t r i x Figure 1: Transverse Young’s modulus of the composite. Dependence of the number of iterations atconvergence on the contrast of the elastic moduli of phases ( e ≤ − ). Square array as shown in Figure 5.Spatial resolution 128 ×
128 pixels. Fiber volume fraction 47.5 %. Poisson coefficients ν f = ν m = 0 .
35. Thestiff phase is the fiber. number of pixels C P U ti m e i n s ec ond s ( on C r a y Y M P ) Figure 2: CPU time on one processor of a CRAY YMP as a function of the size N of the problem. The solidline is obtained by linear regression on all points and passes through the origin. a) x (b) Figure 3: Two-phase laminate. ( E = 68 . , ν = 0 . E = 400 GPa , ν = 0 .
23 ). Volumefraction of both phases 50%. Spatial resolution 32 ×
32 pixels. Applied loading: pure shear in the plane( x , x ). (a): map of the local strain field ε . (b) : cut through ε along an arbitrary horizontal line. a) (b) x (c) Figure 4: Circular fiber in a matrix with elastic mismatch between the phases ( E m = 68 . , ν m = 0 . E f = 400 GPa , ν f = 0 .
23 ). Shear loading: E = 0 . , E ij = 0 ∀ ( i, j ) (cid:54) = (1 , ε . (a): analytical solution, (b): numerical simulation. Spatial resolution 1024 × ε at x = 0. Dotted line: analytical solution, solid line: numerical simulation. a) (b) Figure 5: Standard fiber distributions. (a) : square lattice, the unit cell contains one fiber. (b): hexagonallattice, the unit cell contains 1 + 4 × = 2 fibers.Figure 6: Periodic unit cell containing 16 circular fibers randomly placed. .0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 E xx xx Random S q u a r e o S qu a r e o A v e r a g e H e x a g . o H e x a g . o M a t r i x Fibers (a) E xx xx Random Square 0 o Square 45 o AverageHexag. 45 o Hexag. 0 o MatrixFibers (b)
Figure 7: Overall stress-strain response computed with the present method. Volume fraction of fibers:47 . (resp: Square 45 ): fibers placed at the nodes of a square lattice, tension at 0 (resp. 45 ). Hexag. 0 (resp: Hexag. 45 ): fibers placed at the nodes of a hexagonal lattice, tension at 0 (resp. 45 ).a) (d)(b) (e)(c) (f)Figure 8: Two different microstructures ( (a) and (d) ) and the corresponding plastic strain maps. Thematrix is elastic - ideally plastic in (b) and (f). The matrix is elastic plastic with linear hardening in (c)and (f). Transverse uniaxial tension. Overall strain E = 1%. 0% strains are displayed in black, 10%strains (and more) are displayed in white. Straight slip bands can form easily in configuration (a). The slipbands are more tortuous in configuration (d). When the matrix is ideally plastic the overall flow stress ofconfiguration (d) is 6 .
4% higher than the flow stress in configuration (a). .0 0.002 0.004 0.006 0.008 0.01 E f = . % f = . % f = . % f = . % f = . % f=54.6% f=46.9% f=39.2%f=24.4% f=14.8% f=3.1% Figure 9: Overall stress-strain response of fiber reinforced materials at different fiber volume fraction f . 100 penetrable circular fibers with increasing radius in an elastic ideally plastic matrix. ( E m =68 . , ν m = 0 . σ = 68 . M P a and E f = 400 GPa , ν f = 0 .
23 ). Tension at 0 . a) (b) (log) ( M P a ) I r o n S i l v e r b l e n d m a t e r i a l (c) Figure 10: (a): Microstructure of a silver/iron blend material observed by Scanning Electron Microscopy.(b): response of the individual constituents under uniaxial tension. (c): numerical simulation. Uniaxialtension in the horizontal direction. Overall strain E = 3 . . Influence of spatial resolution on the overallYoung’s modulus. Square arrangement Hexagonal arrangementResolution Young’s modulus Error (%) Young’s modulus Error (%)32 129 670. 0.83 140 810. 0.8864 128 400. -0.16 140 200. 0.44128 128 750. 0.12 139 680. 0.07256 128 660. 0.05 139 520. -0.04512 128 600. 0.00 139 580. 0.00Table 2: Square and hexagonal array. Transverse tension at 45 ◦ . Influence of spatial resolution on theoverall Young’s modulus. Square arrangement Hexagonal arrangementResolution Flow stress Error (%) Flow stress Error (%)32 112.39 15.04 88.48 0.6064 107.46 9.99 88.32 0.42128 102.29 4.70 88.10 0.18256 99.65 2.00 88.01 0.07512 98.61 0.93 87.95 0.001024 98.01 0.32 * *2048 97.70 0.00 * *Table 3: Square and hexagonal array. Transverse tension at 0 ◦ . Influence of spatial resolution on the overallflow stress.quare arrangement Hexagonal arrangementResolution Flow stress Error (%) Flow stress Error (%)32 79.558 0.00 79.554 0.0064 79.558 0.00 79.554 0.00128 79.558 0.00 79.554 0.00256 79.558 0.00 79.554 0.00Table 4: Square and hexagonal array. Transverse tension at 45 ◦ . Influence of spatial resolution on theoverall flow stress. Square arrangement Hexagonal arrangementResolution Hardening modulus Error (%) Hardening modulus Error (%)32 14.4 10 ◦ . Influence of spatial resolution on the overallhardening modulus. Square arrangement Hexagonal arrangementResolution Hardening modulus Error (%) Hardening modulus Error (%)32 4.94 10 ◦ . Influence of spatial resolution on theoverall hardening modulus.umber of Number of Young’s modulus standard error onfibers tests mean (GPa) deviation (GPa) mean (%)4 100 143.7 3.9 0.279 50 143.4 3.1 0.3016 40 143.0 2.6 0.2936 25 143.1 1.51 0.2164 27 143.2 1.33 0.19256 10 142.9 0.57 0.13Table 7: Random configurations. Transverse uniaxial tension in the horizontal direction. Influence of thesize of the unit cell on the overall Young’s modulus.Number of Number of Flow stress standard Errorfibers tests mean (MPa) deviation (MPa) on mean (%)4 100 89.54 6.07 0.689 50 88.01 5.04 0.8116 40 87.94 4.99 0.9036 25 88.15 2.17 0.4964 27 88.70 2.07 0.51256 10 88.88 0.64 0.23Table 8: Random configurations. Transverse uniaxial tension in the horizontal direction. Influence of thesize of the unit cell on the overall flow stress.Fiber shape Young’s Modulus Flow stress Hardening modulusmean (MPa) mean (MPa) mean (Mpa)circle 142 260 86.9 9 382triangle 142 250 91.4 10 448ellipse 142 330 88.7 9 180Table 9: Random configurations. Transverse uniaxial tension in the horizontal direction. Effect of theshape of fibers on the effective properties of the composite.iber shape Maximal length Minimal length Average lengthcircle 1 1 1triangle 1.35 1.17 1.29ellipse 1.72 0.52 0.97Table 10: Projections along different angles of fibers with the same surface s = 1. References [1] H. Moulinec and P. Suquet, A fast numerical method for computing the linear and nonlinear propertiesof composites, C. R. Acad. Sc. Paris II 318 (1994) 1417–1423.[2] H. Moulinec and P. Suquet, A FFT-based numerical method for computing the mechanical properties ofcomposites from images of their microstructure, in: R. Pyrz, ed., Microstructure-Property Interactionsin Composite Materials (Kluwer Academic Pub., Dordrecht, 1995) 235–246.[3] D.F. Adams and D.R. Doner, Transverse normal loading of a unidirectional composite, J. Comp. Mat.1 (1967) 152-164.[4] T. Christman, A. Needleman, and S. Suresh. An experimental and numerical study of deformation inmetal-ceramic composites, Acta Metall. Mater. 37 (1989) 3029–3050.[5] V. Tvergaard, Analysis of tensile properties for a whisker–reinforced metal-matrix composite, ActaMetall. Mater. 38 (1990) 185–194.[6] J.C. Michel and P. Suquet, On the strength of composite materials: variational bounds and numericalaspects, in: C. Mota-Soares and M.P. Bendsoe, eds, Topology Design of Structures (Kluwer AcademicPub., Dordrecht, 1993) 355–374.[7] J.R. Brockenborough, S. Suresh, and H.A. Wienecke, Deformation of metal-matrix composites withcontinuous fibers: geometrical effects of fiber distribution and shape, Acta Metall. Mater. 39 (1991)735–752.[8] H.J. B¨ohm, F.G. Rammerstoffer, and E. Weissenbeck, Some simple models for micromechanical inves-tigations of fiber arrangements in MMCs, Comput. Mat. Sc. 1 (1993) 177–194.[9] T. Nakamura and S. Suresh, Effects of thermal residual stresses and fiber packing on deformation ofmetal-matrix composites, Acta Metall. Mater. 41 (1993) 1665–1681.[10] Ch. Dietrich, M.H. Poech, H.F. Fischmeister, and S. Schmauder, Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study, Comp.Mater. Sc.1