A stochastic multiscale approach for numerical modeling of complex materials – Application to uniaxial cyclic response of concrete
Abstract.
In complex materials, numerous intertwined phenomena underlie the overall response at macroscale. These phenomena can pertain to different engineering fields (mechanical, chemical, electrical), occur at different scales, can appear as uncertain, and are nonlinear. Interacting with complex materials thus calls for developing nonlinear computational approaches where multiscale techniques that grasp key phenomena at the relevant scale need to be mingled with stochastic methods accounting for uncertainties. In this chapter, we develop such a computational approach for modeling the mechanical response of a representative volume of concrete in uniaxial cyclic loading. A mesoscale is defined such that it represents an equivalent heterogeneous medium: nonlinear local response is modeled in the framework of Thermodynamics with Internal Variables; spatial variability of the local response is represented by correlated random vector fields generated with the Spectral Representation Method. Macroscale response is recovered through standard homogenization procedure from Micromechanics and shows salient features of the uniaxial cyclic response of concrete that are not explicitly modeled at mesoscale.
Pierre Jehel^{1,2}^{*}^{*}*Corresponding author: pierre.jehel[at]centralesupelec.fr
MSSMat, CNRS, CentraleSupélec, Université ParisSaclay, Grande voie des Vignes, 92290 ChâtenayMalabry, France
Department of Civil Engineering and Engineering Mechanics, Columbia University, 630 SW Mudd, 500 West 120th Street, New York, NY, 10027, USA
Manuscript form of the book chapter published online (Feb. 15, 2016):
http://link.springer.com/book/10.1007/9783319279961
©Springer International Publishing Switzerland 2016
A. Ibrahimbegovic (ed.), Computational Methods for Solids and Fluids
Computational Methods in Applied Sciences 41, chap. 6, pp. 123–160.
DOI: 10.1007/9783319279961_6
1. Introduction
Widelyused materials in engineering practice such as polymer, composite, steel, concrete, are characterized by engineering parameters for design purposes, while these latter homogeneous macroscopic mechanical properties actually result from heterogeneous structures at lower scales. Material can be qualified as complex as their macroscopic behavior result from numerous multiscale intertwined phenomena that have nonlinear and uncertain evolution throughout loading history. Modifications in the underlying structures of this category of materials can result in dramatic changes in mechanical behavior at the relevant macroscopic scale for engineering applications. Microcracks coalescence in the constitutive material of a structure challenges its capacity for meeting the performance level targeted during its design process. Alkaliaggregate reaction in concrete microscopic structure can lead to hazardous loss of bearing capacity in reinforced concrete structures. Accounting for phenomena at lower scales to reliably predict macroscopic response of heterogeneous structures is one of the challenges numerical multiscale simulation techniques have been developed for over the past ([1, 18, 8, 15] among many others).
In continuum mechanics, explicitly accounting for relevant mechanisms and structures in heterogeneous scales underlying macroscopic scale provides the rationale for representing characteristic features of homogenized material behavior laws at macroscale that can then be used for engineering design. Micromechanics has been developed to extract macroscopic local continuum properties from microscopically heterogeneous media through the concept of Representative Volume Element (RVE). An RVE for a material point at macroscale is statistically representative of the microscopic structure in a neighborhood of this point [23]. Also, Thermodynamics with Internal Variables provides a robust framework for modeling material response at macroscale according to a set of internal variables that carry information about the history of evolution mechanisms at lower scales without explicitly representing them [9, 21]. Other strategies to derive macroscopic mechanical properties of heterogeneous materials have been developed based on the introduction of a mesoscale, that is a scale that bridges the micro and macroscales. In [24], heterogeneities are represented by random fields introduced at a mesoscale, which defines socalled Statistical Volume Elements that tends to become RVEs as mesoscale grows; effective properties at macroscale are retrieved according to two hierarchies of scaledependent bounds obtained from either homogenous displacements or homogenous tensions applied on the boundary of the mesoscale. In [2], a mesoscale is explicitly constructed for representing the macroscopic behavior of heterogeneous quasibrittle materials. This mesoscale consists of a 3D finite element mesh composed of truss elements cut by inclusions. Truss element kinematics is enriched to account for discontinuities in the strain field due to the presence of inclusions along truss elements as well as discontinuities in the displacement field to account for possible cracks in the matrix, in the inclusions, or at their interface. With the improvement of computational ressources, stochastic homogenization of random heterogeneous media can now be achieved without introducing a mesoscale. In [5], an efficient numerical strategy is presented to obtain effective tensors of random materials by coupling random microstructures to tentative effective models within the Arlequin framework for model superposition [1]. In [30], microstructures composed of a medium with randomly distributed inclusions of random shapes are generated and their behaviors are simulated with the extended finite element method (XFEM); homogenized properties at macroscale are then derived through the computation of mean response using Monte Carlo simulations.
In the work presented in this chapter, we focus on the numerical representation of the homogenized onedimensional response of a concrete specimen in cyclic compressive loading, as it can be observed in lab tests. Figure 1 illustrates the main features of such an homogenized response: a backbone curve (dashed line) that is a nonlinear strain hardening phase () followed by a strain softening phase where strength degradation is observed; unloadingreloading cycles show that stiffness decreases while loading increases, hysteresis loops are generated. This typical response is observed at macroscale and results from numerous underlying mechanisms of physical or chemical nature at many different scales. For designing concrete structures, an equivalent homogeneous concrete model is sought, which has to represent concrete mechanical behavior in different loading conditions while accounting for mechanisms at lower scales [37]. Heterogeneities can be observed in concrete at different scales: aggregates of different sizes are distributed in a cement paste; the socalled interfacial transition zone where the aggregates are bound to the cement paste plays a key role in the concrete mechanical properties [38]; cement paste is composed of water, voids and of the products of the complete or partial hydration of the clinker particles, which generates a microscopic structure composed of numerous intertwined phases.
This chapter presents the basic ingredients of a stochastic multiscale approach developed to represent the macroscopic compressive cyclic response of a concrete specimen while attempting not to sacrifice too much of the complexity of this material. To that aim, two scales are considered: the macroscale where an equivalent homogenous concrete model capable of representing the main features that are shown in figure 1 is retrieved, and a mesoscale where heterogeneous local nonlinear response is assumed. Local response at mesoscale is modeled in the framework of Thermodynamics with Internal Variables and is seen as the homogenized response of mechanisms that occur at the micro or nano underlying scales. Spatial variability at mesoscale is introduced using stochastic vector fields. Homogenized macroscopic response is recovered using standard averaging method from micromechanics.
The chapter is organized as follow. In the next section, the ingredients of the proposed stochastic multiscale modeling are presented. First, the averaging method for computing the homogenized model response at macroscale is recalled. Then, the model of the mechanical local behavior of a material point at mesoscale is constructed. Finally, the Spectral Representation Method for generating stochastic vector fields that model heterogeneity at mesoscale is presented. In a third section, the numerical implementation of the approach in the framework of the finite element method is detailed. Before the conclusion, numerical applications are presented to demonstrate the capability of the proposed approach i) for yield homogeneous material behavior at macroscale without stochastic homogenization and ii) for representing salient features of macroscopic 1D concrete response in uniaxial cyclic compressive loading.
2. Multiscale stochastic approach for modeling concrete
Figure 2 presents the three following concepts, which further developments are based on:

