Relevance of a mesoscopic modeling for the coupling between creep and damage in concrete

In its service-life concrete is loaded and delayed strains appear due to creep phenomenon. Some theories suggest that micro-cracks nucleate and grow when concrete is submitted to a high sustained loading, thereby contributing to the weakening of concrete. Thus, it is important to understand the interaction between the viscoelastic deformation and damage in order to design reliable civil engineering structures. Several creep-damage theoretical models have been proposed in the literature. However, most of these models are based on empirical relations applied at the macroscopic scale. Coupling between creep and damage is mostly realized by adding some parameters to take into account the microstructure effects. In the authors’ opinion, the microstructure effects can be modeled by taking into account the effective interactions between the concrete matrix and the inclusions. In this paper, a viscoelastic model is combined with an isotropic damage model. The material volume is modeled by a Digital Concrete Model which takes into account the “real” aggregate size distribution of concrete. The results show that stresses are induced by strain incompatibilities between the matrix and aggregates at mesoscale under creep and lead to cracking.

of studies, the mechanisms involved during creep of concrete are still not well known and more particularly tensile creep which is a major issue especially in the case of special concrete structures. Three main mechanisms were recognized for tensile creep (Garas 2009;Bissonnette et al. 2007): water seepage, viscous shear and microcracking. The effects of tensile creep could be divided into beneficial and detrimental effects due, on the one hand, to the relaxation of the stresses generated by autogenous shrinkage and thermal dilation at early ages (Bissonnette and Pigeon 1995;Altoubat and Lange 2003) and, on the other hand, to micro-cracking (Cook and Haque 1974). Deformation and fracture of concrete are associated with very complicated progressive failure and it is usually accepted that the failure process under a sustained load is associated with the development and growth of microcracking (Rossi et al. 1994;Bazant and Li 1997;Barpi and Valente 2005). Thus, understanding the behavior of concrete requires the detailed assessment of creep and the part of damage. The creep strains are directly related to stresses and micro-cracks in concrete. For stress levels below 40 % of the compressive strength, creep strains are proportional to stress and the creep coefficient can be described using a formulation independent of the stress. For high stress levels, this linearity is lost and the creep coefficient is no longer independent of the stress. Finally, for stress levels above 70 % of the maximal strength, cracks propagate rapidly to failure (tertiary creep).
The modeling of the interaction between the viscoelasticity and damage has been performed by coupling linear creep models with damage models (Mazzotti and Savoia 2003;Omar 2004;Challamel et al. 2005;Reviron et al. 2007;Benboudjema and Torrenti 2008;Torrenti et al. 2008;Ozbolt 2001), or by using plasticity (Ruiz et al. 2007) and viscoelasticviscoplastic (Berthollet 2003) models in series in order to reproduce the non-linear behavior of concrete. These approaches have shown an interest because they permit to define directly the damage threshold by taking into account the influence of the viscoelastic behavior. These models suggest an extension for the writing of the equivalent strain to obtain a triaxial state by equivalence to a uniaxial state and are defined by the positive eigenvalues of the strain tensor. The creep strains are a part of the total strains and they are weighted by introducing an intensity factor β in the equivalent strain (Mazzotti and Savoia 2003): where . + represents the positive part of the strain fields, ε e the elastic strain tensor and ε v the viscoelastic strain tensor, while β is used to control the viscoelastic effects on the matrix damage and can be characterized by experimental creep tests.
With the same idea of introducing damage due to creep, some models were developed by considering the effect of creep on crack opening and softening mechanisms. Three approaches were considered: the concept of activation energy (Bazant 1995;Bazant and Li 1997), the introduction of creep models in a cohesive model to control the crack width (Carpinteri et al. 1997;Barpi and Valente 2005) and finally the combination of micromechanical models for the softening behavior of concrete and creep models (Barpi and Valente 2004). However, Le et al. (2007) have shown how the microscopic retardation and relaxation times of the viscoelastic matrix of a heterogeneous material control the creep behavior of the homogenized material.
Most of the macroscopic models require the calibration of a large number of parameters or the use of artificial parameters for weighting the viscoelasticity proportion in the damage evolution. Artificial parameters are generally used to take into account implicitly the interaction between the phases of concrete under loading. In the authors' opinion, this interaction can be modeled with a good consideration of the microstructure and a good approach to creep mechanisms involved. Mesoscopic modeling for concrete has shown a particular interest in the analysis of interactions between the cementitious matrix (cement paste, mortar) and aggregates (Grondin et al. 2007;Dupray et al. 2009;Lopez et al. 2008;Grassl and Jirasek 2010;Nguyen et al. 2011). The multi-scales approach was found very useful for the evaluation of material characteristics affected by the material characteristics of the components. Thus, the mesoscale modeling presents many advantages in the understanding of the fracture process and the corresponding effect of concrete heterogeneities.
In this work, the behavior of concrete is modeled at the mesoscopic scale by coupling a linear viscoelastic model, defined by Kelvin-Voigt chains, and a damage model, based on the microplane theory (Fichant et al. 1997). The model is implemented in the finite element code Cast3M. The volume of concrete is defined by a matrix consisting of the mortar, in which aggregates of different diameters are placed, from the largest to the smallest, according to a random process. This approach requires defining two behavior laws for the mortar and aggregates. The mechanical properties (failure and creep properties) of mortar are here experimentally characterized by direct tensile tests performed using an original device developed for this study. This device is then used to validate the model by comparing the numerical results with the creep and failure tensile tests made on concrete. Also, comparisons of the model calculations with experimental measurements made on three-point bending creep tests in a later study ) are performed. In order to quantify the part of damage occurring in concrete under creep, the acoustic emission (AE) technique is used (Denarié et al. 2006;Rossi et al. 2012;Saliba et al. 2012). It is a nondestructive testing method increasingly used in the determination of structural changes in relation with local phenomena (Mihashi and Nomura 1996;Otsuka and Date 2000;Haidar et al. 2005;Granger et al. 2007). It allows to locate events related to the formation and evolution of damage in the microstructure, to assess the state of the material during loading and thus could be used to validate numerical models.
In the first section, the non-linear viscoelastic model is presented. Then, the model parameters are determined by experimental tests on mortar. Finally, the numerical results obtained for direct tensile creep and flexural creep at the mesoscopic scale are compared to experimental measurements.

