Oyunchimeg Shagdar
email: oyunchimeg.shagdar@inria.fr

Evaluation of Synchronous and Asynchronous Reactive Distributed Congestion Control Algorithms for the ITS G5 Vehicular Systems

Keywords: IEEE 802.11p, channel congestion, distributed congestion control (DCC), simulation IEEE 802.11p, la congestion de canal, congestion control distribué, simulation

is the technology dedicated to vehicular communications to support road safety, efficiency, and comfort applications. A large number of research activities have been carried out to study the characteristics of the IEEE 802.11p. The key weakness of the IEEE 802.11p is the channel congestion issue, where the wireless channel gets saturated when the road density increases. The European Telecommunications Standardization Institute (ETSI) is in the progress of studying the channel congestion problem and proposed so-called Reactive Distributed Congestion Control (DCC) algorithm as a solution to the congestion issue. In this report we investigate the impacts of the Reactive DCC mechanism in comparison to the conventional IEEE 802.11p with no congestion control. Our study shows that the Reactive DCC scheme creates oscillation on channel load that consequently degrades communication performance. The results reveal that the channel load oscillation is due to the fact that in the Reactive DCC, the individual CAM (Cooperative Awareness Message) controllers react to the channel congestion in a synchronized manner. To reduce the oscillation, in this report we propose a simple extension to Reactive DCC, Asynchronous Reactive DCC, in which the individual CAM controllers adopt randomized rate setting, which can significantly reduce the oscillation and improve the network performance.

Introduction

Wireless communication is expected to play an important role for road safety, efficiency, and comfort of road users [START_REF]Vehicular Communications; Basic Set of applications[END_REF]. To support such ITS applications the IEEE 802.11p [2] (ETSI ITS G5 [3] for the European usage) is standardized for V2X communications using the 5.9 GHz frequency bands. The IEEE 802.11p has been the focus of a great number of R&D activities, and its applicability to road safety and efficiency applications have been tested in some projects. The key weakness of the IEEE 802.11p is the channel congestion problem, where channel is saturated when the number of the 802.11pequipped vehicles is large. This problem is obviously due to the limited resource at the 5.9 GHz band, but also because all the vehicles are expected to periodically broadcast CAMs, which are needed for collision avoidance but tend to load the wireless channel.

A number of DCC algorithms such as Reactive DCC [4], LIMERIC [5] and AIMD adaptive control [6] have been proposed. The key differences lie in their ways of controlling the communication parameters.

Having the CAM generation rate as the control parameter, the reactive DCC controls the rate following a parameter table; LIMERIC controls the rate based on a linear adaptive algorithm, and AIMD algorithm control the rate in a similar manner as TCP.

On the other hand, channel busyness ratio (CBR), which is the ratio of the time the channel perceived as busy to the monitoring interval, is the commonly agreed metric used to characterize channel load. Since the wireless channel is shared by the ITS-S that are in the vicinity of each other, CBR monitored at such ITS-Ss take similar values. As a consequence, the ITS-Ss may take synchronized reactions to the channel load, e.g., the ITS-Ss reduce/increase the transmission rate at around the same time. The first contribution of this work is thus to study such a synchronized DCC behavior observed in reactive DCC algorithm. We pay an attention on the following different possible reactions of the CAM generator, which is responsible for adjusting the message generation rate as a means of DCC. Timer handling: In general, a transmission of a CAM is triggered by a timer, which is set to the CAM interval. Hence, upon being informed with a new CBR value (at an arbitrary point of time), the CAM generator may i) wait the expiration of the on-going timer and set the timer to the new CAM interval or ii) cancel the on-going timer and set it to the new CAM interval. We call the former and latter behaviors as Wait-and-Go and Cancel-and-Go.

Interval setting: As mentioned above, CBR measured for the shared channel may lead to the situation where the nearby ITS-Ss increase/decrease the CAM interval at around the same time. This is especially true for the reactive DCC algorithm, which controls the rate following a table. Therefore, one can think of avoiding such a synchronized behavior by applying random intervals. Hence, we can imagine 2 possible behaviors: upon reception of a new CBR value, the CAM generator sets the message generation interval to i) the value (say new_CAM_interval) provided by the table or ii) a random value (e.g., taken from the range [0, new_CAM_interval]) for the first packet and then follows the table. We call the former and latter behaviors as Synchronized and Unsynchronized. The first contribution of this work is hence, to study and compare the performances of these different versions of reactive DCC to understand the synchronization issue and the underlying reasons.

Second contribution of this work is to have a close look to channel load characterization. While it is commonly agreed that CBR should be monitored over a certain interval (e.g., 100 ms), it is not clear if channel load should be characterized with the current value of CBR or it should also consider the past CBR values. To study this aspect, we define channel load (CL) as follows.

𝐶𝐿 ! = 1 -𝛼 ×𝐶𝐿 !!! + 𝛼×𝐶𝐵𝑅 ! (1)
Here, CBRn is CBR measured at the nth monitoring interval, CLn is the channel load calculated upon measurement of CBRn. As can be seen in ( 1), the weight factor, α, plays the key role for defining whether the channel load should consider only the last CBR or should pay an attention on its history.

Obviously by choosing α=1, channel load is characterized by the "current" channel condition. In our study, we evaluate the performances of reactive DCC for different values of α.

The third contribution of this work is to study the DCC performance in road systems, which consist of ITS-Ss with different levels of sensing capability. Specifically, we consider that ITS-Ss sense the wireless channel at different levels, and as consequence, they perceive CL differently and react differently. To realize this study, ITS-Ss in the simulations are provided with random sensitivity offset values in the range of [-6, +6] dBm.

To summarize, the contributions of this work are as follows

Used Simulation Tools

The work is carried out using the open discrete event simulation environment NS3 (version 3.21) [START_REF]NS3, network simulator[END_REF],

and the traffic simulator SUMO (version 0.22) [START_REF] Sumo | Simulation of Urban mobility[END_REF]. The key simulation modules, which are relevant to this work, are illustrated in FIGURE 1, where the modules written in red are newly developed software.

NS3

The latest stable version of NS3 at the time of writing, NS-3.21, is used in this work. Among a number of new functionalities, it includes the WAVE system, which has the IEEE 802.11p (ITS-G5).

Channel load measurement

Radio Propagation module

𝐶𝐵𝑅 = ! !"#$ ! !!"#$!% (2)
NS3 mobility module is responsible for mobility of ITS-Ss, and it is the interface of NS3 with the SUMO traffic simulator.

SUMO

The SUMO traffic simulator is used to generate road network and traffic following user-specified scenarios. The outputs of the traffic simulator are converted in a file format readable by the mobility module of the NS3 simulator.

Used Simulation Tools

Unless otherwise noted, the communication and road parameters take the values listed in The parameter table of the reactive DCC algorithm is shown in Table 2.

Simulation Scenarios

The work is carried out for homogenous highway scenarios. RSUs are installed in every 100 m on the middle lane.

The scenario consists of sparse, medium, dense, and extreme density classes; the density parameters are listed in Table 4. Density parameter for homogenous highway scenario..

Performance Metrics

As described in Section 1 this work makes three contributions 1) study on synchronization issue of DCC control, 2) study on channel load characterization, and 3) study on non-identical sensing capabilities.

Following metrics are used for performance investigations. This scenario represents a typical highway, with 2 directions and 3 lanes in each direction. Even though the average vehicular density should be kept as in Table 18, there are also extra RSUs (ITS-S), which are located on the middle lane and are used to extract statistic (CL, IRT, P tx ) in constant and static locations. These RSU never transmit and, therefore, do not participate to the congestion level. The 1D highway scenario is illustrated in Figure 26 and the specific parameters are listed on Table 19. 27, where the specific configuration parameters are listed in Table 20.

• Packet delivery ratio (PDR): the ratio of the number of received packets over the number of transmitted (generated) packets. PDR is measured at individual ITS-S (vehicles and RSUs)

targeting CAMs transmitted by each mobile ITS-S (i.e., vehicles).

• Packet Inter-Reception time (PIR): time gap between consecutive CAM messages. PIR is measured at individual ITS-Ss for received CAMs from each mobile ITS-S.

• Number of transmissions: the total number of CAM transmissions is counted for 20 milliseconds of time bins.

• CBR: the average CBR is calculated for 20 milliseconds of time bins.

• Jain's fairness index is calculated for the total number of transmissions from individual mobile ITS-Ss.

Study on DCC Synchronization Issue

This section evaluates the performances of the four different versions of Reactive DCC: DccReactive-1, DccReactive-2, DccReactive-3, and DccReactive-4 for homogeneous static highway scenario. The performances of these mechanisms are compared against that of DccOff, which is the ITS-G5 system without distributed congestion control.

Evaluation of Packet Delivery Ratio

Summary

The table below lists the maximum positive and negative performance differences between the individual DccReactive and DccOff for the homogenous highway scenario. The table shows that in terms of PDR, DccReactive-4 (Cancel-and-Go & Unsynchronized) shows the best performance; but in terms of PIR the DccReactive-3 (Wait-and-Go & Unsynchronized) outperforms the other mechanisms.

Since PDR is the ratio of the number of the received and the transmitted packets, a large PDR can be obtained by aggressively reducing the number of transmissions, i.e., CAM message generations. This is not always a good solution for safety applications. Indeed FIGURE 11-FIGURE 14 clearly show that the number of transmissions of DccReactive-4 is significantly smaller than that of DccReactive-3, explaining why DccReactive-4's PDR is better than that of DccReactive-3. On the other hand, PIR is the length of time during which the receiver node does not receive data from the transmitter node. For CAM packets, this obviously implies the time gap during which the receiver ITS-S does not have information about the transmitter ITS-S. Therefore PIR is one of the key parameters that can determine whether the V2X communications protocol can support such safety applications. To this reason, by paying more attention on PIR performances, and we conclude that DccReactive-3, Wait-and-Go & Unsynchronized, mechanism is the best approach among the four versions of DccReactive.

Study on Channel Load Characterization

In this section, we study the impact of the weight factor a of the channel load defined in [START_REF]Vehicular Communications; Basic Set of applications[END_REF]. The study is made based on the performance investigations of DccReactive-1 and -3 for homogeneous static highway scenario. Note that we omit the Cancel-and-Go mechanisms (DccReactive-2 and -4), because our results presented in Section 4 show that Cancel-and-Go mechanisms show degraded performance compared to the Wait-and-Go schemes. For simplicity, we now call DccReactive-1 as SyncDccReactive and DccReactive-3 as UnsyncDccReactive.