Actual heterogeneous medium (Amesoscale): Concrete is made of aggregates distributed in a cement paste. Aggregates, cement and interface between both of them exhibit different mechanical responses. In the cement paste, micro and nanostructures also exist.

Equivalent heterogeneous medium (Emesoscale): The proposed approach does not consist in explicitly generating a multiphase medium with random distribution of aggregates of random geometry in a cement paste with known mechanical behavior for each phase. The approach followed here consists in generating a random medium at each point of which the mechanical response obeys a prescribed behavior that has uncertain parameters and that is the homogenized response of mixtures of aggregates and cement where mechanisms at lower scales are also involved but not explicitly modeled.

Equivalent homogeneous medium (macroscale): Homogenization of Emesoscale yields homogenized homogeneous concrete response. It will be shown in the numerical applications that one realization only of the random Emesoscale can be sufficient to retrieve homogenized properties at macroscale.
2.1. Homogenized material behavior at macroscale
We consider a material elementary domain (ED) that occupies a spatial domain . The boundary of the ED has outward normal , tension can be imposed on the part of the boundary while displacement can be imposed on , where and . There are no external forces other than applied on the ED and no dynamic effects are considered either. Then, the displacement vector field , and the strain and stress tensor fields, and , satisfy at any pseudotime :
(1)  
is the symmetric part of the gradient tensor , the superscript denoting the transpose operation. In the set of equations above, small strains are assumed and behavior law can be nonlinear.
We classically assume that any macroscopic quantity is connected to its Emesoscopic counterpart through domain averaging over the ED:
(2) 
is the measure of the spatial domain occupied by the ED centered at material point of the macroscale, and denotes a material point of the Emesoscale.
In all what follows, we will assume linear displacements imposed all over the boundary of :
(3) 
Hence, and . With this assumption, it can be shown (see e.g. [23, Chap. 1] or [39]) that:
(4) 
and also, because it is assumed there is no external forces applied on :
(5) 
where are the tension forces developed over .
Note that other boundary conditions could be considered. In any case, it is in general not possible to derive the strain or stress fields at Emesoscale from the macroscopic quantities, and consequently simplifying assumptions as in (3) are made. Whether it is displacements or forces that are imposed on , and whether these latter conditions are linear or periodic, this can influence the homogenized macroscopic response of the ED. However, this is out of the scope of this work where the consequences of assuming linear displacements imposed on will not be discussed.
With the boundary conditions (3) applied on , Hill’s lemma can be proved:
(6) 
which means that the medium recovered at macroscale through homogenization is energetically equivalent to the heterogeneous medium considered at Emesoscale.
The possibly nonlinear material response at mesoscale is expressed as:
(7) 
where is the tangent modulus at Emesoscale and the superimposed dot denotes partial derivative with respect to pseudotime. Thanks to Hill’s lemma, we then have the following two equivalent definitions for the tangent modulus at the homogenized macroscale:
(8) 
2.2. Material behavior law at Emesoscale
We assume a coupled damageplasticity model to be suitable for representing material response at any material point at Emesoscale (see figure 2). This choice is motivated by the fact that concrete Amesoscale is composed of both a ductile cement matrix that can be represented by a plastic model, and brittle aggregates that are confined in the cement paste and whose compressive response can be more realistically represented by a damage model. Hereafter, we develop a model in a way that allows for explicitly controlling the coupling of damage and plasticity. Indeed, as illustrated in figure 3, the response of material points at mesoscale can be either better represented by a damage model alone, or a plasticity model alone, or by the appropriate coupling of both models. This is developed in the framework of thermodynamics with internal variables [9, 21] where the internal variables carry the history of irreversible mechanisms occurring in the material at lower (micro and nano) scales.
2.2.1. Basic ingredients
The three basic ingredients for developing this model of local behavior law at Emesoscale are as follows:

