Characterization of random stress fields obtained from polycrystalline aggregate calculations using multi-scale stochastic finite elements

Abstract

The spatial variability of stress fields resulting from polycrystalline aggregate calculations involving random grain geometry and crystal orientations is investigated. A periodogram-based method is proposed to identify the properties of homogeneous Gaussian random fields (power spectral density and related covariance structure). Based on a set of finite element polycrystalline aggregate calculations the properties of the maximal principal stress field are identified. Two cases are considered, using either a fixed or random grain geometry. The stability of the method w.r.t the number of samples and the load level (up to 3.5% macroscopic deformation) is investigated.

Keywords

polycrystalline aggregates / crystal plasticity / random fields / spatial variability / correlation structure

Introduction

In pressurized water reactors of nuclear plants, the pressure vessel constitutes one element of the second safety barrier between the radioactive fuel rods and the external environment. It is made of 16MND5 (A508) steel which is forged and welded. In case of operating accidents such as LOCA (loss of coolant accident), the pressure vessel is subjected to a pressurized thermal shock due to fast injection of cold water into the primary circuit. If some defects (e.g., cracks) were present in the vessel wall this may lead to crack initiation and propagation and to the brittle fracture of the vessel. The detailed study of the embrittlement of 16MND5 steel under irradiation is thus a great concern for electrical companies such as EDF.

The brittle fracture behavior of the 16MND5 steel has been thoroughly studied in the last decade using the local approach of fracture theory [ 1] and the so-called Beremin model [ 2], which assumes that cleavage is controlled by the propagation of the weakest link among a population of pre-existing micro-defects in the material. This approach has been recently coupled with polycrystalline aggregates simulations [ 3, 4].

The main idea is to model a material representative volume element (RVE) as a polycrystalline synthetic aggregate and compute the stress field under given load conditions. As a post-processing a statistical distribution of defects (carbides) is sampled over the volume. In each Gauss point of the finite element mesh the cleavage criterion is attained somewhere along the load path if a) the equivalent plastic strain has attained some threshold (cleavage initiation) and b) a Griffith-like criterion applied to the size of the carbide in this Gauss point is reached (cleavage propagation). Within the weakest link theory the failure of a single critical carbide induces the failure of the RVE.

From a single RVE simulation (i.e., a single stress field) various distributions of carbides are drawn, each realization leading to a maximal principal stress associated to failure. Then the distribution of these quantities is fitted using a Weibull law [ 4]. In such an approach, the current practice of computational micromechanics assumes that the RVE is large enough to represent the behavior of the material so that a single polycrystalline analysis is carried out (the large CPU required by polycrystalline simulations also favors the use of a single simulation). However it is believed that numerous parameters such as grain geometry and orientation may influence the stress field and thus the final result.

The connection between micromechanics and stochastic methods has been given much attention in the past few years, as shown in Graham-Brady et al. [ 5]; Stefanou [ 6]; Xu and Chen [ 7]. Many papers are devoted to developing probabilistic models for reproducing a random microstructure, e.g., Graham and Baxter [ 8]; Liu et al. [ 9]; Chung et al. [ 10]; Chakraborty and Rahman [ 11]. The specific representation of polycrystalline microstructures has been addressed in Arwade and Grigoriu [ 12]; Grigoriu [ 13]; Li et al. [ 14]; Kouchmeshky and Zabaras [ 15] among others. The propagation of the uncertainty on the microstructure through a micromechanical model in order to study the variability of the resulting strain and stresses has not been addressed much though (see e.g., Kouchmeshky and Zabaras [ 16]).

In this paper it is proposed to identify the properties of a stress random field resulting from the progressive loading of a polycrystalline aggregate. More precisely, assuming that the stress random field is Gaussian, a procedure called periodogram method is devised, which allows one to identify the correlation structure of the resulting stress field.

The paper is organized as follows: in Section 2 basics of Gaussian random fields are recalled and the periodogram method is presented (Dang et al. [ 17]). The polycrystalline aggregate computational model is detailed in Section 3. The methodology for identifying the correlation structure of the resulting stress field is presented in Section 4. Two application cases are then investigated in Sections 5 and 6, namely an aggregate with fixed grain boundaries and random crystallographic orientations (Section 5) and an aggregate with both random geometry and orientations (Section 6). The variance of the resulting stress field as well as the spatial covariance function and its correlation lengths is investigated in details. The properties of the identified random fields will be used in a forthcoming study in the context of the local approach to fracture, as explained above.

Inference of the properties of a Gaussian random field

In this section an identification method called periodogram is presented, which uses an estimator of the Power Spectral Density (PSD) in order to identify the correlation structure of a Gaussian homogeneous random field. Based on original developments by Stoica and Moses [ 18] and Li [ 19] for unidimensional fields, it has been extended to two-dimensional cases by Dang et al. [ 17]. As it relies upon the use of the Fast Fourier Transform (FFT) its computational efficiency is remarkable.

Definitions

A Gaussian random field Z(x, ω) is completely defined by its mean value µ(x), its standard deviation σ(x) and its auto-covariance function C(x, x'). It is said homogeneous if the mean value µ(x) and the standard deviation σ(x) are constant in the domain of definition of x and the auto-covariance function C(x, x') only depends on the shift h = x - x'. Let us introduce the n-th statistical moment {\mathit{m}}_z^n and the spatial average {\mathit{m}}_{\mathit{V}}^n :

(1) m_Z^n=E\left[Z^n\left(x_{\mathrm{0}},\omega \right)\right]={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}^n\left(x_{\mathrm{0}},\omega \right)f_z\left(z,x_{\mathrm{0}}\right)\mathrm{d}z},

(2) m_{\mathit{V}}^n=\underset{\mathit{V}\to \infty }{\lim }{\displaystyle \frac{\mathrm{1}}{\mathit{V}}}{\displaystyle \underset{\mathit{V}}{\int }{z}^n\left(x,{\omega }_{\mathrm{0}}\right)\mathrm{d}x}\mathrm{.}