The non-linear viscoelastic model
A volume V of concrete is formed by two media: a matrix defined by the medium V m and inclusions defined by the medium V i . With the objective to model the damage due to creep of concrete, phases have a damage viscoelastic behavior law. A constant load F is applied on one of the surface boundary of V (Γ 1 ) with the unit normal vector n. This load implies local displacements fields u(y), local strain fields ε(y) and local stress fields σ (y) in each point y of V . The non-linear viscoelastic problem is written as follows: div σ (y) = 0 ∀y ∈ V (2) ε(y) = 1 2 ∇u(y) + t ∇u(y) ∀y ∈ V (4) where C ∼ (y, ε(y)) is the secant stiffness tensor depending on the local strains and damage, ε v (y) the local viscoelastic strain fields and Γ 2 the boundary of V on which macroscopic displacements conditions are applied. The chosen damage model is the isotropic model developed by Fichant et al. (1997). It allows to represent the unilateral effect and to obtain objective results independently of the length of the finite elements by controlling the fracture energy. This model is a simplified version of the microplan model of Bazant and Ozbolt (1990) and is based on the relation between the total stress and the effective stressσ (y) of the material defined by: where C ∼ 0 (y) is the initial stiffness of the material phases considered isotropic and linearly elastic, and C ∼ −1 (y, ε(y)) the compliance tensor of the damaged material. For the isotropic version of the model, the relation between the effective stress and the total stress is given by (Fichant et al. 1999): where d is the scalar value of the isotropic damage that depends only on the equivalent strain calculated according to the elastic strain tensor ε e (= ε(y) − ε v (y)): The evolution of the damage variable, due to external mechanical loads, is an exponential form: where B t represents a damage parameter to control the slope of the strain softening constitutive relation in function of the width h of the element, and ε d0 is the strain threshold. Damage increases when the equivalent strain ε eq is higher than the threshold strain ε d0 . The fracture energy G f is given by (Matallah et al. 2010): where h is the element size and f t the tensile strength.
The viscoelastic strain is defined by several Kelvin-Voigt chains. With those elements, the stress history does not have to be stored for the calculation of creep strain . For a Kelvin-Voigt unit i, the basic strain evolution is given by: whereε v (t) represents the derivative of the elementary basic creep strains, k v the stiffness and η v the viscosity of the Kelvin-Voigt unit. The creep process occurs only in the undamaged part of the material and the effective stresses are linearized for each time step. The local total viscoelastic strains are obtained by solving analytically the differential equations and are expressed as: where ε n v is the viscoelastic strain vector at time-step number n and a v , b v , c v depend upon material parameters (Reviron et al. 2007). The total viscoelastic strain is deduced from the sum of all the elementary viscoelastic strains: The local algorithm is summarized in Fig. 1.