Simulation Results

Summary

Table 6 and Table 7 list the weight factor, α (see [START_REF]Vehicular Communications; Basic Set of applications[END_REF]), which corresponds to the shortest PIR for SyncDccReactive (DccReactive_1) and UnsyncDccReactive (DccReactive_3) mechanisms. The difference of PIR between the minimum value and that for α =1 is indicated for the cases, where α is not 1. Table 6 clearly shows that for the synchronized DCC system, the best PIR performances never achieved when α =1; the difference between the minimum PIR and that for α =1 is over 10 milliseconds.

This implies that it is difficult to characterize channel load by only the current CBR. In contrast, as Table 7 shows that when the system is unsynchronized, the best PIR is achieved when α =1, indicating that channel load can be characterized by only the current CBR if the system is unsynchronized. Note that α =0.2 provides the shortest PIR for small transmitter and receiver distances in the extreme scenario for UnsyncDccReactive, the PIR difference between the minimum value that that for α =1 is very small (4ms).

Evaluation of Synchronous and Asynchronous Reactive Distributed Congestion Control Algorithms for the ITS G5 Vehicular Systems 29

Study on Non-Identical Sensing Capabilities

In this section, we study the impact of non-identical sensing capabilities. Based on the previous work, this study targets DccReactive-3 (Wait-and-Go & Unsynchronized) and DccOff mechanisms.

Simulation Results

Figure 31 -Figure 32 compare PDR of DccReactive and DccOff mechanisms for the cases where ITS-Ss have identical and non-identical sensing capabilities, respectively. The horizontal axis is the distance between the transmitters and the receivers. As can be seen in Figure 31, if the ITS-Ss have identical capabilities, the average PDR is stable when the distance between the transmitter and receiver is below 420 meters for both DccOff and DccReactive mechanisms. On the other hand, when the system consists of non-identical ITS-Ss (Figure 32), the stable distance is up to 250 meters for the sparse scenario and shorter for the medium and high density scenarios. In contrast, in the case of non-identical ITS-Ss, the fairness index degrades, and performance degradation is significant for DccReactive.

Summary

The results clearly show that when the ITS-Ss have non-identical sensing capabilities, in general, the average communication range tends to be reduced for both DccOff and DccReactive schemes. The unfair sensing capabilities result in unfair transmission behaviors. The unfairness issue is more significant for DccReactive than for DccOff.

Inria

Conclusion

In this work, we studied the following issues targeting reactive dynamic DCC algorithm.

• Synchronization • Channel load characterization • Non-identical receiver capability Following conclusions were drawn:

• It is very important to provide a solution to avoid synchronized DCC behaviour among ITS-Ss. If a careful attention is given on this issue, the simple reactive DCC algorithm can perform better than DccOff. In the case of rate adaptation, introducing a random rate seems to be a good solution. • If the road traffic is sparse, the reactive DCC algorithm tends to show poorer performance than DccOff.

• Cancelling timer for the CAM generator seems to be not necessary. • If the system is unsynchronized, it seems that only the current CBR can be a good indicator of the channel load.

However, if the system is synchronized, it is necessary to pay attention on CBR for longer interval. • If the system consists of ITS-Ss with heterogeneous channel sensing capability, non-negligible negative impact can be expected in terms of communications range and fairness. • The fairness issue caused by non-identical sensing capabilities is more significant for DCC-enabled system.

Evaluation of Synchronous and Asynchronous Reactive Distributed Congestion Control Algorithms for the ITS G5 Vehicular Systems 35

RT N° 462 Bibliography

  Considering the above-mentioned behaviors of the CAM generator, we obtain the following 4 different versions of Reactive DCC: • DccReactive-1: Wait-and-Go & Synchronized • DccReactive-2: Cancel-and-Go & Synchronized • DccReactive-3: Wait-and-Go & Unsynchronized • DccReactive-4: Cancel-and-Go & Unsynchronized

FIGURE 1 :

 1 FIGURE 1: Simulators and the key modules relevant to the work.

FIGURE 2 :Table 4 .

 24 FIGURE 2: Illustration of a homogenous highway scenario.

Figure 26 :

 26 Figure 26: Illustration of a dense highway scenario, where the measuring RSUs are uniformly distributed every 100 m in the middle lane

FIGURE 3 -

 3 FIGURE 3-FIGURE6plot the average packet delivery ratio (PDR) of the reactive DCC mechanisms in contrast to that of DccOff. The horizontal axis is the distance between the receivers and the transmitters. First of all, we note that DccOff shows an excellent PDR performance in the sparse scenario (see FIGURE3), where obviously the channel is not congested. The channel congestion becomes an issue for medium, dense and extreme density classes, where PDR degrades down to 10% in DccOff. DccReactive mechanisms show better PDR than DccOff. The PDR improvement is much more significant for unsynchronized DCC schemes (DccReactive-3 and -4) than for synchronized scheme (DccReactive-1 and -2). For timer handling, Cancel-and-Go schemes show poorer performances (DccReactive-2 in comparison to DccReactive-1 and DccReactive-4 in comparison to DccReactive-3).

FIGURE 3 :

 3 FIGURE 3: Comparison of PDR performances of DccReactive-1 and DCC-Off.

FIGURE 4 :

 4 FIGURE 4: Comparison of PDR performances of DccReactive-2 and DCC-Off.

FIGURE 5 :

 5 FIGURE 5: Comparison of PDR performances of DccReactive-3 and DCC-Off.

FIGURE 6 :

 6 FIGURE 6: Comparison of PDR performances of DccReactive-4 and DCC-Off.

FIGURE 7 -

 7 FIGURE 7 -FIGURE 10 plot the average packet inter-reception time (PIR) of the reactive DCC mechanisms in contrast to that of DccOff. Similar to the PDR case, DccOff shows an excellent PIR performance in the sparse scenario, but the performance largely degrades for higher density classes and it can exceed 1 second in the extreme density class. The reactive DCC mechanisms show better or worse PIR performances, depending especially on if the mechanism has synchronized or unsynchronized behavior. Both the synchronized schemes, DccReactive-1 and DccReactive-2, in general, show poorer performance w.r.t DccOff, except that the case of DccReactive-1 (Wait-and-Go) in the extreme density class. On the other hand, the unsynchronized schemes, DccReactive-3 and DccReactive-4, provide improved performances for dense and extreme classes. The performance improvement is significant for the DccReactive-3 (Wait-and-go & Unsynchronized) and the performance degradation is significant for DccReactive-2 (Cancel-and-Go & Synchronized).

FIGURE 7 :

 7 FIGURE 7: Comparison of PIR performances of DccReactive-1 and DCC-Off.

FIGURE 13 :

 13 FIGURE 13: The distribution of the number of transmissions for DccReactive-3 and DCC-Off schemes (during a 5-second interval for dense scenario).

FIGURE 14 :

 14 FIGURE 14: The distribution of the number of transmissions for DccReactive-4 and DCC-Off schemes (during a 5-second interval for dense scenario).

FIGURE 15 : 4 . 4

 1544 FIGURE 15: Average CBR during a 5-seconds of time interval for DccReactive-1 and DCC-OFF schemes.

FIGURE 20 :

 20 FIGURE 20: The setting (orange) and the actual (blues) values of the CAM interval and the measured CBR at a randomly selected ITS-S in the dense scenario for DccReactive-2.

FIGURE 21 :

 21 FIGURE 21: The setting (orange) and the actual (blues) values of the CAM interval and the measured CBR at a randomly selected ITS-S in the dense scenario for DccReactive-3.

FIGURE 23 -

 23 FIGURE 23-FIGURE 26 plot the PIR performances of SyncDccReactive scheme for different density scenarios. Similarly, FIGURE 27-FIGURE 30 plot the PIR performances of UnsyncDccReactive schemes for different density scenarios. Each graph has several curves for different transmitter and receiver distance ranges. Specifically, the curve for "d" indicates the results obtained for transmitters and receivers that are at the distance [d-20, d+20] m from each other. The horizontal axis is a (in percentage), and hence we are interested in a, which provides the smallest (or relatively small) PIR value. As can be seen in FIGURE23-FIGURE26, no particular value of that provides satisfying performance can be found for SyncDccReactive for all the scenarios and distance ranges. We can say however, a =1, where the algorithm considers only the last CBR value, tends to lead to poorer performances for SynDccReactive. In contrast, for UnsyncDccReactive (see FIGURE 27-FIGURE30), a =1 provides the best performance almost all the scenarios and distance ranges.

FIGURE 23 :

 23 FIGURE 23: PIR of SyncDccReactive for different values of the weight factor (α) in the sparse scenario.

FIGURE 25 :

 25 FIGURE 25: PIR of SyncDccReactive for different values of the weight factor (α) in the dense scenario.

FIGURE 26 :

 26 FIGURE 26: PIR of SyncDccReactive for different values of the weight factor (α) in the extreme scenario.

FIGURE 27 :

 27 FIGURE 27: PIR of UnsyncDccReactive for different values of the weight factor (α) in the sparse scenario.

FIGURE 29 :

 29 FIGURE 29: PIR of UnsyncDccReactive for different values of the weight factor (α) in the dense scenario.

FIGURE 30 :

 30 FIGURE 30: PIR of UnsyncDccReactive for different values of the weight factor (α) in the extreme scenario.

Figure 33 -

 33 Figure 34 compare PIR of DccReactive and DccOff mechanisms for the cases where ITS-Ss have identical and non-identical sensing capabilities, respectively. A similar observation as for the case of PDR can be made for PIR, which can reach 1.6 seconds for DccOff for non-identical sensing capability. In Figure 35 -Figure 36, Jain's fairness index is calculated targeting the number of transmissions for DccReactive and DccOff mechanisms for the cases where ITS-Ss have identical and non-identical sensing capabilities, respectively. The horizontal axis shows the density classes: 100m, 45m, 20m, and 10m represent the sparse, medium, dense, and extreme classes. In the case of the identical ITS-Ss, the fairness index is 100% for DccOff regardless of density class, and it is slightly lower for DccReactive.