The field is said ergodic if its ensemble statistics is equal to the spatial average, i.e., m_Z^n=m_{\mathit{V}}^n [ 20]. Several popular covariance models for two-dimensional homogeneous random fields are presented in Table 1. In this table, σ is the constant standard deviation of the field, h₁, h₂ are the components of the shift h in the two directions, l₁, l₂ are the correlation lengths in the two directions. Gaussian and exponential models are plotted in Fig. 1 for the sake of illustration. Note that we call correlation length the parameters that appear in the definition of the covariance functions. This is not to be confused with the scale of fluctuation [ 21], which combine both the shape of the covariance function and the lengths l₁, l₂. In one dimension, denoting by ρ(x; l) the autocorrelation function, the scale of fluctuation may be defined by:

\mathrm{2}{l}_{\mathrm{c}}={\displaystyle {\int }_{-\infty }^{\infty }\rho \left(x\mathrm{;}l\right)\mathrm{d}x},

which reduces to l_c = l for the exponential correlation function and l_{\mathrm{c}}=\sqrt{\mathrm{\pi }}/\mathrm{2}\approx \mathrm{0.886}l for the Gaussian case. Similar expression are available in two and three dimensions, see e.g., Xu and Chen [ 7].

The power spectral density (PSD) of the random field is the Fourier transform of its covariance function as a result of the Wiener-Khintchine relationship [ 22]. The following relationships hold:

(3) S\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}C\left(h_{\mathrm{1}},h_{\mathrm{2}}\right)e^{-i\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{h}_{\mathrm{1}}}{e}^{-i\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{h}_{\mathrm{2}}}d{h}_{\mathrm{1}}d{h}_{\mathrm{2}}}},

(4) C\left(h_1,h_2\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}S\left(f_1,f_2\right)e^{i2\mathrm{\pi }{f}_1{h}_1}{e}^{i2\mathrm{\pi }{f}_2{h}_2}d{f}_{1}d{f}_2}}.

The PSD of the Gaussian and exponential covariance models are presented in Table 1.

Identification of a PSD

Empirical periodogram

One considers an ergodic homogeneous random field Z(x₁, x₂), (x₁, x₂) ∈ D₁× D₂ ⊂ {ℝ}^{\mathrm{2}} for which a single realization z(x₁, x₂) is available. If the random field was defined over an infinite domain, the classical estimation of the covariance function would be:

(5) \widehat{C}\left(h_{\mathrm{1}},h_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}MN}}{\displaystyle \sum _{n=-N}^{N-\mathrm{1}}{\displaystyle \sum _{m=-M}^{M-\mathrm{1}}Z\left(x_{\mathrm{1}n}+h_{\mathrm{1}},x_{\mathrm{2}m}+h_{\mathrm{2}}\right)Z\left(x_{\mathrm{1}n},x_{\mathrm{2}m}\right)}}.

By definition, the Fourier transform of the covariance estimation is an estimation of the PSD.

(6) \widehat{S}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}MN}}\tilde{Z}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)\tilde{Z}*\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}MN}}{\left|\tilde{Z}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)\right|}^{\mathrm{2}},

where |.| denotes the modulus operator.

In practice, the problem is to estimate the periodogram from a limited amount of data gathered on N × M grid {z(x_1i, x_2j), i = 1, …, N; j = 1, …, M}. Due to symmetry, the covariance estimation in Eq. (5) is recast as follows:

(7) \widehat{C}\left(h_{1}k,h_{2l}\right)=\frac{1}{NM}{\displaystyle \sum _{n=1}^{N-k}{\displaystyle \sum _{m=1}^{M-l}Z\left(x_{1n}+h_{1k},x_{2m}+h_{2l}\right)Z\left(x_{1n},x_{2m}\right)}}.

By taking the expectation of the above equation one gets:

(8) E[\widehat{C}\left(h_{\mathrm{1}k},h_{\mathrm{2}l}\right)]={\displaystyle \frac{N-k}{N}}{\displaystyle \frac{M-l}{M}}E\left[Z\left(x_{\mathrm{1}n}+h_{\mathrm{1}k},x_{\mathrm{2}m}+h_{\mathrm{2}l}\right)Z\left(x_{\mathrm{1}n},x_{\mathrm{2}m}\right)\right]={\displaystyle \frac{N-k}{N}}{\displaystyle \frac{M-l}{M}}C\left(h_{\mathrm{1}},h_{\mathrm{2}}\right).

The latter equation exhibits some bias term between the expectation of the estimator and the covariance function C(h₁, h₂). Using the symmetry of the covariance function, one can write:

(9) E\left[\widehat{C}\left(h_{\mathrm{1}k},h_{\mathrm{2}l}\right)\right]=w_B\left(k,l\right)C\left(h_{\mathrm{1}},h_{\mathrm{2}}\right),

where w_B(k, l) is the triangle window, also known as the Bartlett window (Fig. 2):

(10) w_B\left(k,l\right)=\{\begin{array}{cc}{\displaystyle \frac{N-\left|k\right|}{N}}{\displaystyle \frac{M-\left|l\right|}{M}},\hfill & \mathrm{i}\mathrm{f}\left|\mathit{k}\right|\le N\mathrm{;}\left|l\right|\le M\hfill \\ \mathrm{0},\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}.

Consequently the expectation of the periodogram estimation becomes:

(11) E\left[\widehat{S}\left(f_{\mathrm{1}k},f_{\mathrm{2}l}\right)\right]=F\left\{E\left[\widehat{C}\left(h_{\mathrm{1}k},h_{\mathrm{2}l}\right)\right]\right\}=F\left\{w_B\left(k,l\right)C\left(h_{\mathrm{1}},h_{\mathrm{2}}\right)\right\}=W_B(f_{\mathrm{1}},f_{\mathrm{2}})\ast S\left(f_{\mathrm{1}},f_{\mathrm{2}}\right),

where F and W_B(f₁, f₂) respectively denote the 2D Fourier transform operator and the Fourier transform of the Bartlett window and * denotes the convolution product. This window tends to a Dirac pulse when N, M tend to infinity and w_B tends to a unit constant. Thus the periodogram estimation is asymptotically unbiased. However it is not consistent since its variance does not tend to zero [ 22]. Furthermore using this window leads to a convolution product which introduces additional computational burden. Hence in practice, the modified periodogram presented in the next section is used to estimate the PSD of the random field.

Modified periodogram

The modified periodogram is built up in order to avoid the convolution product with the transformed window W_B(f₁, f₂) in Eq. (11). In this approach, the data is multiplied directly with the window w(x, y) before the Fourier transform is carried out. It aims at filtering the data to limit the influence of long distance terms and to focus on the information given by the short distance terms. This leads to the following estimate:

(12) \widehat{S}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{NMU}}{\left|F\left\{z\left(x_{\mathrm{1}},x_{\mathrm{2}}\right).w\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)\right\}\right|}^{\mathrm{2}},

where U is the energy of the window calculated by:

(13) U={\displaystyle \frac{\mathrm{1}}{D_{\mathrm{1}}{D}_{\mathrm{2}}}}{\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=\mathrm{1}}^M{w}^{\mathrm{2}}\left(x_{\mathrm{1}i},x_{\mathrm{2}j}\right)}},

and D₁, D₂ denote the size of the two-dimensional domain D₁×D₂. Various window functions are proposed in Preumont [ 22], see Table 2. In this paper we will use mainly the Blackman window (Fig. 2).

Average modified periodogram

As shown in Section 2.2.1, the estimation of the periodogram is asymptotically unbiased, however not consistent since its variance does not tend to zero when N, M tend to infinity. The averaging of the modified periodogram will solve this problem. Assume that L realizations of the random field are available. For each realization z^l(x₁; x₂), one calculates the periodogram as in Eq. (12):

(14) {\widehat{S}}^l\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{NMU}}{\left|F\left\{z^l\left(x_{\mathrm{1}},x_{\mathrm{2}}\right).w\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)\right\}\right|}^{\mathrm{2}},

with 1≤l≤L. Then one calculates the average periodogram by:

(15) \overline{S}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{L}}{\displaystyle \sum _{l=\mathrm{1}}^L{\widehat{S}}^l\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)}.

Therefore the variance of the average periodogram is:

(16) \mathrm{V}\mathrm{a}\mathrm{r}\left[\overline{S}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)\right]={\displaystyle \frac{\mathrm{1}}{L}}\mathrm{V}\mathrm{a}\mathrm{r}\left[\widehat{S}\left(f_{\mathrm{1}},f_{\mathrm{2}}\right)\right].

It is then obvious that this variance tends to zero when L tends to infinity, making the “average modified periodogram” approach more robust.

Final algorithm for PSD estimation

As a summary, the algorithm to estimate the PSD of a random field from L realizations may be decomposed into the four following steps:

1) multiplication of each realization by a selected window, e.g., the Blackman window (see Table 2);

2) computation of 2D Fourier transform of the product of the current realization by the filtering window;

3) computation of the modulus of the result to obtain the PSD estimation of each realization;

4) averaging of the L PSD estimations.

Once the empirical periodogram has been computed, a theoretical periodogram is selected (e.g., Gaussian, exponential, etc., see Table 1) and fitted using a least-square procedure [ 23]. In case of multiple potential forms for the theoretical periodogram the best fitting is selected according to the smallest residual.

Modeling polycrystalline aggregates

In this section the computational mechanical model used in this study is presented. It simulates a tensile test on a bidimensional polycrystalline aggregate under plane strain conditions. The various ingredients are discussed, namely:

– the microstructure of the material and its synthetic representation;

– the material constitutive law;

– the boundary conditions applied onto the aggregate;

– the mesh used in the finite element simulation.

Material characterization

The material is a 16MND5 ferritic steel with a granular microstructure. The ferrite has a body centered cubic (BCC) structure. Three families of slip system should be taken into account, namely {110} 〈111〉, {112} 〈111〉, {123} 〈111〉. However, following Franciosi [ 24] it is assumed that the glides on the plane 123 are a succession of micro-glides on the planes 110, 112. This leads to consider only the two first families, which yields 24 slip systems by symmetry.

Crystal plasticity

The model for crystal plasticity chosen in this work has been originally formulated in Meric and Cailletaud [ 25] within the small strain framework. The total strain rate {\dot{ϵ}}_{ij} is classically decomposed as the sum of the elastic strain rate {\dot{ϵ}}_{ij}^e and plastic strain rate {\dot{ϵ}}_{ij}^P :

(17) {\dot{ϵ}}_{ij}={\dot{ϵ}}_{ij}^e+{\dot{ϵ}}_{ij}^P.

The elastic part follows the Hooke’s law and the plastic part is calculated from the shear strain rates of the 24 active slip systems:

(18) {\dot{ϵ}}_{ij}^P={\displaystyle \sum _{g=\mathrm{1}}^{\mathrm{24}}{\dot{\gamma }}^g{R}_{ij}^g},

where {\dot{\mathrm{\gamma }}}^g is the shear strain rate of the slip system g and R_{ij}^g is the Schmid factor which presents the geometrical projection tensor. The latter is calculated from the normal vector to the gliding plane n and the direction of gliding m.

(19) R_{ij}^g={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{m}}_{\mathit{i}}{\mathit{n}}_j+m_j{n}_i\right).

The Resolved Shear Stress (RSS) τ^g of the slip system g is the projection of the stress tensor via the Schmid factor.

(20) {\tau }^g=R_{ij}^g{\sigma }_{ij}.

The shear strain rates {\dot{\gamma }}^g of each slip system g are the internal variable that describes plasticity. The evolution of these variables depends on the difference between the RSS {\tau }^g and the actual critical RSS {\tau }_c^g in an elastoviscoplastic setting:

(21) {\dot{\gamma }}^g={\left({\displaystyle \frac{{\tau }^g-{\tau }_c^g}{K}}\right)}^{\mathit{n}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\tau }^g\right),

where K and n are material constants, and sign(a) = a/|a| if a ≠ 0 and 0 otherwise. Note that this formula corresponds to an elastoviscoplastic constitutive law but the viscous effect will be negligible if sufficiently large values of K and n are selected. Its power form allows one to automatically detect the active slip systems. All the systems are considered active but the slip rate is significant only if the RSS τ^g is much higher than the critical RSS {\tau }_c^g . This procedure allows one to numerically smooth the elastoplastic constitutive law.

The critical RSS {\tau }_c^g evolves according to the following isotropic hardening law:

(22) {\tau }_c^g={\tau }_{c\mathrm{0}}^g+Q^g{\displaystyle \sum _{s=\mathrm{1}}^{\mathrm{24}}{h}^{gs}\left(\mathrm{1}-e^{-b^g{\gamma }_{cum}^s}\right)},