Total deformation is split into damage () and plastic () parts:
(9) 
Stored energy function is defined as:
(10) with the fourthorder damage compliance tensor. and are the internal variables that drive the evolution of the material. Also, denoting by the elasticity tensor, we set initially, as the material is undamaged, . The elements of are parameters of the model.

A criterium function is introduced as:
(11) It defines the limit states between the states where there is no evolution of the internal variables () and those where there is evolution (). The socalled yield stress is a scalar parameter.
2.2.2. Material dissipation and state equation
Then, the material dissipation reads:
(12)  
should be nonnegative to comply with the principle of thermodynamics. In case there is no evolution of the internal variables, that is for loading steps that do not generate any change of state in the material, there is no evolution of the internal variables: and the process is assumed to be nondissipative, that is is null. According to equation (12), it then comes the state equation:
(13) 
Equation (13) is to this damage model what the more classical constitutive relation is to linear elasticity model.
Introducing this latter state equation into equation (12), we can rewrite:
(14) 
from where we define and .
2.2.3. Evolution of the internal variables
Following what has been done to derive the equations of mechanical models with plasticity solely [12], the evolution of the internal variables is obtained appealing to the principle of maximum dissipation. Accordingly, among all the admissible stresses, that is such that , it is those that maximize the material dissipation that have to be retained. This can be cast into a minimization problem with constraint [19]. Lagrange multiplier method can be used to solve it with the socalled Lagrangian reading:
(15)  
Here, we have split the total Lagrange multiplier so that with two Lagrange multipliers defined as and where is to be taken in the range . is a damageplasticity coupling parameter: if , and there is plasticity evolution only; if , only damage evolves in the material; and for any other inbetween, there is coupled evolution of both damage and plasticity.
In turn, the Lagrangian is also split into damage and plasticity parts:
(16) 
Both parts have to be minimized to ensure the total Lagrangian is minimum. The KuhnTucker optimality conditions associated to these minimization problems result in:
(17) 
Setting , this leads to the following equations of evolution of the internal variables:
(18)  
(19) 
Besides, this minimizing problem also yields the following socalled loading/unloading conditions:
(20) 
2.2.4. Damage and plasticity multipliers
In the case or , there is damage or plasticity evolution and, according to (20), as to remain null during the process so that the stresses remain admissible. We thus have the consistency condition that can be rewritten as:
(21) 
Remarking from (13) that and using equations (9), (18) and (19), we have:
(22) 
Then, assuming , the consistency condition (21) is satisfied when if:
(23) 
Or, with the damage and plasticity multipliers and :
(24) 
2.2.5. Tangent modulus
The tangent modulus at mesoscale is a fourthorder tensor that has been defined in (7) such that . Assuming that , the inverse of , exists (, where is the identity fourthorder tensor), and reminding equations (22) and (23), we have:
(25) 
To sum up, the proposed material model at Emesoscale is based on a set of internal variables that consists of the damage compliance tensor and the plastic deformation tensor . Besides, the model is parameterized by the elasticity tensor , the stress threshold above which damage or plastic evolution occurs, and the damageplasticity coupling coefficient : if , there is no plastic evolution and the material can only damage, while if , there is no damage evolution and the material is perfectly plastic. Figure 4 shows material constitutive behavior at two different material points and of the Emesoscale where the parameters take different values: parameters, and consequently local response, vary over the domain due to heterogeneities at Emesoscale.
2.3. Stochastic modeling of heterogeneous Emesoscale
Spatial variability at Emesoscale of a set of parameters over is conveyed through stochastic modeling: it is assumed that the fluctuations of correlated stochastic fields can describe the actual material heterogeneous mesostructure (Amesoscale). Thus, we introduce the probability space where is the sample space containing all the possible outcomes from the observation of the random phenomenon that is studied; is the algebra associated with ; is a probability measure. A real parameter taking values in is then considered as the realization of a random variable : . A random variable can be completely defined by its cumulative distribution function: or, when a probability density function (PDF) exists: .
Before we go on with the definition of stochastic fields, we recall some basic definitions for random variables. The mean and the variance of a random variable are defined as:
(26)  
(27) 
where is the socalled mathematical expectation and , . Also, and being two random variables, the covariance is defined as:
(28) 
where is the joint PDF of and , with . We also introduce the correlation, which is defined as:
(29) 
And finally the following correlation coefficient will also be used later on:
(30) 
2.3.1. Random vector fields for modelind heterogeneous mesostructure
It is assumed that the heterogeneity of the parameters , and of the model developed above for representing material response at Emesoscale over a concrete elementary domain (ED) can be represented as the realization of a random vector field. A random vector is a vector of random variables. Let be a spatial domain of dimension ; this can be a volume (), an area () or a length (). A random vector field over is a collection of random vectors indexed by the position . For any fixed , any component of , , is a random variable. In the case of random vector fields, we have the following definitions of the mean, the autocorrelation and crosscorrelation functions respectively:
(31)  
(32)  
(33) 
where is the separation distance vector between two points of .
To fully characterize random vector fields, we need the marginal and joint PDFs of all possible combinations of random variables . From a practical point of view, this implies that many concrete ED have to be considered, that for each of them the parameters of interest have to be identified at many points all over , so that these PDFs can be empirically constructed. If gathering such a huge amount of information was needed, the usefulness of using random vector field for modeling heterogeneity in concrete structure at mesoscale would be questionable. Therefore, we will make assumptions on the structure of the random vector field that would justify the efficiency of the proposed approach.
Emesoscale construction will rely on the two following assumptions:

Random fields will be generated as Gaussian. This means that random field is fully characterized by both its mean function and autocorrelation function . Nevertheless, nonGaussian random field can then be obtained through nonlinear translation of Gaussian field, which will be discussed in section 2.3.3.

Random fields are jointly homogeneous. This means that their mean function is independent of the position and that auto and crosscorrelation functions depend on the separation distance only:
(34) (35) Note that efficient techniques can be used to account for heterogeneity in the random field (see e.g. [25]).
Also, we will consider the possible ergodicity of the generated random fields in mean and correlation functions. One realization of such an ergodic random vector field contains all the statistical information needed to retrieve the first two moments: means and correlation functions can be computed as spatial averages.
2.3.2. Spectral representation of homogeneous Gaussian random vector fields
We present here the Spectral Representation Method for generating standard homogenous Gaussian fields [33, 34, 6, 7, 26]. Note that considering the Gaussian random fields to be standard, that is with zero mean and unit variance, does not introduce any loss of generality because nonstandard Gaussian fields can always be retrieved through linear transformation.
The basic ingredient is the definition of a target correlation matrix, which can be built from experimental observations for instance:
(36) 
Superscript has been added to highlight that these functions are target correlations, which should be retrieved in the statistical analysis of the generated random fields.
According to WienerKhinchin theorem, power spectral density functions , , and crossspectral density functions , , are the Fourier transform of the corresponding correlation functions:
(37)  
(38) 
where is the wave number vector, is the scalar product of the two vectors and , and is the imaginary unit. Power spectral density functions are by definition real functions of while crossspectral functions can be complex functions of . It has been shown in [32] that matrix is Hermitian and semidefinite positive, which implies that it can be decomposed using Cholesky’s method as:
(39) 
where denotes the complex conjugate and is a lower triangular matrix. The diagonal elements are real and nonnegative functions of while the offdiagonal elements can be complex:
(40) 
Then, the ^{th} component of a realization of a VD homogeneous standard Gaussian stochastic vector field with crossspectral density matrix reads:
(41) 
for and ,…, . In equation (2.3.2), the following notation has been introduced:
(42) 
Wave numbers increments are defined as:
(43) 
where ’s are socalled cutoff wave numbers such that can be assumed to be null for any . Also, are the different vectors composed of ’s and ’s where for all . For instance, for , there are the following different such arrangements: , , and , so that , , and . are independent sequences of independent random phase angles drawn at any wave number from a uniform distribution in the range .
Random fields generated with relation (2.3.2) are periodic along the axes, , with period:
(44) 
Also, the values of the field are bounded according to:
(45) 
It has been shown that the random fields generated according to equation (2.3.2) have the following properties [33, 34, 6, 7, 26]:

They tend to be standard Gaussian as , ; rate of convergence is investigated in [33].

They ensemble auto and crosscorrelations are identical to the target functions.

Each realization is ergodic in mean and correlation (spatial mean and correlation over domain are equal to ensemble mean and correlation) when the size of the spatial domain tends to be infinite in every directions.

Each realization is ergodic in mean as (see equation (44)).
For properties 3 and 4 to be true, this further condition has to be satisified:
, , as any of the , , is equal to zero.

As discussed in appendix 1, the random fiels generated from equation (2.3.2) are not ergodic in correlation as . However, using properly defined wavenumber shifts [7, 26], ergodicity in correlation is recovered on a finite domain as the spatial correlations are calculated over a domain of size . In this case, the wave number vector introduced in equation (42) is modified so as it also depends on the index , as follows:
(46)
Besides, as wavenumber shifts are introduced, the condition that functions be equal to zero as any can be removed for properties 3 and 4 to be valid.
2.3.3. Translation to nonGaussian stochastic vector fields
The approach presented above generates zeromean unitvariance homogeneous Gaussian stochastic fields , , with crosscorrelation matrix . homogeneous nonGaussian stochastic translation fields can be obtained from their Gaussian counterparts . The translation fields are defined by the following memoryless – meaning that the outputs at any point do not depend on the inputs at any other point – mapping:
(47) 
where is the standard Gaussian cumulative density function (CDF) of the random variables , the inverse of the marginal CDF of the nonGaussian random variables , and .
Then, the components of the nonGaussian correlation matrix can be computed as:
(48)  
where denotes the standard Gaussian joint PDF of the two random variables and . Note that in the case of standard Gaussian distribution, we have (see equations (29) and (30)).
In practice, one is interested in generating realizations of nonGaussian random fields with targeted marginal PDF and targeted crosscorrelation matrix . To that purpose, the crosscorrelation matrix of the underlying standard Gaussian fields has to be determined. We recall (see equations (29) and (30)) that:
(49) 
Suppose that, , we calculate from relations (2.3.3) and (49) the two quantities and . Following [11, 10], if the functions all fall in the range , , then equation (2.3.3) can be analytically or numerically inverted to calculate a unique . Besides, it must be verified that the matrix really is a correlation matrix, namely that the autocorrelation functions , , as well as the correlation matrix are positive semidefinite for every separation distance .
Inverting relation (2.3.3) is not always possible, and when it is not, crosscorrelation matrix and marginal CDFs , , are said to be “incompatible” [31]. In this case, the iterative method presented in [31] can be implemented (see also [3]). With this method, the nonGaussian CDFs are taken as and the correlation functions of the underlying standard Gaussian fields are iteratively modified until the correlation fonctions of the translated fields are sufficiently close to the targets: .
3. Numerical implementation
3.1. Random vector fields generation using FFT
For numerical implementation, equation (2.3.2) is rewritten as:
(50)  
(51)  
(52) 
for , with and for all , where is the real part of the complex number , and where we introduced:
(53) 
Relation (51) can be numerically computed in an efficient way using fast Fourier transform (FFT) algorithm.
The random fields are generated over a spatial period setting
(54) 
with to avoid aliasing. Introducing definitions (42) and (53), along with (44) and (54) in (51), and reminding that for any , that is , it comes:
(55) 
where a sequence of Fourier or inverse Fourier transforms, according to the sign of , can be recognized.
Note that in the case wavenumber shifts are applied, the equations in this section 3.1 can be straightforwardly adapted by introducing equation (46) instead of (43) for the wave numbers vector. The side effect is that the periods over which random fields are generated are elongated as:
(56) 
and with the random fields generated over the grid (compare to (53)):
(57) 
3.2. Material response at mesoscale
The components of the elasticity tensor , damageplasticity ratio and yield stress are parameters of the material model introduced in section 2.2. These parameters are realizations of random variables over the ED , according to the random vector fields generated as presented in the previous section. For the sake of readability, reference to the spatial position () and to the random experiment () are dropped in this section.
3.2.1. Discrete evolution equations
We introduce the discrete process for the pseudotime: with , , and the pseudotime increment .
The numerical integration of the evolution of the internal variables over the process is performed using the unconditionally stable backward Euler time integration scheme. Accordingly, the evolution of the internal variables (see equations (18) and (19)) are implemented as:
(58)  
(59) 
where .
Besides, considering equation (13), we have for the stress tensor:
(60)  
(61) 
Finally, the tangent modulus in equation (25) is computed as:
(62) 
3.2.2. Solution procedure
The problem to be solved at any material point of the mesoscale, reads:
Given , find such that .
This is solved using a socalled returnmapping algorithm (see e.g. [35, 13]) where a trial state is first considered and followed by a corrective step if required:

Trial state:
It is assumed that there is no inelastic evolution due to deformation increment , that is . Accordingly, the internal variables remain unchanged: and . The trial stress along with the trial criterium function can then be computed as:(63) (64) The admissibility of this trail state then has to be checked:

If , the trial state is admissible and the local variables are updated accordingly: , , . Besides, the tangent modulus is: .

If