FIGURE 31 :

 31 FIGURE 31: Comparison of PDR of DccReactive and DccOff schemes for the systems with ITS-Ss with identical receiver capabilities.

FIGURE 33 :

 33 FIGURE 33: Comparison of PIR of DccReactive and DccOff schemes for the systems with ITS-Ss with identical receiver capabilities.

FIGURE 34 :

 34 FIGURE 34: Comparison of PIR of DccReactive and DccOff schemes for the systems with ITS-Ss with non-identical receiver capabilities.

NS3 SUMO Demand generation Road network generation SUMO Mobility Module CAM generator DCC rate adaptation IEEE 802.11p module

  The system follows the TCP/IP communication architecture. The key software components used in our simulations are a CAM generator, UDP/IP, IEEE 802.11p, radio propagation, and mobility modules.

	The CAM generator is a newly developed module, which takes position and mobility information from
	the mobility module and periodically generates CAMs. The module is implemented with DCC rate
	adaptation algorithms. This work focuses on the reactive DCC algorithm. When the reactive DCC
	module is provided with a CL value (see (1)), it adjusts the CAM generation interval following the
	parameter table.
	Messages generated at the CAM generator processed by the UDP and IP modules, and received at the
	IEEE 802.11p MAC. While, BTP/ GeoNetworking protocols are standardized in ETSI, utilizing UDP/IP
	is equivalent to utilizing BTP/Geonetworking for the objective of studying channel congestion caused
	by 1-hop broadcast messages (CAM). It should be noted that since the header lengths of UDP/IP and
	BTP/Geonetworking are different, the necessary message length adjustment is made at the CAM
	generator such that the length of the frames transmitted on the wireless channel have the same length to
	the case using BTP/Geonetworking.
	The PHY layer of NS3 is extended with a CBR monitoring functionality, which monitors the channel
	activities and calculates CL. Since NS3 is an event-based simulator, the CBR monitoring module
	exploits the event notifications installed in NS3. In addition, the module holds a timer and calculates
	CBR in every Tmonitor interval following (2). It should be mentioned that timer setting is made
	independently at each ITS-S, and hence the CL notifications to the CAM generator is not synchronized
	among the individual ITS-Ss.

Table 1 Table 1 .

 11 Simulation Parameters.

	Parameters	Value
	Communication
	CAM default Tx rate	10 Hz
	CAM message size	400 Bytes
	Tx Power	23 dBm
	ED threshold	-95 dBm
	EDCA Queue / TC	1 DENM / 3CAM
	Modulation scheme	QPSK ½ 6 Mbps
	Antenna pattern	Omnidirectional, gain = 1dBi
	Access technology	ITS G5A
	ITS G5 Channel	CCA
	Fading model	LogDistance, exponent 2
	Road network
	Lane width	3 m
	Lanes in-flow	3
	Lanes contra-flow	3
	DCC parameters
	CBR monitor interval (Tmonitor)	100 ms
	(see (1))	1

Table 2 . DCC Reactive Parameters.

 2 Evaluation of Synchronous and Asynchronous Reactive Distributed Congestion Control Algorithms for the ITS G5 Vehicular Systems 9

	States	CL(%)	Toff
	Relaxed	0%≤CL<19%	60
	Active_1	19%≤CL<27%	100
	Active_2	27%≤CL<35%	180
	Active_3	35%≤CL<43%	260
	Active_4	43%≤CL<51%	340
	Active_5	51%≤CL<59%	420
	Restricted	CLR59%	460

  Table 3 lists the road configuration. As shown in the table and illustrated in FIGURE 2: Illustration of a homogenous highway scenario.

Table 3 . Road Configuration.

 3

	Class	Inter-vehicle distance
	Highway length	1000 m
	Lanes/Directions	3 lanes / 2 directions
	RSU inter-location	100 m
	Vehicle size	2m x 5m

Table 18 : Scenario Parameter for Scalability Test

 18

	38	ETSI TR 101 612 V1.1.1 (2014-09)

Table 19 : Specific highway configurations for scalability tests

 19

		Parameter	Value	Default
		Highway Length	1 000 m to 50 000 m	10 000 m (1 000 m if
				static)
		RSU Inter-Location 50 m to 500 m	100 m
		Vehicle size	2 m × 5 m	2 m × 5 m
		Flow density class Sparse/Medium/Dense Dense
		Contra-flow density	As Flow	Dense
		class
	7.3.1.3	2D Parking lot

In this scenario, vehicles are homogeneously distributed in a 2D space. Accordingly, this scenario uses a homogeneous 2D vehicular density as indicated in Table

18

and is adapted to fit any 2D simulation area. 7.3.2 Heterogeneous scenarios 7.3.2.1 Heterogeneous highway In this scenario, the same average density classes as indicated in Table 18 are kept as much as possible. But as vehicles move, a limited heterogeneity in the local vehicular density may be observed. It corresponds to a real highway environment and is illustrated in Figure

Table 5 . Maximum performance improvement/deterioration of DccReactive schemes w.r.t DccOff. Algorithms PDR difference (PDRDccReactive -PDRDccOff) PIR difference (PIRDccReactive - PIRDccOff) Max improvement (Positive difference) [%] Max deterioration (Negative difference) [%]

 5

	Max	Max
	improvement	deterioration
	[s]	[s]

Table 6 . The weight factor, α, which corresponds to the shortest PIR for SyncDccReactive.

 6

	Scenario		Transmitter and Receiver distance range [m]
		40 m	100	200	300	400
	Medium	0.7 (10ms	0.4 (12ms	0.4 (13ms	0.2 (16ms	0.3 (20ms shorter
		shorter than	shorter	shorter than	shorter than	than 𝑃𝐼𝑅(𝛼 = 1) )
		𝑃𝐼𝑅(𝛼 = 1) )	than	𝑃𝐼𝑅(𝛼 = 1) )	𝑃𝐼𝑅(𝛼 = 1) )
			𝑃𝐼𝑅(𝛼 = 1)
			)
	Dense	0.2 (10ms	0.1 (29ms	0.7 (47ms	0.7 (59ms	0.7 (66ms shorter
		shorter than	shorter	shorter than	shorter than	than 𝑃𝐼𝑅(𝛼 = 1) )
		𝑃𝐼𝑅(𝛼 = 1) )	than	𝑃𝐼𝑅(𝛼 = 1) )	𝑃𝐼𝑅(𝛼 = 1) )
			𝑃𝐼𝑅(𝛼 = 1)
			)
	Extreme	0.1 (61ms	0.9 (12ms	0.9 (16ms	0.9 (17ms	0.9 (17ms shorter
		shorter than	shorter	shorter than	shorter than	than 𝑃𝐼𝑅(𝛼 = 1) )
		𝑃𝐼𝑅(𝛼 = 1) )	than	𝑃𝐼𝑅(𝛼 = 1) )	𝑃𝐼𝑅(𝛼 = 1) )
			𝑃𝐼𝑅(𝛼 = 1)
			)

Table 7 . The weight factor, α, which corresponds to the shortest PIR for UnsyncDccReactive

 7

	Scenario		Transmitter and Receiver distance range [m]
		40	100	200	300	400
	Medium	1	1	1	1	1
	Dense	1	1	1	1	1
	Extreme	0.2 (4 ms	0.2 (4 ms	1	1	1
		shorter than	shorter than
		that of 𝛼 = 1)	that of 𝛼 =
			1)

Channel Load Evaluation


F Lavergne

K Sab

J Sanahuja

M Bornert

C Toulemonde

Investigation of the effect of aggregates'morphology on concrete creep properties by numerical simulations

Keywords: Modeling, Creep, Particle Size Distribution, Microstructure, Interfacial Transition Zone

. It is shown that this feature of concrete's microstructure can explain differences between creep behaviors.

Introduction

Concrete structures evolves due to time-dependent phenomena such as shrinkage and creep. Predicting basic creep is necessary to estimate the time-dependent prestress losses in thick structures such as nuclear containment buildings. For this purpose, experimental investigations were conducted by Granger [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF] in order to assess the basic creep of six concretes with various formulations.

Postprint of Cement Concrete Research doi:10.1016/j.cemconres.2015.01.003

According to this study, concretes with apparently similar formulation may exhibit different long term behaviors. For instance, the concrete from Flamanville and the one from Paluel have very similar formulation regarding cement composition, water to cement ratio and grain size distribution as shown in table [Tab. 1]. Fillers were added to the concrete from Paluel so that the total amount of fine aggregates (smaller than 100 µm) was almost equal in these two concretes. However, the major difference between the two concretes is the nature of the aggregates : crushed granite (granodiorite) for Flamanville and semi-crushed sand-lime river gravel for Paluel. Actually, the time-dependent strains of these two concretes are very different as shown in figure [Fig. 1] : The average strain measured on concrete samples from Flamanville after three years of the same compressive loading was twice as much as the one measured on samples from Paluel. Since tests revealed comparable mechanical properties for aggregates, Granger [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF] also measured the basic creep strains of the cementitious matrices without finding noticeable deviation. He concluded that the differences of overall creep properties were due to the aggregates (shape, mineralogy or interface with the cement paste).

0MPa-52MPa

Table 1: Formulation of two concretes of Granger [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF] Because experimental results are difficult to obtain, a natural way to perform parametric studies is to rely on models. The purpose of this paper is to numerically investigate the effect of size distribution and shape of inclusions on the long term creep of a matrix-inclusion material where the viscoelastic phase is the matrix. The matrix has the same behavior as a cement paste, the inclusions are perfectly bonded to the matrix and they have the elastic properties of aggregates.

One could use X-ray tomography images as inputs of 3D FFT-based compu-2 Granger [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF]. Although the formulations of concretes from Flamanville and Paluel were comparable, the basic creep strain (10 -6 MPa -1 vs days) proved different. The compressive stress was 12 MPa and the concrete was loaded at t0 = 28 days.