where {\gamma }_{cum}^s={\int }_{t_0}^t{\displaystyle \left|{\dot{\gamma }}^S\right|dt} . The exponential term presents the hardening saturation in the material when the accumulated slip is high. {\tau }_{c\mathrm{0}}^g is the initial critical RSS on the considered system g. Q^g and b^g are parameters which depend on the material. h^gs is the hardening matrix of size 24 × 24 whose component h^gs presents the hardening effect of the system g on the system s. In the present work, one considers only two families of slip systems named 110〈111〉, 112〈111〉. Thus the hardening matrix h^gs is completely defined by four coefficients h₁, h₂, h₃, h₄ only. The values of these coefficients and this matrix are presented in Mathieu [ 26]. All the parameters describing crystal plasticity for 16MND5 steel are gathered in Table 3.

Microstructure and boundary conditions

The construction of the aggregate is based on the Voronoi polyhedra model [ 27], generated in this work with the Quickhull algorithm [ 28]. The geometry of the resulting synthetic aggregate, which is a simplified representation of the real microstructure of the 16MND5 steel, is shown in Fig. 3. It corresponds to a square of size 1,000 (this is a relative length which shall be scaled with a real length depending on the grain size). Grain boundaries are considered as perfect interfaces. Note that more detailed models of grains have been proposed recently using so-called Laguerre tessellations [ 29] in order to better fit the observed distributions of grain size, see e.g., Zhang et al. [ 30]; Leonardia et al. [ 31].

The same crystallographic orientation, defined by the three Euler angles {\phi }_{\mathrm{1}} , \varphi , {\phi }_{\mathrm{2}} , is randomly assigned to all integration points inside each individual grain using a uniform distribution. In Fig. 3(a), the color of each grain corresponds to a given crystallographic orientation. The mesh is generated by the BLSURF algorithm [ 32] of the Salome software (http://www.salome-platform.org). The mesh of the generated specimen is presented in Fig. 3(b). The finite elements are quadratic 6-node triangles with three integration points.

The boundary conditions applied onto the aggregate are sketched in Fig. 4. The lower surface is blocked along the Y direction. The displacements DX = DY = 0 are blocked at the origin of the coordinate system (lower left corner). On the upper surface, an homogeneous displacement is applied by steps in the Y direction up to a macroscopic strain equal to 3.5%. The computation is carried out using the open source finite element software Code_Aster (http://www.code-aster.org).

The computational cost for such a non linear analysis is high. The number of degrees of freedom of the finite element model is 33, 885. A parallel computing method based on sub-domain decomposition is used. One simulation of a full tensile test up to 3.5% strain requires about 2h computation time when distributed over four processors.

Results

In this section, we present the result of the simulation of a tensile test on the 2D aggregate at different scales. We define the mean stress and strain tensor calculated in a volume V by:

(23) \mathbf{\Sigma }={\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \underset{V}{\int }\sigma \mathrm{d}V},

(24) \mathit{E}={\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \underset{V}{\int }ϵ\mathrm{d}V}.

Figure 5 shows the macroscopic strain/stress curve. It is observed that Σ_XX = 0 as expected whereas the uniaxial behavior shows a global elastoplastic behavior.

At the mesoscopic scale one can observe the mean strain-stress relationship in each grain as shown in Fig. 6. Because of the different crystallographic orientations in each grain, the mean elastoplastic behavior is different from grain to grain. Furthermore, whereas the mean stress Σ_XX calculated in over the whole specimen is zero, the mean values calculated in each single grain are scattered around zero. This observation shows the first scale of heterogeneity of the material.

The microscopic behavior of a single grain (Grain #24, see tag in Fig. 3) is finally studied. The mean behavior and the strain-stress relationship at each node of this grain are plotted in Fig. 7 for four levels of macroscopic strain, namely E_YY = 0.15%, 0.65%, 1.5%, and 3.5%.

In this figure the blue point represents the stress field within the grain for a macroscopic strain level E_YY = 0.15%. This single point shows that the stress field is homogeneous within the grain in the elastic domain. The red points represent the stress values σ_YY in each node of the grain at E_YY = 0.65% macroscopic strain. One observes that the mean strain calculated for this single grain is 0.85% and the maximal strain value ϵ_YY in a specific node may attain about 2%. Similar effects are observed at other levels of macroscopic strain, which show the heterogeneity of the strain and stress fields at the microscopic scale. It is observed that the scattering around the mean curve increases with the macroscopic strain. Indeed, for the final loading step corresponding to E_YY= 3.5% the mean strain in the grain is about 4.54%, while the local strain varies form 2.4 to 9%.

Identification of the maximal principal stress field

In this section the method developed in Section 2 is applied to the identification of the properties of the random stress field in polycristalline aggregate calculations. More specifically the maximal principal stress field σ_I that is computed from repeated polycrystalline simulations is considered. Throughout the paper this stress field is considered Gaussian. This is a strong assumption which shall be considered as a first approximation. Indeed the maximal principal stress is positive in nature under the uniaxial loading that is considered and a Gaussian assumption cannot totally fit this feature. Yet it is believed that the results obtained in terms of the description of the spatial variability (covariance functions), which is the main outcome of the paper, will not be strongly influenced by this assumption. Note that methods for identifying the properties of non Gaussian random fields have been recently developed, see e.g., Perrin et al. [ 33].

Finite element calculations and projection

The maximal principal stress field is assumed to be Gaussian and homogeneous (the latter assumption will be empirically checked as shown in the sequel). The periodogram method is applied using K = 35 realizations of stress fields, i.e., 35 full elastoplastic analysis of aggregates up to a macroscopic strain of 3.5%. The identification is carried out successively at various levels of the macroscopic strain. Two cases are considered:

– Case #1: the grains geometry is the same for all the finite element calculations. Only the crystallographic orientations are varying from one calculation to the other.

– Case #2: both the grains geometry and the crystallographic orientations vary.

The input data of the identification problem is the maximal principal stress field σ_I obtained from the finite element calculations. As the periodogram method is based on a regular sampling of the random field under consideration, the brute result (i.e., the maximal principal stress at the nodes of the mesh) has to be projected onto a regular grid. This operation is carried out using internal routines of Code Aster. Note that a slice of width 100 (i.e., 10% of total size) is discarded along the edges of the aggregate in order to avoid the effect of boundary conditions on the computed stress field, as suggested in Mathieu [ 26]. A typical maximal principal stress field is shown in Fig. 8.

Check of the homogeneity of the field

As it was described in Section 2 the periodogram method assumes that the random field under consideration is homogeneous. From the available realizations SIG_i(x, y), i = 1,…,35 one first checks empirically this assumption using the following approach:

1) The ensemble mean and variance of the field is computed point-by-point throughout the grid for an increasing number of realizations K = 2,…,35:

(25) {\mu }_K\left(x,y\right)={\displaystyle \frac{\mathrm{1}}{K}}{\displaystyle \sum _{i=\mathrm{1}}^{K}SI{G}_i\left(x,y\right)},

(26) {\sigma }_K^{\mathrm{2}}\left(x,y\right)={\displaystyle \frac{\mathrm{1}}{K-\mathrm{1}}}{\displaystyle \sum _{i=\mathrm{1}}^K{\left(SI{G}_i\left(x,y\right)-\mu \left(x,y\right)\right)}^{\mathrm{2}}},

If the field is homogeneous these quantities should tend to constant values that are independent from the position (x, y) when K tends to infinity.

2) To measure the magnitude of the spatial fluctuation of the latter, the spatial average and spatial variance of a realization of a field Z(x, y) sampled onto a N × M grid is defined by:

(27) {\overline{\mu }}_Z={\displaystyle \frac{\mathrm{1}}{NM}}{\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=\mathrm{1}}^{M}Z\left(x_i,y_j\right)}},