Digital concrete model
The Digital Concrete Model has been developed with the aim to have a 'realistic' representation of cement-based materials by taking into account the random size distribution of heterogeneities (Mounajed 2002;Grondin et al. 2007). So, for previous applications on the behavior of concrete, the concrete microstructure was represented as a multi-phase material with successions of three material phases (cement paste, pores, aggregate of various sizes) spatially distributed in a random way. Each phase is characterized by a set of physical and geometrical parameters: the volume fraction, the unit size (aggregate or pore diameter), physico-mechanical properties. A specific algorithm has been developed to make a spatial and random distribution of these phases on the basis of a finite element grid.
First of all, all grid elements have mortar properties. Then, aggregates are placed in the grid from the biggest to the smallest, according to the aggregate size distribution, such that no overlapping is obtained with other particles already placed. For this, the distance between the center of gravity of the new inclusion (with a radius of R i2 ) and the inclusion already validated (with a radius of R i1 ) should be superior to the distance d min defined as the mean radius of the inclusions: The properties of the grid elements located in an aggregate area are those of aggregate. The number of inclusions is given by: where V tot is the total volume of elements and V i the volume of an inclusion. Details are given in Mounajed (2002), Grondin et al. (2007). The presented model is implemented in the finite element code Cast3M (Verpaux et al. 1988). Mechanical and basic creep parameters are identified from an experimental campaign.

Materials properties
Concrete and mortar specimens were mixed with Portland cement CPA-CEMII 42.5, crushed limestone aggregate distributed in fine sand, with a maximum size of 5 mm and a density of 2570 kg/m 3 and crushed gravel of size 5 to 12.5 mm with a density of 2620 kg/m 3 . 1.9 kg/m 3 of polycarboxylic ether based superplasticizer agent has been added for improving the workability of concrete. The same constituents were used for mortar. Table 1 shows the mix quantities of constituent materials for mortar and concrete. The mixtures were characterized by a water-to-cement ratio of 0.56. Three-point bending tests were realized on concrete beams with the dimensions 100 × 200 × 800 mm 3 (Saliba et al. 2010. Tensile tests were realized on cylinder specimens φ 110 × 450 mm 2 .

Experimental procedures
All tests are performed in a climate-controlled chamber at 50 % of relative humidity and temperature of 20 • C. For basic creep tests, the exchange of moisture is prevented by a double layer of self-adhesive aluminum paper.  The flexural creep tests are performed on frames ( Fig. 2(b)) with a capacity ranging from 5 to 50 KN (Omar et al. 2009). The load is applied by gravity with a weight and counterweight system and the displacement is measured at midspan. Three-point bending creep tests are realized on notched concrete specimens loaded at 70 and 85 % of the maximal load (F max ). F max is measured at the age of 28 days in fracture tests controlled with the crack mouth opening displacement (CMOD).
Tensile tests are conducted in a 100-kN capacity electromechanical machine ( Fig. 2(a)). Two thick steel plates are fixed on the top and the bottom of the specimens using a high modulus and strength epoxy. The adhesive used is a methacrylate adhesive obtained by mixing powder Plex 7742 (80 %) and liquid Pleximon 801 (20 %). In order to improve cohesion between the steel plates and the specimen, the latter is subjected to a slight compressive load (0.2 kN) by means of the jack. Two rings are fixed on the specimen with three screws where three LVDT gauges with an accuracy of ±1.5 µm are placed at 120°to each other. The relative displacement between the two sections is measured on a base of 280 mm at the central zone. A general view of the experimental setup is provided in Fig. 2(a). The nuclei of each sensor are suspended on three rods attached to the upper ring to overcome the bending effects that can occur if the load is not centered properly. To avoid bending due to defects in the geometry of mechanical parts and the cylinder, a ball joint in the lower attachment is designed to compensate any defects. Tensile tests were realized 48 hours later, the time for which the glue is totally dry. The load was applied with a constant rate of 0.2 kN/s, so that the maximal strength occurred in about 1 minute. During each test, load, LVDT displacement and load cell displacement are measured and recorded up to final failure with a data acquisition system. Mortar and concrete specimens are also loaded in creep at 70 and 85 % of their maximal load.
The AE method is used to follow the crack development in the microstructure during the tensile creep tests. The AE system comprises an eight-channel AE Win system, a generalpurpose interface bus (PCI-DISP4) and a PC for data storage analysis. Eight piezoelectric transducers (resonant frequency of 150 kHz) are used to convert the mechanical waves to electrical signals. They are placed on the specimen with silicon grease as the coupling agent (Fig. 3). The recorded AE amplitudes range from 0 to 100 dB. The detected signals are amplified with a 40-dB gain differential amplifier. In order to overcome the background noise, the signal detection threshold is set at a value of about 30 dB (value adjusted before every test) slightly above the measured background noise (RILEM TC212-ACD). The acquisition parameters are set as follows: peak definition time (PDT) = 100 µs, hit definition  For this analysis, the effective velocity is assumed to be a constant equal to 3800 m/s, for the analysis of AE source locations, even though there may be some variability depending on the wave propagation path. A high-pass filter with a cut-off frequency of 20 kHz, and a low-pass filter with a cut-off frequency of 400 KHz are used to eliminate mechanical and electromagnetic disturbances. Signal descriptors such as rise time, counts, energy, duration, amplitude, average frequency and counts to peak are captured and calculated by AEwin system.

Measurement of the concrete components properties
Concrete is modeled by a two-phase material considering aggregates embedded in the mortar matrix. So, the mechanical properties of mortar and aggregates have to be determined separately.
Two mortar specimens are submitted to failure at the age of two months to measure the maximum load (considered in creep tests). In this study, it is considered that the matrix governs the viscoelastic behavior of concrete. So, basic creep parameters are determined by tensile creep tests on mortar (Saliba et al., submitted). A constant load is maintained for 15 days at 70 and 85 % of the maximal strength. The displacement rate is very fast in the first few days of loading (primary creep) and then stabilizes (secondary creep). The Kelvin-Voigt parameters of the viscoelastic model are then calibrated to calculate the total strain evolution of mortar, performed for a uniform mortar volume (Fig. 4). Three Kelvin-Voigt Viscoelastic parameters  chains have been retained corresponding respectively to the characteristic times of 0.1, 1 and 10 days. Creep strain for aggregates is considered equal to zero and the mechanical characteristics are obtained from the literature. A Young's modulus of 60 GPa and a tensile strength of 6 MPa are considered (Granger 1996). The direct tensile test does not allow measuring the post-peak regime, thus the fracture energy (G f ) is taken from experimental flexural tests (Saliba 2012). All the parameters used at the mesoscopic scale are summarized in Table 2.

Application on concrete at the mesoscale
The considered mesoscale for modeling the behavior of concrete is the scale at which the material can be observed as a set of coarse aggregates embedded in a mortar matrix. Here, coarse aggregates are inclusions of a size greater than 5 mm following the experimental aggregate size distribution while mortar matrix is a mixture of finer aggregates and the cement paste (Table 3). The aggregates volume represents 37 % of the total volume of concrete. The stability of results depending on the size of the specimen and the maximal diameter of the inclusion was performed in a recent study (Grondin et al. 2007) and was respected here with dimensions equal to four times the highest inclusion diameter (L/D max > 4) and finite element size equal to 0.8 times the smallest inclusion diameter.

Influence of the random aggregates distribution
Direct tensile tests are carried out on three concrete specimens to measure the tensile strength (σ max ) of concrete, the Young's modulus and the maximum deformation at failure (ε max ) ( Table 4). Figure 5 shows that the crack in concrete can be localized arbitrarily in different positions along the length of the specimens. Does the aggregates distribution in concrete influence its failure behavior? Could the same result be obtained with the mesoscopic modeling? Grondin et al. (2007Grondin et al. ( , 2011 showed that the random distribution, if it satisfies the size conditions of the representative elementary volume, does not influence the macroscopic behavior law. But they did not study the influence on the cracking process. In order to answer these questions, calculations are then performed with the digital concrete model in direct tension for three mesh generations representing three different random aggregates distributions in the volume. The displacement is calculated between two points on the vertical axis at 14 cm around the center of the specimen similarly to the experimental measurements. Figure 6 shows a good correlation between the calculated load-displacement evolution and the experimental measurements. These results show that the mechanical parameters selected for mortar and aggregates are valid and that the mesoscopic modeling can reproduce well the concrete failure. The tensile strength varies slightly with the random dis- Fig. 7 Local damage at the mesoscopic scale in concrete at the failure time tribution of the aggregates. However, no effect was observed on the pre-peak part of the load-displacement curves where specimens were loaded in creep. The localization of the damage fields obtained with the mesoscopic modeling shows different cracking evolutions for the three meshes with different positions along the specimen (Fig. 7). The random aggregates distribution would have a non-negligible effect on the stress concentration that can move the position of the macro-crack as observed in experimental tests (Fig. 5). The digital concrete model permits also to simulate the cracking process. Damage begins to develop in a distributed manner perpendicularly to the axis of loading at the mortar-aggregates interface. Then damage increases with loading and develops in the matrix to localize at the peak in a single macro-crack.

Creep behavior of concrete in tension
The imposed constant load is equal to 70 and 85 % of the maximal strength of concrete determined by the fracture tests. Numerical simulations reproduce well the experimental creep tests (Fig. 8). The numerical creep displacement obtained with the digital concrete model is lower than creep displacement for mortar specimens which is expected as the aggregates do not creep and restrain the viscoelastic strain of the matrix. In experiments, specimens loaded at 85 % fail after few days. However, the crack is localized at the end of the specimen which may be due to stress concentration at both extremities. In addition, while concrete specimens failed, mortar specimens did not show any damage. Other calculations have been made for higher sustained loads and the tertiary creep was obtained. So, the tertiary creep could be attributed to a high interaction between the viscoelastic behavior of the matrix with the rigid behavior of aggregates. To obtain the same results as experiments for the test at 85 % of the maximal strength, more fine studies are needed. The presence of the interfacial transition zone, which is characterized by weak mechanical properties, could be simulated. Also, the matrix could be represented by the cement paste and sand grains could be placed with aggregates to have more interfaces which contribute to increase damage. Figure 9 shows the damage localization and the crack opening for concrete specimens loaded at 70 and 85 % after 15 days. Here the crack opening is calculated with the numerical method suggested by Matallah et al. (2010). Specimens loaded at 70 % show little damage localized at the mortar-aggregates interface. The same result is obtained for the specimen loaded at 85 %; however, the damage intensity is more important. This damage is due to the load application at the beginning and increases with creep. These results are confirmed by the AE measurements (Fig. 9). The number of hits increases highly during the first period which corresponds to the instantaneous deformation in correlation with the load application. Then, the number of hits increases slightly with time for the specimen loaded at 70 % of the maximal strength during secondary creep. For the specimen loaded at 85 %, an exponential increase of the number of hits is observed during tertiary creep reached due to the high initial damage at loading. AE hits are related to the release of elastic energy similarly to the strain energies released by damaged elements. So, it could be reasonable to assume, as an approximation, that the number of AE hits is proportional to the number of damaged elements (Zhu and Tang 2002;Zhu et al. 2010).

Structural response on a three-point bending test
Numerical simulations are also performed on concrete submitted to the three-point bending test and results are compared to experimental results analyzed in a previous study (Saliba et al. 2010. The loading is applied as an incremental vertical displacement of a rigid plate (linearly elastic law) fixed at the top middle of the beam (Fig. 10).
A regular mesoscopic mesh related to the middle part of the concrete specimen is generated and two homogeneous concrete blocks are attached to the left and right ends of the beam with progressively larger mesh to avoid stress concentration. The numerical simulations were conducted under plane stress condition. The parameters determined by the direct tensile tests are used for the mesoscopic mesh. For the homogeneous part, the parameters are determined by the direct tensile tests on concrete specimens. A fracture test is performed to determine the maximal flexural strength and to validate the choice of the mechanical properties defined for mortar and aggregates (Fig. 11). The results reproduced well the experimental load-CMOD curves.   The evolution of the local damage and crack opening is also visualized in the material microstructure (Fig. 12). The damage is localized at the front of the notch at the beginning of the test. Then cracks begin to join up together and cracks propagate in the mortar matrix exhibiting different realistic features as crack bridging and branching. The proposed model is also suitable for the computations of creep damage in bending. Creep compliance is assumed equal in compression and in tension which is in accordance with experimental evidences for mature concrete (Brooks and Neville 1977). To model the long-term creep, a fourth Kelvin-Voigt chain has been added with a characteristic time of 100 days. No tensile tests on mortars were performed for long time periods. So, the fourth chain was calibrated on the experimental flexural tests. Three calculations were performed with three different stiffness values for the fourth chain k v 4 = [98.1; 30; 25] GPa to show its influence on the creep displacement. A constant load is applied to the plate at 70 and 85 % of the maximal strength and deflection-time curves are plotted (Fig. 13). The kinetics and amplitude of basic creep displacement increase when k v 4 decreases. For k v 4 = 98.1 GPa the kinetic of creep follows the kinetic developed during the third phase modeled with k v 3 . So, the damage progression is not sufficient to increase the displacement as obtained in experiments. Results obtained for the two other values of k v 4 are better and results for the value 25 GPa are in good agreement with the experimental ones. So, it appears that the creep displacement depends highly on the choice of the viscoelastic parameters of the matrix. To define optimized values for all chains, more tensile creep tests have to be performed on mortar for long time periods.
The macroscopic behavior of concrete under creep loading is characterized by a progressive degradation of the material stiffness. This behavior is caused by the growth and coalescence of micro-cracks at the mortar-aggregates interfaces due to strain incompatibilities  (Fig. 14). That means that the mechanical properties of concrete are gradually and locally degraded under creep. Further, it is assumed that the concrete specimen collapses when a long continuous crack develops at the bottom of the specimen and at the same time the maximum axial strains of the elements along the main crack are higher than the corresponding ultimate strain.

Conclusions
A new model for the coupling of the viscoelastic behavior and damage is developed at the mesoscopic scale to assess creep of concrete. The parameters adopted for the simulations are identified by experimental direct tensile tests on mortar and concrete. The influence of the microstructure is modeled by taking into account the effective interactions between the concrete matrix and inclusions. This model allows understanding the physical mechanisms behind the failure of concrete under constant loading. Under tensile creep, damage induced at the moment of load application increases due to strain incompatibilities between mortars and aggregates and causes a strength decrease. Similar observations are obtained with the acoustic emission technique used during creep tests. The obtained results for the flexural creep tests show damage localization and propagation at the notch.
This work has allowed highlighting the nucleation of micro-cracks in the concrete microstructure under a sustained load which contribute to the weakness of concrete observed in a recent work . A work is in course to follow the acoustic events during the flexural creep of concrete and mortar. A fine analysis will be proposed to distinguish the type of damage according to the acoustic signature (RILEM TC212-ACD 2010). These experimental results will be presented in the next paper.