tations [START_REF] Escoda | Estimation of local stresses and elastic properties of a mortar sample by {FFT} computation of fields on a 3d image[END_REF] . This procedure has two drawbacks: the image must be segmented in order to differentiate the phases, and boundary effects may occur since the sample is not the unit cell of a periodic material. Filtering strategies have recently been defined to mitigate the effect of such boundary effects and get a macroscopically consistent estimate of the overall behaviour [START_REF] Bignonnet | Macroscopically consistent non-local modeling of heterogeneous media[END_REF]. This study is focused on numerically generated periodic microstructures so as to tightly control the volume fraction of inclusions, their sizes and their shapes. The generation of several microstructures with prescribed grain size distribution and spherical or polyhedral shape of inclusions is made by dense sphere packing or Random Sequential Adsorption algorithm. The 3D numerical homogenization procedure for viscoelastic microstructures proposed by Šmilauer and Bažant [START_REF] Šmilauer | Identification of viscoelastic c-s-h behavior in mature cement paste by fft-based homogenization method[END_REF] is used to compute the time-dependent response to a constant applied load. The Interface Transition Zone is described, the model of Nadeau [START_REF] Nadeau | A multiscale model for effective moduli of concrete incorporating {ITZ} water-cement ratio gradients, aggregate size distributions, and entrapped voids[END_REF] is recalled and the assumption to extend it to viscoelasticity are exposed. Then, the obtained results are compared to the predictions of mean-field homogenization.

In the first part of the paper, the microstructure generator is described, the 3D numerical homogenization method is introduced, some numerical tests are performed to validate and assess the performances of the simulations, the case of an Interfacial Transition Zone is exposed and its effect on the creep strain is numerically estimated. In the second part, the micromechanics methods are recalled and their predictions are compared to the numerical results. The paper ends with a discussion and a conclusion.

Simulation of viscoelastic matrix-inclusion microstructures

Generation of matrix-inclusion microstructures

In this section, the microstructure generator and the 3D numerical method for computing the overall viscoelastic response are presented. A matrix-inclusion microstructure is generated according to given volume fraction, sieve curve and shape of inclusions.

Scope statement

Some features are to be enforced in order to obtain a valuable matrix-inclusion microstructure generator :

• Inclusions are placed in a cubic cell and they should not overlap.

• The microstructure should be periodic.

• A tight control of the volume fraction, the size distribution and the shape of inclusions are required.

• High volume fractions of inclusions (> 50%) are needed to represent realistic microstructures of concrete.

• The simulated material should be isotropic.

• The generator should run as fast as possible.

At least two methods may fulfill these requirements. The Random Sequential Adsorption (RSA) algorithm [START_REF] Feder | Random sequential adsorption[END_REF] is the simplest method to design and implement. Its ability to handle convex polyhedra makes it useful to assess the effect of shape of inclusions. A high volume fraction of inclusions (up to 68%) may be reached for wide sieve curves. On the other hand, this method fails to simulate monodisperse grain size distribution for volume fraction greater than 30%.

The other algorithm is the one of Lubashevky and Stillinger [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] which is dedicated to pack efficiently sets of ellipsoidal inclusions. In the case of monodisperse spheres, volume fraction up to 63% can be reached. To the authors' knowledge, there is no extension of this method to polyhedral inclusions.

Determination of the number of inclusions and their volumes

Given the sieve curve (volume fraction vs particle diameter), the volume fraction of inclusions f g and the dimension l of the cubic cell, the starting point, common to both methods, is the generation of the random number of the inclusions in the cubic cell and their volumes. The sieve curve is first approximated by a piecewise linear function of the equivalent diameters on twenty classes, each class representing 5% of the total volume of inclusions. The equivalent diameter of an inclusion, d, is given by d = Let F i (d) be the probability to pick in the class i an inclusion of diameter lower than d. Since the sieve function is piecewise linear, F i (d) π 6 d 3 must be uniform. Taking into account F i (d i+1 ) = 0 and F i (d i ) = 1 , we find that F i is given by

d 20 = d min d 0 = D max d 3 d 2 V 0%-5% V 95%-100%

Volume fraction vs particle diameter

Particle diameter (mm) volume fraction

F i (d) = d -2 -d -2 i+1 d -2 i -d -2 i+1
Hence, d picked is generated by :