(28) {\overline{\sigma }}_Z^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{NM}}{\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=\mathrm{1}}^M{\left(Z\left(x_i,y_j\right)-{\overline{\mu }}_Z\right)}^{\mathrm{2}}}},

whereas the associated “spatial” coefficient of variation is defined by:

(29) C{V}_Z={\displaystyle \frac{{\overline{\sigma }}_Z}{{\overline{\mu }}_Z}}.

3) The spatial coefficient of variation of the ensemble mean and variance (Eqs. (25) and (26)) are computed and plotted as a function of K. If the underlying random field is homogeneous it is expected that the curves of CV_µK and C{V}_{{\sigma }_K^{\mathrm{2}}} converge to zero.

Choice of theoretical periodograms and fitting

From a visual inspection of the obtained empirical periodograms it appears that a Gaussian or an exponential model of periodogram such as those presented in Table 1 may be consistent with the data. However it appeared in the various analyses that the peak of the periodogram is not always in zero. An initial frequency is thus introduced which shifts the theoretical periodogram. Finally, due to lack of fitting of the simple periodograms (e.g., Gaussian and exponential), a combination thereof is also fitted. The most general model finally reads:

(30) S\left(f_x,f_y\right)=\mathrm{}{\sigma }_{\mathrm{1}}^{\mathrm{2}}\mathrm{\pi }{l}_{x\mathrm{1}}\exp \left[{\mathrm{\pi }}^{\mathrm{2}}{l}_{x\mathrm{1}}^{\mathrm{2}}{\left(f_x-f_{x\mathrm{0}}^{\left(\mathrm{1}\right)}\right)}^{\mathrm{2}}\right]+{\sigma }_{\mathrm{2}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}{l}_{x\mathrm{2}}}{\mathrm{1}+\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}{l}_{x\mathrm{2}}^{\mathrm{2}}{\left(f_x-f_{x\mathrm{0}}^{\left(\mathrm{2}\right)}\right)}^{\mathrm{2}}}}\mathrm{}{\displaystyle \frac{\mathrm{2}{l}_{y\mathrm{2}}}{\mathrm{1}+\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}{l}_{y\mathrm{2}}^{\mathrm{2}}{\left(f_y-f_{y\mathrm{0}}^{\left(\mathrm{2}\right)}\right)}^{\mathrm{2}}}},

where l_x1 , l_y1 , l_x2, l_y2 are correlation lengths in each direction X and Y (aniso-tropic field) for each component (1) (Gaussian part) and (2) (exponential part). Similarly f_{x\mathrm{0}}^{\left(\mathrm{1}\right)} , f_{y\mathrm{0}}^{\left(\mathrm{1}\right)} , f_{x\mathrm{0}}^{\left(\mathrm{2}\right)} , f_{y\mathrm{0}}^{\left(\mathrm{2}\right)} are initial shift frequencies.

Note that Eq. (30) corresponds only to positive values of f_x, f_y. The periodogram is then extended by symmetry for negative frequencies. In terms of associated covariance models, the linear combination of periodograms leads to a linear combination of covariance models. The initial frequency shift in the periodogram leads to oscillatory cosine terms in the covariance by inverse Fourier transform:

(31) C\left(h_x,h_y\right)={\sigma }_{\mathrm{1}}^{\mathrm{2}}\exp \left[-\left({\displaystyle \frac{h_x^{\mathrm{2}}}{l_{x\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{h_y^{\mathrm{2}}}{l_{y\mathrm{1}}^{\mathrm{2}}}}\right)\right]\cos \left(\mathrm{2}\mathrm{\pi }{f}_{x\mathrm{0}}^{\left(\mathrm{1}\right)}{h}_x\right)\cos \left(\mathrm{2}\mathrm{\pi }{f}_{y\mathrm{0}}^{\left(\mathrm{1}\right)}{h}_y\right)+{\sigma }_{\mathrm{2}}^{\mathrm{2}}\exp \left[-\left({\displaystyle \frac{\left|h_x\right|}{l_{x\mathrm{2}}}}+{\displaystyle \frac{\left|h_y\right|}{l_{y\mathrm{2}}}}\right)\right]\cos \left(\mathrm{2}\mathrm{\pi }{f}_{x\mathrm{0}}^{\left(\mathrm{2}\right)}{h}_x\right)\cos \left(\mathrm{2}\mathrm{\pi }{f}_{y\mathrm{0}}^{\left(\mathrm{2}\right)}{h}_y\right).

To compare the various fittings the least-square residual between the empirical periodogram \overline{S}\left(f_x,f_y\right) (Eq. (15)) and the fitted periodogram S^theor(f_x, f_y) is finally computed. The following non dimensional error estimate is used:

(32) ϵ=\sqrt{{\displaystyle \frac{\mathrm{1}}{NM}}{\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=\mathrm{1}}^M{[\overline{S}\left(f_{xi},f_{yj}\right)-S^{theor}\left(f_{xi},f_{yj}\right)]}^{\mathrm{2}}}}}/\underset{\left(f_x,f_y\right)}{\max }\overline{S}\left(f_x,f_y\right).

Results – Case #1: fixed grain geometry

Check of the homogeneity

First the homogeneity of the maximal principal stress field is checked using the methodology proposed in Section 4.2. Figure 9 shows the evolution of CV_µK and C{V}_{{\sigma }_K^{\mathrm{2}}} . These quantities regularly decrease and it is seen that they would tend to zero if a larger number of realizations was available. This leads to accepting the assumption that the random field is homogeneous since the fluctuations around the constant spatial average tend to zero when K increases.

Identification of periodograms at 3.5% macroscopic strain

The average empirical periodogram obtained from L = 35 realizations of the maximal principal stress field σ_I at 3.5% of macroscopic strain is plotted in Fig. 10(a).

Table 4 presents the results of the fitting of the average empirical periodogram calculated from 35 realizations of the field using three models, namely Gaussian, exponential and a mixed “Gaussian+ exponential” as in Eq. (30)

From the results in Table 4 it appears that the mixed model provides a significantly smaller least-square error than that obtained from the Gaussian and exponential models respectively. The corresponding fitted periodogram is plotted in Fig. 10(b).

To better appreciate the quality of the fitting, two-dimensional cuts of the empirical (resp. fitted) periodogram are given in Figs. 11–13. Figure 11 corresponds to a cut along the X direction for two values of f_y = 0; 0.0013. Figure 12 corresponds to a cut along the Y direction for two values of f_x = 0; 0.0013. Finally Fig. 13 corresponds to a cut along the diagonal f_x = f_y.

From the above figures it appears that the fitting of the empirical periodogram by a mixed model is remarkably accurate. It is interesting to interpret the fitted parameters reported in Table 4. First it is observed that the amplitude of each component of the mixed periodogram is similar since σ₁≈σ₂. The variance of the field is equal to {\sigma }_{\mathrm{1}}^{\mathrm{2}}+{\sigma }_{\mathrm{2}}^{\mathrm{2}}\approx \mathrm{6},\mathrm{309}\mathrm{M}\mathrm{P}{\mathrm{a}}^{} . The associated standard deviation is 79.4 MPa. As the mean principal stress is 720 MPa at 3.5% macroscopic strain, the coefficient of variation of the field is about 11%.

To interpret the correlation length parameters let us define the mean size of a grain S_g for such a two-dimensional aggregate. As the volume of edge length equal to 1,000 corresponds to 100 grains, the equivalent diameter of a single grain reads:

(33) D_g=\sqrt{{\displaystyle \frac{\mathrm{4}}{\mathrm{\pi }}}{S}_g}=\sqrt{{\displaystyle \frac{\mathrm{4}}{\mathrm{\pi }}}{\displaystyle \frac{\mathrm{1000}\times \mathrm{1000}}{\mathrm{100}}}}=\mathrm{112.8.}

Thus the correlation lengths obtained from the fitting vary from 0.55 to 1.3D_g. This shows that the characteristic dimension of the underlying microstructure (i.e., D_g) is of the same order of magnitude as these parameters. In other words the scale of local fluctuation of the stress field is related to the grain size, as heuristically expected. Moreover, it appears that the lengths in the X and Y directions are almost identical. The stress field does not show any significant anisotropy in this case.

Influence of the number of realizations

In this section the stability of the fitted parameters as a function of the number of available realizations K used in the average periodogram method is considered. In practice the procedure applied in the previous paragraph is run using K = 8, 9, …,35 realizations of the stress field. The evolution of the standard deviations (σ₁, σ₂) is shown in Fig. 14. The evolution of the correlation lengths l_(x,y)(1,2) is shown in Fig. 15. The evolution of the initial frequencies f_{{\left(x,y\right)}_{\mathrm{0}}}^{\left(\mathrm{1},\mathrm{2}\right)} is shown in Fig. 16.

From these figures it appears that the fitted parameters tend to a converged value when at least 20 realizations of the stress field are used for their estimation.

Influence of the macroscopic strain level

In this section the evolution of the parameters of the fitted periodograms as a function of the macroscopic strain is investigated. For this purpose the methodology presented in Section 5.2 is applied using the realizations of the maximal principal stress fields corresponding to various levels of the loading curve, i.e., various values of the equivalent macroscopic strain E_YY = 0, …, 3.5%.

The evolution of the standard deviations (σ₁, σ₂) is shown in Fig. 17. The two components of the periodogram (e.g., Gaussian and exponential) contribute for approximately the same proportion to the total variance of the field since the curves are almost superimposed. Note that these standard deviations increase with the applied load in the same way as the mean load curve (Fig. 5).

The evolution of the correlation lengths l_(x,y)(1,2) is shown in Fig. 18. The evolution of the initial frequencies f_{{\left(x,y\right)}_{\mathrm{0}}}^{\left(\mathrm{1},\mathrm{2}\right)} is shown in Fig. 19. It is observed that once plasticity is settled (i.e., once the macroscopic strain EYY is greater than ~ 0.5%) the parameters describing the fluctuations of the maximal principal stress field are almost constant. This conclusion is valid for both the correlation lengths and the initial frequencies. Note that the convergence is faster for the parameters related to the X direction, i.e., the direction that is transverse to the one-dimensional loading. Finally it is also observed that f_{y\mathrm{0}}^{\left(\mathrm{1}\right)} is almost equal to zero whatever the load level, thus the zero value in Table 4.

Results – Case #2: random grain geometry

In this section both the randomness in the grain geometry and in the crystallographic orientations are taken into account. A total number of 35 finite element models are run. In each case, the grain geometry is obtained from a uniform sampling of points from which a Voronoï tessellation is built.

Check of the homogeneity

As in Section 5 the homogeneity of the maximal principal stress field is checked using the methodology proposed in Section 4.2. Figure 20 shows the evolution of CV_µK and C{V}_{{\sigma }_K^{\mathrm{2}}} . These quantities regularly decrease and it is seen that they would tend to zero if a larger number of realizations was available. This leads to accepting the assumption that the random field is homogeneous.

Identification of periodograms at 3.5% macroscopic strain