d picked = F -1 i (R()) = d -2 i+1 + (d -2 i -d -2 i+1 )R() -1 2
where F -1 i is the inverse function of F i and R() is the generator of the uniform random law on [0, 1[

Convex polyhedra and the RSA algorithm

A convex polyhedron is generated aiming at controlling its shape [Fig. 3,4]. A center is generated and a number of points is picked according to a Poisson law (parameter λ > 0). Then, directions are selected on the unit sphere according to the uniform random law, and distances to the center are generated according to a log-normal law (parameters µ l ∈ R, σ l > 0). The polyhedron is defined as the convex hull of those points which is computed through the C++ interface of the Qhull software [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF][START_REF]Qhull software[END_REF].

Each convex polyhedron is generated according to the previously described procedure and is placed into the cubic cell avoiding overlapping with already generated polyhedra : its position is randomly picked in the cubic cell until non-overlapping is achieved. By putting the bigger inclusions first, high volume fraction may be reached. This method has already been used to generate concrete and mortar microstructures with spherical [START_REF] Bentz | A hard core/soft shell microstructural model for studying percolation and transport in three-dimensional composite media[END_REF] or polyhedral inclusions [START_REF] Escoda | Modélisation morphologique et micromécanique 3d de matériaux cimentaires[END_REF]. Unlike [START_REF] Escoda | Modélisation morphologique et micromécanique 3d de matériaux cimentaires[END_REF], the computation of overlapping between polyhedra is performed analytically in the proposed generator. Indeed, overlapping between two polyhedra is computed thanks to the Gilbert-Johnson-Keerthi distance algorithm [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three space[END_REF] as implemented in the Bullet physics library [START_REF]Bullet physics library[END_REF]. This procedure suppresses the exclusion-zone of one voxel adopted in [START_REF] Escoda | Modélisation morphologique et micromécanique 3d de matériaux cimentaires[END_REF] and it speeds up the computations. In addition, in order to limit the potential pair of polyhedra to be tested, the cubic cell is divided into sectors : two polyhedra can potentially overlap if they share the same sector. Moreover, each polyhedra is embedded in a rectangular bounding box because testing overlapping of boxes is much easier than testing overlapping of polyhedra. A polyhedron may be present in different sectors of the cubic cell in order to ensure periodicity.

The maximum reachable volume fraction depends on the size distribution and on inclusions' shape. A wide range of inclusions' sizes allows larger volume fractions. Simulating spherical inclusions is faster than simulating unaligned elongated inclusions with sharp edges. Artificial concrete microstructures with high volume fraction (60 -68%) and polyhedral inclusions can be built with the RSA method for some size distributions [Fig. 3]. In particular, the ratio between the largest inclusions' equivalent diameter to the smallest one must be less than 20. If this ratio is too large, the number of inclusions would become prohibitive.

If this ratio is too small, reaching the given volume fraction would take too long or could not be achieved.

Lubashevsky and Stillinger algorithm for spherical inclusions

The RSA algorithm described above cannot achieve high volume fraction for monodisperse particles. For this reason, we use the Lubashevsky and Still-inger algorithm for generating matrix-inclusion microstructures with monodisperse spherical inclusions.

In our implementation, we use the algorithm described in [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] and [START_REF] Kansal | Computer generation of dense polydisperse sphere packings[END_REF] for simulating polydisperse packing of spherical inclusions. This dynamic, event-driven algorithm enables very high packing ratio of spheres. Seeds are placed in the cubic cell and their growth rate are chosen to comply with the size distribution. Each radius increases proportionally to time t and each center moves with a constant velocity. As spheres collide, the algorithm changes their velocities to avoid overlapping. Thanks to uniform speeds and growth rates, collisions, or events, can be computed analytically. There is no need for time steps since nothing happen until the next event. Therefore the procedure jumps from one event to another and the program runs tremendously fast. This method has been extended to ellipsoids [START_REF] Donev | Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. i. algorithmic details[END_REF][START_REF] Donev | Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles.: Ii. applications to ellipses and ellipsoids[END_REF], but this implementation is limited to spherical inclusions or sets of spheres. This algorithm has already been used to model composite materials in the range of elasticity [START_REF] Ghossein | A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites[END_REF].

Numerical determination of the effective viscoelastic properties of heterogeneous materials 1.2.1. The unit cell problem

The aim of this section is to numerically determine the overall viscoelastic behavior of a matrix-inclusion microstructure where the inclusions are elastic and the matrix viscoelastic.

The strain tensor ε(t) in a viscoelastic material depends on the history of stress tensor σ(t). If the constitutive law is linear, the Boltzmann superposition principle states that the material properties are defined by a compliance function (fourth order tensor), J(t, t ), such that :

ε(t) = t 0 J(t, t ) dσ dt (t )dt
If the elapsed time since loading is the only relevant parameter, the material is non-aging :

J(t, t ) = J(t -t )
The determination of the overall viscoelastic behavior of a periodic microstructure can be obtained by solving the following auxiliary problem on the periodic

unit cell V . div σ(x, t) = 0 x ∈ V ε(x, t) = t 0 J(x, t, t ) dσ dt (x, t )dt x ∈ V ε(x, t) = E(t) + ∇ s u(x, t) x ∈ V u(x, t) periodic x ∈ ∂V σ(x, t) • n(x)
skew -periodic x ∈ ∂V Here E(t) is the time-dependent overall strain, u(t, x) is the displacement field in V ,∇ s u(x, t) is its symmetric gradient, ∂V is the boundary of V and n(x) is the outer normal to ∂V . Actually, E(t) is the volume average of ε(x, t) and we denote by Σ(t) the volume average of σ(x, t).

Numerical discretizations

We use the computational method of Šmilauer and Bažant in [START_REF] Šmilauer | Identification of viscoelastic c-s-h behavior in mature cement paste by fft-based homogenization method[END_REF] developed for cementitious materials. This method which relies on the exponential algorithm [START_REF] Zienkiewicz | A numerical method of visco-elastic stress analysis[END_REF][START_REF] Taylor | Thermomechanical analysis of viscoelastic solids[END_REF][START_REF] Bažant | Creep and shrinkage law for concrete at variable humidity[END_REF] is a time-iteration procedure to solve the viscoelastic problem for the case of steady loads. It features an integration of the constitutive equations on each time step assuming a constant stress rate. It enables the time step to grow exponentially when performing a relaxation (or creep) simulation. This method has recently made use of the FFT algorithm as the solver for the unit cell tangent problem [START_REF] Šmilauer | Identification of viscoelastic c-s-h behavior in mature cement paste by fft-based homogenization method[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF].

Let ∆ε(x) = ε(x, t i+1 ) -ε(x, t i ), ∆σ(x) = σ(x, t i+1 ) -σ(x, t i ) and ∆u(x) = u(x, t i+1 ) -u(x, t i ) be respectively the strain, stress and displacement increment between time t i and t i+1 . These fields are solution of the following unit cell tangent problem : It is assumed that the stress rate dσ dt in the interval [t i : t i+1 ] is constant and that the compliance corresponds to a series of Kelvin chains [Fig. 5] :

div ∆σ(x) = 0 x ∈ V ∆ε(x) = t i+1 t i J(x, t i+1 , t ) dσ dt (x, t )dt + ε 0 (x) x ∈ V ∆ε(x) = ∆E + ∇ s ∆u(x) x ∈ V ∆u(x) periodic x ∈ ∂V ∆σ(x) • n(x) skew -periodic x ∈ ∂V where ε 0 (x) = t i 0 J(x, t i+1 , t ) -J(x, t i , t ) dσ dt (x, t )dt and ∆E = E(t i+1 ) -E(t i ). σ γ k ε C k τ k C k
J(t -t ) = t t n k=1 τ τ k e -τ -t τ k C -1 k dτ + C -1 0
where C k , k = 0, 1, ..., n are n + 1 fourth-order elasticity tensors and τ k , k = 1, 2, ..., n are n characteristic times. The integration of the constitutive equations on the time step gives the following incremental constitutive equation :

∆σ(x) = C i t (x)(∆ε(x) -ε i (x))
where C i t (x) and ε i (x) are explicitly given by :

(C i t ) -1 = C -1 0 + n k=1   1 -τ k 1 -e - t i+1 -t i τ k t i+1 -t i   C -1 k ε i = n k=1 τ k (1 -e - t i+1 -t i τ k )γ i k
where γ i k are n symmetric second-order strains corresponding to internal variables of Kelvin chains [Fig. 5]. They are defined incrementally as γ 0 k = 0 and for i ≥ 0 :

γ i+1 k = 1 -e - t i+1 -t i τ k t i+1 -t i C -1 k ∆σ + e - t i+1 -t i τ k γ i k
Notice that the choice of a series of Kelvin chains combined with the above described exponential algorithm allows for the introduction of these internal variables. Hence, storing these few internal variables enables large computations (FEM or FFT).

To solve the unit cell tangent problem at each time step, a method is required. The FFT algorithm of H. Moulinec and P. Suquet has proved a powerful tool to compute the elastic [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and viscoelastic [START_REF] Idiart | Macroscopic behavior and field fluctuations in viscoplastic composites: Second-order estimates versus full-field simulations[END_REF] response of a periodic microstructure. It takes digital images of materials and prestress as inputs. Both the basic [Alg.1] [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and the accelerated [Alg.2] [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] versions are written for this study. An isotropic reference material C * 0 is to be introduced and its Green operator Γ 0 is applied in the Fourier space. The loading can be either the prescribed macroscopic strain or the prescribed macroscopic stress as explained in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF].

Regarding implementation details, the MPI extension of the FFTW library [START_REF]Fastest fourier transform in the west[END_REF][START_REF] Frigo | The design and implementation of FFTW3[END_REF] and the PETSc library [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF] are combined. The software routinely performs its computation on several processors on shared-memory machines as well as on clusters (Ivanhoe at EDF R&D). Using less than twenty nodes and about 100-200 processes on this cluster, the size of the problem may be increased up to 700 × 700 × 700 for elasticity and 400 × 400 × 400 for viscoelasticity with 8 Kelvin chains and our computations last less than two hours (180 time steps).

Constituents' properties

Long term behavior of concretes were experimentally identified by Granger in [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF]. Characteristic times τ k and stiffness C k of Kelvin chains correspond to a creep test of a concrete (Penly) loaded at 28 days. Eight characteristic times were introduced, one for each decade, from 0.002 days to 20000 days. Properties of the cement paste are assumed to be similar. To model a viscoelastic cimentitious matrix, the elastic and viscoelastic stiffnesses of a concrete were divided by three [Tab. 2]. Its Young modulus is E matrix = 12 GPa, which is comparable to the measured ones of reconstituted cementitious matrices [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF]. A reverse identification of these values could have been performed. All Poisson's ratios are assumed to be 0. [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF] and the behavior of the defined cement paste. All Poisson ratios are assumed to be 0.2.

Studied microstructures

To explore the effect of grain size distribution and shape, the volume fraction was set to 63% and three different microstructures were defined :

• Unimodal : all inclusions are spherical and their diameters are equal. There are 150 inclusions for a volume of 8 mm 3 and their radius is about 0.2 mm.

• Bimodal : all inclusions are spherical. Two diameters coexist : about 0.2 mm and 1 mm.

• Concrete : inclusions are polyhedra and their grain size distribution is described on [Fig. 3]. Regarding their shape, crushed inclusions correspond to convex hulls of about λ = 12 points whose distances to center are picked according to a log-normal law (µ l = 0.08, σ l = 0.02) before scaling.

With these parameters, the polyhedral inclusions are similar to inclusions used in concrete mixes. Indicators describing the inclusions' shape such as elongation ratio and flakiness ratio have been defined and measured on real samples of aggregates [START_REF] Kwan | Particle shape analysis of coarse aggregate using digital image processing[END_REF][START_REF] Astm | Standard practice for description and identification of soils (visual-manual procedure)[END_REF][START_REF] Afnor | Essais pour déterminer les caractéristiques géométriques des granulats -partie 3 : détermination de la forme des granulats -coefficient d'aplatissement[END_REF][START_REF] Hu | Shape characterization of concrete aggregate[END_REF]. For an ellipsoidal inclusion with principal axes of lengths a > b > c, the elongation ratio is a/b and the flakiness ratio is c/b. The average elongation ratio of the polyhedral inclusions used in our simulations is 1.08 and their average flakiness ratio is 0.8. The elongation index and flakiness index correspond to the volume fraction of inclusions featuring respectively an elongation ratio greater than 1.8 and a flakiness ratio lower than 0.6. For normal mix design, the combined (flakiness + elongation) index for coarse aggregates must be limited to 25%, for workability reasons [START_REF] Gambhir | Concrete Technology: Theory and Practice[END_REF]. In our simulations, the elongation index is less than 1% and the flakiness index is 11%.

Examples of these microstructures are displayed on [Fig. 9].

Validation of the numerical procedure

To ensure the reliability of the method and its implementation, some basic tests were performed.

Time Discretization

Is the response of a homogeneous microstructure precisely computed ? The integration of the constitutive law on the time step may induce a numerical error. In the case of Kelvin chain, the result of a relaxation test is estimated by the Laplace-Carson transform and it serves as a reference to compute the relative numerical error for the cement paste (Fig. 6). It should be emphasized that the long term stress of a series of Kelvin chain under a constant strain is steady. Therefore the numerical error due to time discretization vanishes. Different values of the rate b = t i+2 -t i+1 t i+1 -t i of the geometric growth of the time steps were tested. In all our simulations, we chose the geometric growth b = 1.12. Therefore, 180 unit cell tangent problems have to be solved to cover the range between t-t 0 = 5×10 -6 days and t -t 0 = 21900 days where t 0 is the time of the initial loading.

Space discretization

In the case of heterogeneous microstructures, the discretization of the FFT grid also triggers numerical errors :

• The stress and strain field are stored per voxel which induces numerical errors. See [START_REF] Brisard | Fft-based methods for the mechanics of composites: A general variational framework[END_REF] for further discussions about this issue.

• The heterogeneous microstructure is projected on a grid. Therefore, there is a loss of information about the local behavior of the material. A critical comparison of numerical methods working on conforming meshes or Cartesian grids was performed in the range of thermal conduction and linear elasticity at high phase contrast [START_REF] Dunant | A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials[END_REF]. FFT solvers lead to very good approximations of overall properties since extremely fine meshes can be considered.

In order to lower the error due to the projection on a grid, for each voxel, a local volume fraction based on 64 sensing points is computed and a Reuss-like constitutive law is computed and assigned to each voxel. Therefore, the microstructure is made of black voxels (pure matrix), white voxels (pure inclusions) and gray voxels (composite). To estimate the error due to discretization, the response of the unimodal microstructure with a volume fraction of 55% is computed at various resolutions [Fig. 7 8]. The use of a low contrast between the Young modulus of the phases (< 20), the ability of parallel FFT to deal with large grids and the choice of the Reuss bound for gray voxel lead to a mesh-independent overall response.

Representative Volume Element

It is well-known that the asymptotic overall response should not depend on the generated sample neither on the size of the unit cell which must be large enough to be representative of the microstructure ( [START_REF] Sab | On the homogenization and the simulation of random materials[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] among others). For each microstructure defined in section 1.2.4, different size of the unit cell and many samples were generated. Numerical hydrostatic and shear creep tests are performed.

Two criteria were defined to ensure that the unit cell is a representative volume element :

• For a hydrostatic creep test, the unit-cell hydrostatic compliance at time t is defined as 1 3K(t) = tr(ε(t)) tr(σ) and for a shear creep test, the unit-cell shear 13 compliance is 1 2G (t) = εxy(t)

σxy The corresponding relative standard deviation of the compliance is defined as the ratio of the standard deviation to the average of the unit-cell compliance over the samples that were generated. It represents the statistical dispersion of the estimated overall time-dependent strains and it is expected to decrease as the size of the unit cell increases.

• An anisotropic indicator a is defined to check that the unit-cell compliance is isotropic. For a hydrostatic creep test, it is defined as

a = (ε dev xx ) 2 + (ε dev yy ) 2 + (ε dev zz ) 2 + 2ε 2 xz + 2ε 2 yz + 2ε 2 xy (tr(ε)/3) 2
where ε dev is the deviatoric part of the overall strain. For a shear creep test, σ xy being the direction of loading, it writes :

a = ε 2 xx + ε 2 yy + ε 2 zz + 2ε 2 xz + 2ε 2 yz ε 2
xy a is positive and a is null if and only if the overall strain corresponds to the response of an isotropic material. 14 The results of these tests are described in [Tab. 3]. The relative standard deviation of the estimated macroscopic compliances are very low and the overall response is isotropic. These features are obtained for cubic unit cells of concrete of 4 cm×4 cm×4 cm size, the largest inclusion diameter being less than 1.25 cm and the constituents' properties being those described in section 1.2.3. In the case of bimodal microstructures, a 2mm unit cell, twice as large as the largest aggregates, was sufficient to retrieve a precise estimate of the overall strain. The frame of periodic homogenization allows precise results on small unit cells in the range of elasticity [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF][START_REF] Sab | Periodization of random media and representative volume element size for linear composites[END_REF]. Gal and Kryvoruck have performed 3D elastic computations on 50mm large samples of concrete, the maximum diameter of aggregates being 15mm [START_REF] Gal | Meso-scale analysis of {FRC} using a two-step homogenization approach[END_REF] . Our simulations indicate that the size of the RVE remains small in the range of linear viscoelasticity, as long as the elastic contrast between phases is limited, while it may enlarge for phenomena like plasticity [START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF] or fatigue [START_REF] Lachihab | Aggregate composites: a contact based modeling[END_REF].

Influence of grain size distribution and shape

Responses to hydrostatic and shear creep tests are computed, the matrix being made of cement paste. The microstructures are those defined in section 1.2.4 on page 11. The volume fraction of inclusions is set at 63%.

The overall time dependent strains remain close, even for very different microstructures[Fig. 9]. The sieve curve and the shape of inclusions have little effect on the overall viscoelastic behavior of the matrix-inclusion composite. Therefore, the only relevant parameter is the volume fraction of inclusions.

Thanks to 3D computations, stress concentrations can be estimated. a shear creep test σ xy = 1MPa, the instantaneous local shear stress is above 1.37MPa in 5% of the microstructure. After 60 years, the local shear stress is above 1.71MPa in 5% of the microstructure : for prescribed constant average stress, the stress distribution in the microstructure evolves with time. The matrix tends to relax while the stress increases in inclusions [Fig. 10]. The concrete and the bimodal microstructures cannot be distinguished based on the probability distributions of the stress component σ xy during a shear creep test [Fig. 10].

The limited heterogeneity between the inclusions and the matrix may be one of the reasons for the overall time dependent strains being so close. On the first time step, the elastic loading, the inclusions are 5 times stiffer than the matrix. On the last time step, 60 years after loading, the tangent behavior of the inclusions is 14 times stiffer than the one of the matrix.

As long as concrete is modeled as a matrix-inclusion microstructure, the discrepancies between creep strains of similar concretes cannot be attributed to the sieve curve or to the shape of inclusions. The description of concrete's microstructure must be refined.

Considering the cementitious matrix as homogeneous seems arguable as there is no strong scale separation in concrete. Considering a 40 mm-wide sample of concrete and a 400 3 discretization, each voxel is 100 µm wide. At this scale, SEM images [START_REF] Scrivener | The interfacial transition zone (itz) between cement paste and aggregate in concrete[END_REF][START_REF] Diamond | The {ITZ} in concrete -a different view based on image analysis and {SEM} observations[END_REF][START_REF] Hussin | Petrography evidence of the interfacial transition zone (itz) in the normal strength concrete containing granitic and limestone aggregates[END_REF] clearly show that the behavior of the matrix may differ from one voxel to another due to fine aggregates, bundled partly unhydrated cement grains or a variable porosity [Fig. 11]. The assumption of a matrix-inclusion material might also be questionable for concrete since there sometimes exists an Interfacial Transition Zone between the inclusions and the cement paste. This feature is studied in the next section.

Introducing an Interfacial Transition Zone

The Interfacial Transition Zone is a ≈ 20 µm-thick layer of the cement paste that is more porous than the rest of the cement paste and has a different chemical composition [START_REF] Scrivener | The interfacial transition zone (itz) between cement paste and aggregate in concrete[END_REF]. It is due to the imperfect packing of cement grains (< 50 µm) on aggregates known as wall effect and it evolves during hydration. This feature may be observed on SEM images [START_REF] Scrivener | The interfacial transition zone (itz) between cement paste and aggregate in concrete[END_REF][START_REF] Diamond | The {ITZ} in concrete -a different view based on image analysis and {SEM} observations[END_REF][START_REF] Hussin | Petrography evidence of the interfacial transition zone (itz) in the normal strength concrete containing granitic and limestone aggregates[END_REF] and it might affect the mechanical properties [START_REF] Scrivener | The interfacial transition zone (itz) between cement paste and aggregate in concrete[END_REF][START_REF] Neubauer | A three-phase model of the elastic and shrinkage properties of mortar[END_REF]. Garboczi and Berryman have successfully compared their analytical model to numerical elastic computations in the range of elasticity, the ITZ being introduced as effective layers around particles [START_REF] Garboczi | Elastic moduli of a material containing composite inclusions: effective medium theory and finite element computations[END_REF]. For the elastic 3D computations of Gal and Kryvoruck [START_REF] Gal | Meso-scale analysis of {FRC} using a two-step homogenization approach[END_REF], this porous zone is merged with the aggregates. This procedure increases the volume fraction of inclusions, but it decreases their elastic stiffness, especially for small aggregates if the thickness of the ITZ is taken as uniform. For the 2D finite element computations of Grondin [START_REF] Grondin | How to consider the interfacial transition zones in the finite element modelling of concrete?[END_REF], an effective mixed interphase around each aggregate is formed by the ITZ and by a volume fraction of the bulk cement paste. Micromechanics methods may also take account of the ITZ as an effective layer surrounding inclusions [START_REF] Nadeau | A multiscale model for effective moduli of concrete incorporating {ITZ} water-cement ratio gradients, aggregate size distributions, and entrapped voids[END_REF][START_REF] Garboczi | Elastic moduli of a material containing composite inclusions: effective medium theory and finite element computations[END_REF] or as an imperfect bonding between inclusions and matrix [START_REF] Vandamme | Nanoindentation investigation of creep properties of calcium silicate hydrates[END_REF].

Representing the ITZ would not be possible a concrete scale since the size of the voxel would be 100 µm. Multiscale approaches may deliver a consistent way to set the behavior of composite voxels [START_REF] Bignonnet | Macroscopically consistent non-local modeling of heterogeneous media[END_REF], since the finer details of the microstructure cannot be accurately described at concrete scale.

To perform realistic homogenization computation, a 4 mm-wide sample of a mortar was considered, so that 10µm details such as an ITZ can be pictured . The model of Nadeau [START_REF] Nadeau | A multiscale model for effective moduli of concrete incorporating {ITZ} water-cement ratio gradients, aggregate size distributions, and entrapped voids[END_REF] described hereafter was used to estimate the local behavior as a function of the distance to the nearest inclusion. Let δ be the thickness of the ITZ.

• The local volume fraction of cement α c (d i ) depends on the distance to the nearest inclusion

d i α c (d i ) =      ᾱc 1 -0.5 δ -d i δ 2 d i < δ ᾱc d i > δ
where ᾱc is the volume fraction of cement far from the inclusions. As described below, it is computed for each microstructure so that the overall water to cement ratio w/c is equal to 0.48.

Let f (d i )δd be the volume fraction of the cementitious matrix located at a distance of the nearest inclusion between d i and d i +δd. f (d i ) is numerically estimated for each microstructure. The following equality stands :

w/c = ∞ 0 α w (x)f (x)dx ρ c ∞ 0 α c (x)f (x)dx
where ρ c ≈ 3.15 is the density of cement and α w (d i ) is the local volume fraction of water. If cement and water are the only constituents of the cement paste, α w (d i ) + α c (d i ) = 1 for all d i and :

∞ 0 α c (x)f (x)dx = 1 1 + ρ c w/c and finally : ᾱc = (1 + ρ c w/c)( δ 0 (1 -0.5 δ -x δ 2 )f (x)dx + ∞ δ f (x)dx) -1
• The elastic behavior C 0 as a function of the local water to cement ratio, equal to αw ρcαc , is provided by fit on experimental tests, as reported by Nadeau [START_REF] Nadeau | Water-cement ratio gradients in mortars and corresponding effective elastic properties[END_REF]. To set viscoelastic parameters, two assumptions were made.

2. The effect of water to cement ratio on the elastic and viscoelastic parameter are identical :

C k = A k C 0 ( αw ρcαc )
where A k is a constant. For instance, if increasing the water to cement ratio from 0.48 to 0.55 decreases the elastic Young Modulus (C 0 ) by 30%, then the drop of all viscoelastic modulus of Kelvin chains C k , k ≥ 1 is also 30%.

The assumption C k = A k C 0 ( αw ρcαc ) holds true if the cement paste is considered as a porous material which solid phase is a non aging linear viscoelastic material with a uniform Poisson's ratio. Indeed, in this case, at the scale of the porous material, the local compliance writes :

J(t, x) = J(0, x)f (t)
where f (t) is a scalar function of time t. The determination of the overall response of the cement paste as a function of porosity necessitates the resolution of a homogenization problem on the unit cell representing the porous material.

In the Laplace-Carson space (Appendix B), the constitutive equation for this unit cell problem becomes :

J * (p, x) = J(0, x)f * (p)
Let ε(x) be the elastic solution of the unit cell problem for prescribed overall strain E. Then, ε(x) is also the solution of the unit problem with prescribed overall strain E for all p in the Laplace Carson space. Therefore :

J hom, * (p) = J hom (0)f * (p)
Inverting the Laplace Carson transform gives :

J hom (t) = J hom (0)f (t)
This demonstrates that if the cement paste is a two-phase porous material, the local compliance having a uniform in time Poisson's ratio, then the cement paste has a uniform in time Poisson's ratio and the effect of time and microstructure are uncoupled.

If the effect of changing the water to cement ratio were limited to a change of porosity in a two-phase porous material, the consequence on the viscoelastic behavior would be similar to the one described above. But it is known that the microstructure of the cement paste includes porosity, elastic unhydrated cement, elastic or viscoelastic products of hydration and other additives [START_REF] Diamond | The microstructure of cement paste and concrete--a visual primer[END_REF]. The higher the water to cement ratio is, the higher the hydration degree becomes [START_REF] Cook | Mercury porosimetry of hardened cement pastes[END_REF] and the hydrated phases are assumed to be responsible for the viscoelastic behavior of cement paste [START_REF] Šmilauer | Identification of viscoelastic c-s-h behavior in mature cement paste by fft-based homogenization method[END_REF]. Therefore, the assumption

C k = A k C 0 ( α w ρ c α c ) probably
underestimates the long term effect of an ITZ. Numerical tests were performed on a 4 mm wide cubic cell depicting the microstructure of a mortar [Fig. 12]. At this scale, a voxel is approximately 10 µm wide and the ITZ can be introduced. Various thickness δ between 0 µm (no ITZ) and 50 µm were tested. The larger the ITZ is, the lower the elastic and viscoelastic stiffnesses are. The relative decreases of identified viscoelastic 21

Lamé coefficient λ (GPa) at elastic time step, t -t0 = 0, for δ = 0µm, 20µm and 50µm. )

t-t 0 (days) itz=0µm itz=10µm itz=20µm itz=50µm
Figure 12: Top : an ITZ can be pictured on the microstructure of a mortar, depending on its thickness δ. Bottom left : the local water to cement ratio is displayed as a function of the distance to the inclusion. The corresponding elastic Young modulus and the Young modulus of a 10 µm voxel overlapping with inclusions are also displayed. Bottom right : the numerical solution to uniaxial creep test depends on the thickness of the ITZ. Both the elastic response and the the long term strain are modified by the introduction of an ITZ.

stiffnesses are twice larger than the one of the elastic stiffness [Tab. 4]. Moreover, according to the previous paragraph, it is likely that this effect is underestimated by our computations.

Therefore, the presence of an ITZ is a plausible explanation of discrepancies between creep strains of various concretes.

Comparison of numerical and micromechanics methods

Semi-analytic methods

Linear viscoelastic homogenization

Micromechanics methods remain practical to upscale the overall viscoelastic properties. In the frame of non-aging creep, the Laplace-Carson transform turns the time-depend problem into several elastic problems thanks to the correspondence principle [START_REF] Lee | Stress analysis for linear viscoelastic materials[END_REF][START_REF] Mandel | Cours de mécanique des milieux continus[END_REF]. In this space, the Hashin-Shtrikman bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF], the Mori-Tanaka model [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Fisher | Viscoelastic interphases in polymer-matrix composites: theoretical models and finite-element analysis[END_REF], the three-phase model of Christensen & Lo [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF] or the n + 1-phase spherical model of Hervé & Zaoui [START_REF] Herve | n-layered inclusion-based micromechanical modelling[END_REF] deliver estimates of the macroscopic response. Coming back to the real space requires inverting the Laplace-Carson transform with numerical tools such as the Gaver-Stehfest algorithm [START_REF] Stehfest | Algorithm 368: Numerical inversion of laplace transforms [d5[END_REF][START_REF] Whitt | A unified framework for numerically inverting laplace transforms[END_REF] or the collocation method [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF][START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF] enization models may also provide the time-dependent response through time stepping procedures [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF][START_REF] Vu | A self-consistent estimate for linear viscoelastic polycrystals with internal variables inferred from the collocation method[END_REF][START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF]. The analytic method combining the Hashin-Shtrikman bound and the Laplace Carson transform is described in Appendix B.

2.1.2. The tri-sphere model : a scalar parameter to account for the morphology of the microstructure Compared to the Hashin-Shtrikman model, the tri-sphere model of de Larrard & Le Roy [START_REF] De Larrard | Structures granulaires et formulation des bétons -Concrete Mixture -Proportionning -A scientific approach[END_REF][START_REF] Roy | Déformations instantanées et différées des bétons à hautes performances[END_REF] introduces an additional parameter g * to account for the morphology of the microstructure. The cementitious matrix is divided in two parts. A minimal part fills the gaps between the aggregates at maximum packing density g * . The rest of the matrix ensures good rheological properties.

The corresponding homogenization procedure has two steps. First, the minimal part is merged with the inclusions using the upper Hashin-Shtrikman bound ( volume fraction g * ). Then, the rest of the matrix is introduced thank to a lower Hashin-Shtrikman bound [Fig. 13].

The close packing volume fraction g * arises from the compressible packing model of de Larrard [START_REF] De Larrard | Structures granulaires et formulation des bétons -Concrete Mixture -Proportionning -A scientific approach[END_REF], which extends the linear packing model of Stovall [START_REF] Stovall | Linear packing density model of grain mixtures[END_REF]. It depends on the close packing density Φ * i and volume fraction of solid material y i of each class i of the sieve curve.

• For unimodal spherical inclusions, g * = Φ * = π 3 √ 2 ≈ 0.74 • For bimodal spherical inclusions and if there exists no interaction between the two classes :

g * = min Φ * 1 -y small , Φ * 1 -(1 -Φ * )(1 -y small )
At packing limits, if y small < 0.21, large inclusions dominate and small inclusions fill in the gaps. On the opposite, if y small > 0.21, large inclusions are scarce and small inclusions are tightly packed. The bimodal distribution used in our numerical simulations has a volume fraction of solid material of small inclusions of y small = 0.5.

• A more precise knowledge of interactions between classes is required to address continuous inclusions' size distributions and irregular shapes. Typical values of g * for concretes are between 80% and 90%.

n + 1-phase spherical model to account for the Interfacial Transition Zone

The n + 1-phase spherical model of Hervé & Zaoui [START_REF] Herve | n-layered inclusion-based micromechanical modelling[END_REF][Fig. 16] extends the three-phase model of Christensen & Lo [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF]. It has been used by Nadeau [START_REF] Nadeau | A multiscale model for effective moduli of concrete incorporating {ITZ} water-cement ratio gradients, aggregate size distributions, and entrapped voids[END_REF] to upscale the elastic properties of mortars with an ITZ. The ITZ is depicted as successive layers of cement paste having different mechanical properties. This method is combined with the Laplace-Carson transform to estimate the time dependent strain of such a material.

In our case, the inputs are the same as in the numerical method : the distribution of the cement paste as a function of the distance to the nearest inclusion f (d i ) is provided by the microstructure. The link between d i , the water to cement ratio and the viscoelastic properties is described in section 1.2.8. Each layer has a width of 1µm and the corresponding volume fraction of the microstructure is estimated according to its discretized image. These volume fractions are then used as inputs for the n + 1-phase spherical model in the Laplace Carson space.

The sieve curve and the shape of inclusions change the output of this micromechanics method since they affect the distribution of the cement paste as a function of the distance to the nearest inclusion f (d i ) and the local volume fraction of cement α c (d i ).

Outputs of these models are to be compared to full field numerical results. The volume fraction of inclusions being set to 63%, a small deviation from the Hashin-Shtrikman bound in the Laplace-Carson space occurs at long term [Fig. 9]. The heterogeneity between the inclusions and the matrix may explain this feature. On the first time step, the elastic loading, the inclusions are 5 times stiffer than the matrix. On the last time step, 60 years after loading, the tangent behavior of the inclusions is 14 times stiffer than the one of the matrix.

Hashin-Shtrikman, tri-sphere model and numerical estimates

The creep strains estimated by the tri-sphere model are close to our numerical results for the bimodal and concrete microstructures. As g * is large, this model delivers estimates close to the Hashin-Shtrikman ones. But estimated time-dependent strains are different for the unimodal microstructure. Though the tri-sphere model anticipates a lower creep strain due to g * = 0.74, the numerical results for unimodal microstructures remain close to the other ones [Fig. 15]. Since concrete materials exhibit large g * , the three methods are relevant to produce estimates of the overall macroscopic strain.

2.2.3. Interfacial Transition Zone : n + 1 phase spherical model and numerical estimates Numerical uniaxial creep tests are performed on the mortar with ITZ defined in section 1.2.8 on page 21. There is a difference between the estimated overall strain estimated by the numerical method and the one produced by the n + 1phase spherical model : the later seems to slightly underestimate the effect of the ITZ [Fig. 16].

Inputs are similar up to the mechanical computations. Therefore, representing the ITZ as successive layers of coating around the inclusions is not precise for this mortar. At this scale, with this volume fraction of inclusions (50%) and the given sieve curve, a 10µm-wide ITZ is a percolated phase. Moreover, stress concentrations occur between close inclusions, that is in the ITZ. It could explain why the effect of the ITZ is larger than expected by the semi-analytical model.

A two-scale homogenization method was defined by merging the layers of the ITZ and the matrix first according to the Reuss bound. Then, the inclusions are introduced through a lower Hashin-Shtrikman bound. Even if this crude method performs slightly better than the n + 1-phase spherical model on this mortar, it still underestimates the time dependent strain. The fact that the numerical estimates are softer than the response of this crude model is a clue than the scale of the ITZ and the scale of the inclusions cannot be separated.

Discussions

There are several assumptions to question :

• Are the inclusions perfectly bonded to the matrix ?

• Is viscoelasticity the right frame to model the basic creep of concrete ?

The long-term discrepancies between these two concretes [Tab. 1] may also be partly due to the chemical reaction between calcium carbonate CaCO 3 from the limestone filler at Paluel and tricalciium aluminate C 3 A from the clinker. For the concrete of Paluel, using the model defined by de Larrard [START_REF] De Larrard | Structures granulaires et formulation des bétons -Concrete Mixture -Proportionning -A scientific approach[END_REF], the long term effect of the additional binder on compressive strength is comparable to 30kg/m 3 of cement par cubic meter of concrete.

It is well known that the time dependent behavior of concrete goes beyond viscoelasticity. The time-depend strain may be separated in different components : endogenous shrinkage, drying shrinkage, basic creep and drying creep. This divide-and-conquer strategy is the frame of EU standards Eurocode 2 and models of Granger [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF], Benboudjema [START_REF] Benboudjema | Modelisation des déformations différées du béton sous sollicitation biaxiales. application aux enceintes de confinement de batiments réacteurs des centrales nucléaires[END_REF], de Larrard [START_REF] De Larrard | Structures granulaires et formulation des bétons -Concrete Mixture -Proportionning -A scientific approach[END_REF] and B3 model of Bažant [START_REF]Creep and shrinkage prediction model for analysis and design of concrete structuresmodel b3[END_REF] among others. It eases the introduction of moisture content, structural effects and temperature effects. Rossi pledges instead for a unified strategy in which basic creep is a self-drying shrinkage under stress [START_REF] Rossi | Basic creep behavior of concretes investigation of the physical mechanisms by using acoustic emission[END_REF]. The assumption is that the creation of microcracks during a creep step (constant load level imposed during the creep test) generates water transfers which induces some additional self-drying shrinkage. This model is consistent with records of acoustic emissions during compressive creep tests [START_REF] Rossi | Investigation of the basic creep of concrete by acoustic emission[END_REF].

The basic creep of concrete is an ageing phenomenon and the current numerical model could be improved by considering ageing linear viscoelasticity. Moreover, additional eigenstrains due to shrinkage (autogeneous or drying) could be introduced, in order to better estimate the local stresses in the matrix and assess the possible occurrence of damage. The introduction of microcracking, or damage, would also help to explain non-linear creep above ≈ 40% of compressive strength. In [START_REF] Rossi | Basic creep behavior of concretes investigation of the physical mechanisms by using acoustic emission[END_REF], stress levels are above ≈ 50% of compressive strength whereas the overall biaxial prestress of containment buildings is below 30% of compressive strength in each direction (σ θθ ≈ 12M P a and σ zz ≈ 8M P a) .

Conclusion

Methods to build various artificial microstructures have been defined and a numerical frame to study the viscoelastic behavior of concrete as a matrixinclusion material is presented. Numerical results do not differ largely from the Hashin-Shtrikman estimates as long as volume fractions remain below 63% and the aspect ratio of inclusions and the contrast between phases remain low. Either the sieve curve and the inclusion shape do not influence creep or the matrixinclusion model is too simple to be accurate. Since concretes with close volume fraction of aggregates exhibit very different long term creep, considering such concretes as a two-phase matrix-inclusion materials to study creep is questionable. A heterogeneous cementitious matrix, interfacial zones or decohesion phenomena might explain such differences.

3D numerical simulations may be performed to assess the correctness of micromechanics models. For instance, it has been shown that numerical estimates of creep strains are different from the one obtained by the micromechanics n + 1 phase spherical model for modeling the Interfacial Transition Zone in a mortar. The n + 1 phase spherical model assumes that the ITZ is an isolated layer of around each inclusion whereas the numerical model accounts for the ITZ being a percolated phase in the mortar.

The gains of 3D numerical simulation are the ability to handle aging behaviors and, to a certain extent, the access to the magnitude of stress concentrations. Large scale computations are required to study concrete or other heterogeneous materials with a large representative volume element, especially if a non-linear phenomenon is introduced at the local scale. div (σ(x, p)) = 0

x ∈ V σ(x, p) = C i ε(p)

x ∈ inclusions σ(p) = (1 + pτ )C m ε(x, p)

x ∈ matrix ε(x, p) = E(p) + ∇ s û(x, p)

x ∈ V û(x, p) periodic x ∈ ∂V σ(x, p) • n(x) skew -periodic x ∈ ∂V The Hashin-Shtrikman bound is an analytical model which provides estimates of the macroscopic response of the elastic material. The Laplace-Carson transform of the macroscopic stress < σ > (p) is computed and the last stage is inverting this transform.

Figure 1 :

 1 Figure1: Results of uniaxial compressive creep tests reported by L. Granger[START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF]. Although the formulations of concretes from Flamanville and Paluel were comparable, the basic creep strain (10 -6 MPa -1 vs days) proved different. The compressive stress was 12 MPa and the concrete was loaded at t0 = 28 days.

3 6 π

 6 v, where v is the volume of the inclusion. Each class must provide a volume V 5% = 0.05 × l 3 × f g of inclusions. On each class i ∈ {0, 1, ..., 19}, the distribution of diameters is such that the approximated sieve curve s(d) is an affine function of d on the interval [d i+1 , d i [ with D max = d 0 > d 1 > ... > d 20 = d min . [Fig. 2]. The classes are filled one after another, starting from the class 0 of largest inclusions [d 1 , d 0 [ to the class 19 of the smallest inclusions [d 20 , d 19 [. While the volume of the generated inclusions in class i is less than V 5% , a new inclusion of diameter d picked is repeatedly generated in this class as explained hereafter.1 0

Figure 2 :

 2 Figure 2: The sieve curve : volume fraction vs particle diameter. This curve is approximated by a piecewise linear function on twenty classes, each class representing 5% of the total volume of inclusions.

Figure 3 :

 3 Figure 3: A section of a 4 cm-wide sample of a concrete made of coarse sand (d > 0.5 mm) and small gravel (d < 13.5 mm) built by the RSA algorithm is shown. There are 28760 polyhedra in the cubic cell and the volume fraction is 68%.

Figure 4 :

 4 Figure 4: On the left, polyhedra are rather round whereas polyhedra feature sharp edges in the middle. The RSA algorithm has been used. The sieve curve of the corresponding mortar is presented on the right.

Figure 5 :

 5 Figure 5: The rheological model is made of a series of Kelvin chains

Figure 6 :

 6 Figure 6: Time discretization. Left : Dilation strain during a hydrostatic creep test on a cement paste. The closed-form response and numerical estimate through the iterative procedure remain close. Right : Relative error on εx during an uniaxial relaxation test on a homogeneous microstructure. The numerical error depends on the base b for the geometric growth of the time step.

Figure 7 :

 7 Figure 7: Top : One microstructure of concrete of size 40 mm made of 63% of inclusions at discretization 100 3 ,200 3 and 400 3 . Bottom : The delayed responses to hydrostatic (left) and shear(right) creep tests are computed.

Figure 8 :

 8 Figure8: On the left : a microstructure made of 55% of spherical inclusions is discretized on a 500 3 grid. On the right : Numerical elastic hydrostatic tests were performed on this microstructure (E inclusions =60 GPa, Ematrix=6 GPa). The smallest voxel size corresponds to a discretization on a 700 3 grid. Assigning the Reuss-like behavior to gray (composite) voxels is less sensitive to discretization that the Voigt-like behavior.

Figure 9 :Figure 10 :

 910 Figure 9: Viscoelastic responses to hydrostatic or shear macroscopic loadings. Various matrixinclusion microstructures with the same inclusions' volume fraction lead to very close timedependent strains. The unimodal microstructures seem to creep slightly less than the other microstructures.

Figure 11 :

 11 Figure 11: A blue 100 µm grid is put on top of a SEM image from Scrivener [45]. The red lines indicate the largest inclusions. The cement paste is clearly heterogeneous over a 10 0µm voxel.

Figure 13 :

 13 Figure 13: Homogenization schemes of the Hashin-Shtrikman model (left) and tri-sphere model (right).

Figure 14 :Figure 15 :

 1415 Figure 14: Viscoelastic responses to a shear loading for different volume fraction of inclusions. Results of crushed inclusions are close to the Hashin-Shtrikman bound in the Laplace Carson space.

Figure 16 :

 16 Figure16: Results of the n + 1 phase spherical model (left)[START_REF] Herve | n-layered inclusion-based micromechanical modelling[END_REF] for a mortar is compared to numerical estimates for different thickness of ITZ δ. The n + 1 phase spherical model underestimates the time dependent strain in case of a creep test (right).

1 K 1 µ

 11 < σ(p) >= Ĉ HS (p) < ε(p) >= Ĉ HS (p)E ĈHS-(p)is isotropic and its bulk modulus and shear modulus are :K m (p) = (1 + pτ )K m ,µ m (p) = (1 + pτ )µ m K HS-(p) = K m (p) + f i i -Km(p) + 3(1-f i ) 3Km(p)+4µm(p) µ HS-(p) = µ m (p) + f i i -µm(p) + 6(Km(p)+2µm(p))(1-f i )5µm(p)(3Km(p)+4µm(p))

Table 2 :

 2 2. The behavior of inclusions is elastic isotropic (Young Modulus E = 60 GPa and Poisson's ratio 0.2 ). Mechanical behavior of a concrete identified on a creep test

		Concrete, L. Granger	cement paste
	τ k (days) Young Modulus E, GPa Young Modulus E, GPa
	Elastic	36.02	12
	0.002	1901.73	634
	0.02	1426.16	475
	0.2	546.23	182
	2	240.54	80.2
	20	112.68	37.6
	200	54.79	18.3
	2000	36.5	12.2
	20000	27.04	9.01

Table 3 :

 3 For each kind of microstructure and different size of the unit cell, many samples are built and numerical hydrostatic and shear creep tests are performed. The relative standard deviation of the estimated macroscopic compliances are very low and the overall response is isotropic.

Table 4 :

 4 Effect of the ITZ on the identified macroscopic Young modulus of C k . The effect of the interface is more important on the long term creep behavior than on the elastic behavior as shown by the relative difference to the no-ITZ case.

The viscoelastic Poisson's ratio of each Kelvin chain is the same as the elastic one.

Appendix A. FFT solvers Algorithm 1 Basic FFT scheme Initial strain field ε 0 (x) and prestress σ 0 (x) are provided Initial strain field ε-1 (ξ) is set while error eq > 10 -7 × error eq,0 do σ i (x) = C(x)ε i (x) + σ 0 (x) σi = F F T (σ) Compute error on equilibrium error eq ε(ξ

Initial strain field ε 0 (x) and prestress σ 0 (x) are provided Initial strain field ε-1 (ξ) is set while error eq > 10 -7 × error eq,0 or error bc > 10 -7 × error bc,0 or error comp > 10 -7 do

Error on boundary conditions : error bc =< e i+1 > -E Error on compatibility : In our case, elastic inclusions (volume fraction f i ) are embedded in a viscoelastic matrix modeled by a single Kelvin chain. The relaxation problem reads :

skew -periodic x ∈ ∂V In the Laplace-Carson space, for each p, this set of equation corresponds to the elastic problem :

Appendix B.2. Inverting the Laplace-Carson transform

Lots of methods are available to invert the Laplace-Carson transform. In this study, the Gaver-Stehfest formula [START_REF] Stehfest | Algorithm 368: Numerical inversion of laplace transforms [d5[END_REF] has been used :

and

Computing the binomial coefficients requires high precision and the long double type (IEEE 754, decimal on 128 bits) provided it. If M is too low, the formula lacks precision [START_REF] Whitt | A unified framework for numerically inverting laplace transforms[END_REF]. If M is too high, small errors on ĝ( k ln (2) t ) may trigger large errors on the outcome. M is set to 7.

Therefore, to estimate the response at time t, about 14 elastic computations are required. This formula does not seem to be practical for FEM since it lacks stability or precision. It is suitable as long as the numerical error in the Laplace-Carson space remains very low. This algorithm is useful when the Hashin-Shtrikman analytical formula or self-consistent estimate are computed in the Laplace-Carson space.