The average empirical periodogram obtained from L = 35 realizations of the maximal principal stress field σ_I at 3.5% of macroscopic strain is plotted in Fig. 22(a). Three types of theoretical periodograms have been fitted as in the previous section, which lead to the conclusion that the mixed model that combines a Gaussian and an exponential component is best suited. The fitted parameters are gathered in Table 5 where the parameters fitted for Case #1 are also recalled for the sake of comparison.

To check the accuracy of the fitting, two-dimensional cuts of the empirical (resp. fitted periodogram) are plotted in Fig. 22 (cut along the X direction), Fig. 23 (cut along the Y direction), Fig. 24 (cut along the diagonal f_x = f_y). Again the fitting is remarkably accurate, meaning that the fitted model of periodogram accurately represents the spatial variability of the maximal principal stress field.

It can be observed from the figures in Table 5 that the fitting is of equal quality in both cases (relative error less than 2×10^-3). As far as the contribution of each component of the periodogram is concerned, the symmetry reported in Case #1 is not existing anymore since the standard deviation of the exponential contribution (σ₂ = 81.6) is much greater than that of the Gaussian part (σ₁ = 35.8). The total variance of the field is 7,940 MPa², corresponding to a standard deviation of 89.1 MPa and a coefficient of variation of 12%. Thus there is a little more scattering in the random stress field obtained in Case #2 when considering both the random grain geometry and orientations.

The correlation lengths associated with the exponential part do not differ much in Case #2 compared to Case #1 (corresponding here to 1.5 to 2.4 D_g). In contrast the correlation lengths related to the Gaussian part are increased, which tends to produce less rapidly varying realizations. This may be explained by the fact that the grain boundaries are “averaged” in Case #2 whereas they were fixed in Case #1. The stress concentrations that are usually observed at the grain boundaries are thus smoothed in Case #2 compared to Case #1.

Influence of the number of realizations

In this section one considers the stability of the fitted parameters as a function of the number of available realizations K used in the average periodogram method. The procedure applied in the previous paragraph is run using K = 8, 9,…,35 realizations of the stress field. The evolution of the standard deviations (σ₁, σ₂) and the initial frequencies f_{\left(x\right)\mathrm{0}}^{\left(\mathrm{1},\mathrm{2}\right)} is shown in Fig. 25 (note that f_{\left(y\right)\mathrm{0}}^{\left(\mathrm{1},\mathrm{2}\right)} = 0 in the present case). The evolution of the correlation lengths l_(x,y)(1,2) is shown in Fig. 26.

From these figures it clearly appears that the fitted parameters are almost constant when the number of realizations of the stress field used in their estimation increases. The minimal number of K = 8 could be used here without significant errors although it is recommended to keep a value of K = 20 as in Case #1 for robustness.

Influence of the macroscopic strain level

Finally the evolution of the parameters of the fitted periodograms as a function of the macroscopic strain E_YY is investigated. For this purpose the identification method is applied using the realizations of the maximal principal stress fields corresponding to various levels of the loading curve, i.e., various values of the equivalent macroscopic strain E_YY = 0, …,3.5%.

The evolution of the two standard deviations look similar to the results obtained in Case #1 (Fig. 27). It is observed that the ratio σ₂/σ₁ is almost constant all along the loading path up to 3.5% strain. As far as the initial frequencies are concerned, there is a complete independance with the load level as soon as E_YY is greater than ~ 0.5%, i.e., when plasticity has settled in the aggregate. The same conclusion can be drawn for the various correlation lengths.

Conclusions

The distribution of stresses in a material at a microscopic scale (where heterogeneities such as grain structures are taken into account) has been given much attention in the context of computational homogenization methods. However the current methods usually stick to a deterministic formulation. Starting from the premise that any representative volume element (such as a polycrystalline aggregate) is a single specific realization of a random quantity, the present paper aims at using methods of computational stochastic mechanics for representing the (random) stress field.

After recalling the basic mathematics of Gaussian random fields, the paper presents a periodogram method for estimating the parameters describing the spatial fluctuation of a random field from a collection of realizations of this field. This method is adapted in two dimensions from well-known techniques originating from signal processing.

The material under consideration, namely the 16MND5 steel used in nuclear pressure vessels is then presented together with a local modeling by polycrystalline finite element calculations. From a collection of 35 realizations of the (maximal principal) stress field, the spatial correlation structure of the latter is identified. By fitting various theoretical periodograms, a mixed model combining a Gaussian and an exponential contribution is retained. These two contributions may be empirically interpreted as follows: The Gaussian part corresponds to the fluctuation from grain to grain; the (less smooth) exponential component corresponds to the sharp grain boundaries stress concentrations.

Two cases are considered, namely a “fixed-geometry” case in which only the crystallographic orientations changes within the 35 realizations (fixed grain boundaries), and a “variable geometry” in which the grain geometry is randomly sampled for each realization. In both cases, a good convergence of the procedure is observed when the number of realizations increases. A set of 20 realizations is recommended, although good results are already obtained for ~8 realizations in Case #2.

Moreover it is shown that the correlation lengths are of the same order of magnitude as the grain size. The initial frequencies that are required for a best fitting of the periodogram and that translate into some kind of spatial periodicity in the covariogram could be explained by spurious edge effects due to the limited size of the aggregate. This should be investigated more in details in further analysis.

Another important result is drawn from the comparison of the fitted parameters at various load levels. Once plasticity is settled within the aggregate, the parameters describing the spatial fluctuations of the field are almost constant. Moreover the variance of the field (sum of the variance of each component of the periodogram) increases proportionally to the mean strain/stress curve, meaning that the coefficient of variation of the stress field is almost constant (around 11% for the fixed geometry and 12% for the variable geometry).

The results presented in this paper should be confirmed by additional investigations under different types of loading (e.g., biaxial loading). The tools that are presented here may be applicable to three-dimensional aggregates and stress fields at a much larger computational cost though. This work is currently in progress.

The identified stress fields may eventually be re-simulated: new realizations of the stress fields are straightforwardly obtained at a low computational cost by random field simulation techniques such as the spectral approach or the circulant embedding method [ 17, 22, 34]. This allow us to apply local approach to fracture analysis (such as that presented in Mathieu et al. [ 26]) for the assessment of the brittle fracture of metallic materials, as shown in Dang et al. [ 35].

References

[1] Tanguy B. Modélisation de l’essai Charpy par l’approche locale de la rupture – application au cas de l’acier 16MND5 dans le domaine de la transition. Ph.D thesis, École des Mines Paris Tech, 2001

[2] Beremin F, Pineau A, Mudry F, Devaux J C, D’Escatha Y, Ledermann P. A local criterion for cleavage fracture of a nuclear pressure vessel steel. Metallurgical Transactions, 1983, 14A(11): 2277–2287

[3] Mathieu J, Berveiller S, Inal K, Diard O. Microscale modelling of cleavage fracture at low temperatures: influence of heterogeneities at the granular scale. Fatigue & Fracture of Engineering Materials & Structures, 2006, 29(9—10): 725–737

[4] Mathieu J, Inal K, Berveiller S, Diard O. A micromechanical interpretation of the temperature dependence of Beremin model parameters for French RPV steel. Journal of Nuclear Materials, 2010, 406(1): 97–112

[5] Graham-Brady L, Arwade S, Corr D, Gutiérrez M, Breysse D, Grigoriu M, Zabaras N. Probability and materials: from nano- to macro-scale – a summary. Probabilistic Engineering Mechanics, 2006, 21(3): 193–199

[6] Stefanou G. The stochastic finite element method: past, present and future. Computer Methods in Applied Mechanics and Engineering, 2009, 198(9-12): 1031–1051

[7] Xu X F, Chen X. Stochastic homogenization of random elastic multi-phase composites and size quantification of representative volume element. Mechanics of Materials, 2009, 41(2): 174–186

[8] Graham L, Baxter S. Simulation of local material properties based on moving-window GMC. Probabilistic Engineering Mechanics, 2001, 16(4): 295–305

[9] Liu H, Arwade S, Igusa T. Random composites characterization using a classifier model. Journal of Engineering Mechanics, 2007, 133(2): 129–140

[10] Chung D B, Gutiérrez M A, Borst R. Object-oriented stochastic finite element analysis of fibre metal laminates. Computer Methods in Applied Mechanics and Engineering, 2005, 194(12—16): 1427–1446

[11] Chakraborty A, Rahman S. Stochastic multiscale models for fracture analysis of functionally graded materials. Engineering Fracture Mechanics, 2008, 75(8): 2062–2086

[12] Arwade S, Grigoriu M. Probabilistic model for polycrystalline microstructures with application to intergranular fracture. Journal of Engineering Mechanics, 2004, 130(9): 997–1005

[13] Grigoriu M. Nearest neighbor probabilistic model for aluminum polycrystals. Journal of Engineering Mechanics, 2010, 136(7): 821–829

[14] Li Z, Wen B, Zabaras N. Computing mechanical response variability of polycrystalline microstructures through dimensionality reduction techniques. Computational Materials Science, 2010, 49(3): 568–581

[15] Kouchmeshky B, Zabaras N. Microstructure model reduction and uncertainty quantification in multiscale deformation processes. Computational Materials Science, 2010, 48(2): 213–227

[16] Kouchmeshky B, Zabaras N. The effect of multiple sources of uncertainty on the convex hull of material properties of polycrystals. Computational Materials Science, 2009, 47(2): 342–352

[17] Dang H X, Sudret B, Berveiller B. Benchmark of random fields simulation methods and links with identification methods. In: Faber M, Köhler J, Nishijima K. eds. Proceedings of the 11th International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP11), Zurich, Switzerland, 2011

[18] Stoica P, Moses R. Introduction to spectral analysis. Prentice Hall, 1997

[19] Li J. Spectral estimation. Dep. of Electrical and Computer Engineering – University of Florida, 2005

[20] Cramer H, Leadbetter M. Stationary and Related Processes. Wiley & Sons, 1967

[21] Vanmarcke E. Random Fields: Analysis and Synthesis. Cambridge, Massachussets: The MIT Press, 1983

[22] Preumont A. Vibrations aléatoires et analyse spectrale. Presses polytechniques et universitaires romandes, 1990

[23] Marquardt D. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 1963, 11(2): 431–441

[24] Franciosi P. The concepts of latent hardening and strain hardening in metallic single crystals. Acta Metallurgica, 1985, 33(9): 1601–1612

[25] Meric L, Cailletaud G. Single crystal modelling for structural calculation: Part 2- F.E. Implementation. Journal of Engineering Materials and Technology, 1991, 113(1): 171–182

[26] Mathieu J. Analyse et modélisation micromécanique du comportement et de la rupture fragile de l’acier 16MND5: prise en compte des hétérogénéités microstructurales. Ph.D thesis, École Nationale Supérieure des Arts et Métiers, 2006

[27] Gilbert E. Random subdivisions of space into crystals. Annals of Mathematical Statistics, 1962, 33(3): 958–972

[28] Barber C, Dobkin D, Huhdanpaa H. The Quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software, 1996, 22(4): 469–483

[29] Lautensack C, Zuyev S. Random Laguerre tessellations. Advances in Applied Probability, 2008, 40(3): 630–650

[30] Zhang P, Balint D, Lin J. Controlled Poisson Voronoi tessellation for virtual grain structure generation: a statistical evaluation. Philosophical Magazine, 2011, 91(36): 4555–4573

[31] Leonardia A, Scardia P, Leonia M. Realistic nano-polycrystalline microstructures: beyond the classical Voronoi tessellation. Philosophical Magazine, 2012, 92(8): 986–1005

[32] Laug P, Borouchaki H. BLSURF-mesh generator for composite parametric surface. Tech. Rep. 0235, INRIA, User’s Manual, 1999

[33] Perrin G, Soize C, Duhamel D, Funfschilling C. A posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations. SIAM/ASA J. Uncertainty Quantification, 2014, 2(1): 745–762

[34] Dang H, Sudret B, Berveiller M, Zeghadi A. Identification of random stress fields from the simulation of polycrystalline aggregates. In: Proceedings of the 1st International Conference Computational Modeling of Fracture and Failure of Materials and Structures (CFRAC). Barcelona, Spain, 2011

[35] Dang H X, Berveiller M, Zeghadi A, Sudret B, Yalamas T. Introducing stress random fields of polycrystalline aggregates into the local approach to fracture. In: Deodatis G. ed. Proceedings of the 11th International Conference Structural Safety and Reliability (ICOSSAR’2013). New York, USA, 2013

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg