Jipeng Fei

Di Han

Xuan Zhang

Ke Li

Nicolas Lavielle

Kai Zhou

Jun Xingli Wang

Yan Tan

Jianwei Zhong

Man Pun Wan

Elyes Nefzaoui

Tarik Bourouina

Shuzhou Li

Bing Feng Ng

Lili Cai

Hong Li

Ultrahigh Passive Cooling Power in Hydrogel with Rationally Designed Optofluidic Properties

The cooling power of a radiative cooler is more than halved in the tropics, e.g., Singapore, because of its harsh weather conditions including high humidity (84% on average), strong downward atmospheric radiation (∼40% higher than elsewhere), abundant rainfall, and intense solar radiation (up to 1200 W/m 2 with ∼58% higher UV irradiation). So far, there has been no report of daytime radiative cooling that well achieves effective subambient cooling. Herein, through integrated passive cooling strategies in a hydrogel with desirable optofluidic properties, we demonstrate stable subambient (4-8 °C) cooling even under the strongest solar radiation in Singapore. The integrated passive cooler achieves an ultrahigh cooling power of ∼350 W/m 2 , 6-10 times higher than a radiative cooler in a tropical climate. An in situ study of radiative cooling with various hydration levels and ambient humidity is conducted to understand the interaction between radiation and evaporative cooling. This work provides insights for the design of an integrated cooler for various climates. KEYWORDS: optofluidic design, passive cooling, radiative cooling, hydrogel, integrated cooling structure impeded by the high relative humidity (RH) and abundant 48 thick clouds, leading to more than halved cooling power, 13 as 49 evidenced by the exceptionally high downwelling radiation 50 from the sky (up to 500 W/m 2 ). As a result, effective cooling 51 power through the atmospheric window is significantly lower 52 (down to <15 W/m 2 ) than its theoretical upper limit of ∼150 53 W/m 2 . 14 Thus far, subambient radiative cooling has hardly 54 been achieved in a tropical climate with >1000 W/m 2 solar 55 irradiance even using the best reported PRC (R solar ≈ 97% and 56 long-wavelength infrared emissivity E LWIR ≈ 96%). 15 57 As a very effective method, evaporative cooling has been 58 employed from ancient times (e.g., an Egyptian qullah) to 59 modern days (e.g., a Zeer fridge in Africa). However, it is far 60 from providing optimal cooling in terms of cooling power per 61 unit of water consumption. Very recently, a combination of 62 evaporative and radiative cooling has been attempted in a

T o achieve a sustainable development, COP26 set the ambitious target of net-zero carbon emission by the middle of the century, which urgently calls for low-carbon technologies. [1][2][3][4] Energy-related carbon emissions totaled 31.5 billion metric tons in 2020. Active cooling consumes a considerable amount of energy, especially in the summer season of most countries and in all seasons of tropical zones such as Singapore. For instance, cooling is responsible for over 60% of the electricity consumption in nonresidential buildings in Singapore, where 95% of the electricity comes from burning natural gas. 5 Thus, developing an energy-saving cooling method is crucial for achieving the COP26 goal of carbon neutrality and sustainability.

Passive cooling, free of energy consumption, can be realized by various strategies such as radiative cooling, evaporative cooling, high solar reflection, heat insulation, etc. 6 Recently, passive radiative coolers (PRCs) have attracted ever-increasing attention, which achieve cooling by emitting long-wave infrared (LWIR) radiation to outer space through the atmospheric window (8-13 μm). [7][8][9][10][11][12] Daytime subambient cooling in many temperate regions is achievable when the PRC has an ultrahigh solar reflectivity (R solar ) of more than 94% to greatly reduce solar heat gain. 7 However, effective thermal radiation to outer space in the tropical region is significantly 88 Herein, we achieved an ultrahigh cooling power of ∼350 W/ 89 m 2 with integrated cooling strategies in a rationally designed hydrogel. We carefully engineered the optofluidic properties of a hydrogel to concurrently achieve excellent solar reflection, heat isolation, LWIR emission, and water evaporation.

Effective subambient passive cooling (up to 6 °C below ambient temperature) is achieved in tropical Singapore, while a better cooling performance (up to 10 °C) is observed in other regions, including Illinois in the USA and Paris in France. Our integrated passive cooler (IPC) maintains a stable temperature adaptively regardless of the fluctuating ambient conditions, attributed to the sensitive self-tuning of evaporative behavior depending on incident energy. An experimental comparison shows that the IPC has 6-10 times higher cooling potential (power output is solar intensity dependent) in a tropical climate than that of a state-of-the-art radiative cooler, while consuming 200 times less water to achieve cooling performance comparable to that of a conventional spray cooling system. Our in situ radiative cooling study under controlled hydration levels and ambient humidity shed light on the interaction between radiative and evaporative cooling in the system. Beyond working as a single composite layer, IPC is proved to have broad engineering potential, serving as a functional composite privileged by its facile fabrication, providing insights for next-generation passive cooler design.

Design of IPC with Desired Optofluidic Properties.

We hypothesize achieving subambient passive cooling by rational integration of various strategies in a single composite material through one-step formation to tackle the issues above. Water evaporation, leading to evaporative cooling, can be 141 triggered by several factors including solar irradiance, ambient 142 temperature, relative humidity, wind, etc. Atmospheric radiation, solar illumination, and parasitic heat transfer can be compensated well by evaporative cooling power generated by water evaporation. Apart from the incident energy, water evaporation continuously takes away internal thermal energy of a passive cooler under the vapor pressure difference between an IPC and the ambient environment. Such a water evaporation driven by unsaturated vapor pressure is inevitable as long as the relative humidity is not 100% regardless of other climate conditions, leading to constant heat reduction. Most important, this environmental-condition-dependent behavior of water evaporation leads to an adaptive cooling performance and thus a steady temperature profile.

Optical and Structural Characterizations of IPC.

An IPC (Experimental Section in the Supporting Information) was synthesized by dispersing BaSO 4 particles (BSPs) into a f1 typical PVA hydrogel matrix. As illustrated in Figure 1a, PVA polymeric chains (gray) tightly stabilize water (light blue) through hydrogen bonding, where BSPs (white) within the polymeric network backscatter UV-vis-NIR light (300-1300 nm, ∼87.5% of solar energy). [START_REF] Ustin | Estimating canopy water content from spectroscopy[END_REF][START_REF] Fei | Mechanistic Investigation of Electrostatic Field-Enhanced Water Evaporation[END_REF] OH bending (water and PVA) and SO 4 vibration endow the IPC with high IR emissivity to facilitate radiative cooling. Compared to bare surfaces (Figure 1b), the IPC-coated surface can effectively minimize incoming solar heat gain. The active cooling demand (ACD) is the energy required to maintain a constant indoor temperature by an active cooling system (e.g., air conditioner), S12). Cooling performances are compared between IPC and conventional spray cooling in Figure S13 (Supplementary Note 10). Though the surficial temperature of bare aluminum (Al) decreased by 25 °C under a continuous water spray (13.8 mL/(m 2 s)), its temperature is still ∼2 °C above that of IPC-2 under high solar irradiance (800-1000 W/m 2 ). Moreover, the fluctuating Al temperature is in obvious contrast to the stable temperature of the IPC. Importantly, the water consumption is >200 times higher than that of IPC-2 (∼0.05 mL/(m 2 s)). It is worth noting that the surface temperature rose immediately after most of the water evaporated, suggesting nonsustainable cooling by spraying water.

IPC Cooling Capability Evaluation. The lower indoor (to exclude radiative cooling) temperature of IPC-2 than that of Pure-Gel (Figure 2e and Figure S14 and Supplementary Note 11) proves the effectiveness of enhanced fluidic properties on evaporative cooling. Apparently, 3000 sun solar illumination (1000 W/m 2 ) heated Pure-Gel to 3.68 °C above ambient temperature. In contrast, IPC-2 remained at a subambient temperature. Notably, the thickness of IPC is negligible in affecting surface temperatures (Figure S15) but is responsible for increasing the length of each cooling cycle (Figures S16 andS17).

An in situ study was conducted to investigate the interaction between radiation and evaporation, which helps to understand the rationality of the multifunction surface design in the IPC. It is clear that continuous water evaporation negligibly affects surficial emissivity within the atmospheric window, indicating the coexistence of excellent radiative and evaporative features (Figure 2f). Correspondingly, a near-steady-state thermal analysis derived net power density (sunlight and environment)

shows that the IPC-2 coating significantly reduces the ACD of rooftop by a >81.5% to maintain a comfortable indoor temperature (25 °C). Compared to the ideal PRC, IPC-2 leads to a lower ACD when it is thicker than 6 mm and possesses the potential to fully minimize the ACD when it gets thicker over a long run (Figure S18a). [START_REF] Li | Full Daytime Sub-ambient Radiative Cooling in Commercial-like Paints with High Figure of Merit[END_REF][START_REF] Wang | A structural polymer for highly efficient all-day passive radiative cooling[END_REF] Heat isolation is enhanced when IPC thickness increases, leading to a more effective rejection of incoming atmospheric thermal energy. Though the ACD is reduced effectively, the cooling efficiency (the ratio of cooling power density to thickness) starts to decrease because thickness plays a negligible role in decreasing solar heat gain. S18a) with 296 acceptable ACD without apparently compromising cooling 297 efficiency (Figure S18b).

298

To further evaluate the cooling capability of the IPC, the 299 theoretical energy output from the thermodynamic model is 300 provided (Figure 2g). Obviously, the evaporative energy 301 consumption varies significantly with solar irradiance, while 302 the radiative energy output slightly changes and remains at a 303 low level due to the unfavorable atmospheric hindrance (high 304 RH). Notably, ∼2-5% solar absorbance totally compensates 305 for the radiative cooling power of the state-of-the-art radiative 306 cooler in the tropics when the solar intensity reaches 1000 W/ 307 m 2 . As such, no subambient passive cooling could be achieved 308 by LWIR radiation alone. Interestingly, the higher and stable 309 cooling power exhibited by the IPC under stronger sunlight is 310 obviously different from that of a conventional radiative cooler, 311 which indicates its high performance and adaptive cooling 312 feature. As such, the IPC automatically tunes the contributions 313 of evaporative cooling to maintain a stable temperature, 314 leading to adaptive passive cooling under various conditions 315 regardless of the varying weather or indoor/outdoor conditions. In comparison with recent reported works, we experimentally achieved a record-high cooling power of approximately 350 W/m 2 on a typical sunny day, which is superior to those of state-of-the -art passive coolers with various designs (Figure 2h,i). [START_REF] Feng | Bilayer porous polymer for efficient passive building cooling[END_REF][START_REF] Li | A tandem radiative/evaporative cooler for weather-insensitive and high-performance daytime passive cooling[END_REF][START_REF] Yao | Integrated radiative and evaporative cooling beyond daytime passive cooling power limit[END_REF][START_REF] Xu | Nanocomposite hydrogel for daytime passive cooling enabled by combined effects of radiative and evaporative cooling[END_REF][START_REF] Yang | Phase Change Material Enhanced Radiative Cooler for Temperature-Adaptive Thermal Regulation[END_REF] In comparison with the reported multilayer structure, direct optical modification of the hydrogel matrix ensured mechanical stability since a interfacial adhesion problem between the radiative and evaporative layers may occur when the hydrogel matrix suffers great water loss.

The higher cooling power of IPC-2 is attributed to a larger evaporative surface without covering a layer of a radiative cooler.

Cooling Performance Test. 24 h indoor cooling (∼3

°C below ambient temperature) was achieved with 4 mm thick IPC-2 (Figure S19), in apparent contrast to a PRC that does not work indoors due to a blocked atmospheric window and indicated the effectiveness of evaporative cooling driven by unsaturated vapor pressure. Strikingly, outdoor cooling with ∼6 °C below ambient temperature in the early afternoon of a sunny day (solar irradiance up to 1140 W/m 2 at RH > 60%) was achieved, in apparent contrast to the near-ambient S20a and Supplementary 339 Note 12). Moreover, a greater contrastive cooling performance 340 was demonstrated between a commercial radiative cooling 341 paint (Nippon Solareflect Si) and the IPC, where the optical 342 performances are quite similar between each other (com-343 parable solar reflectance of ∼88%, Figure S20b-d). 15 This 344 directly reveals the cooling contribution from evaporative 345 behavior. Moreover, the temperature of the cooled substrate 346 remained very stable at 28 ± 0.5 °C, effectively avoiding the 347 drastic temperature variation caused by fluctuating solar heat 348 gain (Figure 3a). It is worth noting that a state-of-the-art PRC 349 shows 2 °C higher temperature on sidewalls than that on a 350 rooftop, which is at near-ambient temperature (Figure 3b) due 351 to the halved radiation angle and increased incident thermal 352 energy from the surroundings. [START_REF] Xu | Broadband directional control of thermal emission[END_REF] In contrast, IPC-2 shows 4-6 353 °C subambient temperature on sidewalls and a rooftop, which 354 are attributed to the nondirectional feature of evaporative 355 behavior. It is proved that a high NIR reflective material for 356 sidewall application leading to inevitable thermal discomfort 357 would be alleviated by the IPC concept that partial solar 358 energy would be directly transferred into inner energy of water 359 vapor through a phase change process due to solar absorption 360 of interfacial water. [START_REF] Yang | Environmental impacts of reflective materials: Is high albedo a 'silver bullet' for mitigating urban heat island? Renew[END_REF] Though water evaporation is strongly 361 affected by relative humidity, nighttime subambient (3-5 °C) 362 cooling was still achieved despite a high RH of ∼80% (Figure 363 S21). To further illustrate the adaptive cooling feature and 364 durability of the IPC cooler, we conducted a 24 h continuous 365 field test, during which various atmospheric conditions 366 appeared including strong/weak sunshine, raining, cloud, and 367 nighttime weather conditions. As shown in Figure 3c, 6 °C 368 cooling was achieved with intense sunlight in the daytime. 369 Even when it was raining, the passive cooler still reached a 370 subambient temperature, and the temperature fluctuation was 371 much smaller than that of the air temperature due to its 372 adaptive cooling mechanism. This indicates that the IPC 373 exhibits superior sustainable cooling performance. PRCs 374 typically appear white to reflect the entire visible solar 375 spectrum and colored PRCs that have broader applications are preferred. Coloring inevitably results in increased solar heat gain that raises surface temperatures to above ambient. [START_REF] Chen | Colored and paintable bilayer coatings with high solar-infrared reflectance for efficient cooling[END_REF][START_REF] Sheng | Colored Radiative Cooler under Optical Tamm Resonance[END_REF][38] In contrast, colored IPC-2 still shows a subambient (∼3.5 °C, Figure 3d and Figure S22) cooling performance under ∼1000 W/m 2 thanks to the rationally integrated cooling strategies.

Since the tropical climate is the most challenging for passive cooling, and our IPC cooler works well in tropical Singapore, we anticipate it would work equally well or better in other climates. To verify this, the cooling performance of the IPC was measured on two different continents (America and Europe). Recorded temperature profiles in Illinois, USA (Figure 3e) and Paris, France (Figure 3f) indicate a subambient cooling performance up to 10 °C with the IPC.

These results suggest a generic cooling performance and excellent adaptivity.

Next, we further evaluated the applicable potential of the IPC design. IPC holds promise as cooling pads due to its nonhazardous nature (Figure S23a) and excellent mechanical flexibility and strength. 2 mm thick IPC-2 can be conformally f4 attached to the skin (Figure 4a and Figure S22b-d The soaked fabric filled with an IPC precursor spontaneously forms a porous network within the framework, holding its original appearance but gaining a cooling capability as an IPC.

As shown, the composited fabric-IPC exhibits great flexibility and antipuncture properties (Figure 4b-iii). Thanks to the abundant hydrogen bonding interactions between IPC and fabric bones, strong adhesion is achieved with supreme mechanical performance, which is further confirmed by the lack of falling off under strong bending and pressing (Figure 4b-iv). The demonstrated composition integrated with a commercial fabric provides better human comfort. Beyond this, we demonstrated an effective cooling performance in a commercial food box, where the IPC serves as an interlayer (replacing a thermal isolative layer) to provide a cooling function (Figure 4c). An obvious temperature reduction of up

Figure 1 .

 1 Figure 1. IPC design and characterizations. (a) Schematic illustration of cooling mechanisms. (b) Illustration of active cooling demand (ACD) difference with bare (left) and IPC-covered (right) concretes. (c) Simulated scattering efficiency (300-1300 nm) with varying BSP sizes. (d) Structural illustration (upper) and scanning electron microscopy image (lower left, scale bar is 1 μm) of IPC and BSP size distribution (lower right). (e) Solar reflectance and (f) transmittance spectra of pure-Gel and IPC-2. (g) Optical (upper) and IR image (lower) of IPC-2 under sunlight taken after 30 min stabilization. (h) LWIR emittance spectrum of IPC-2. The shaded area labels the atmospheric window wavelength range.

Figure 2 .

 2 Figure 2. Passive cooling mechanisms and performance evaluation of IPC. (a) Comparison of water content and bound water ratio within Pure-Gel and IPC-2. (b) Storage and loss moduli of Pure-Gel and IPC-2. (c) DSC thermogram of water evaporation within Pure-Gel and IPC-2. Insets: schematic structures. (d) Raman OH stretching mode of Pure-Gel (upper) and IPC-2 (lower). (e) Indoor evaporative cooling comparison and temperature increase comparison (under solar illumination) between Pure-Gel and IPC-2. (f) In situ emissivity change along continuous water evaporation in a constant environment revealing the possible coexistence of radiative and evaporative cooling within an individual panel. The inset is a schematic illustration of the test bed. (g) Solar intensity dependent output power density by evaporative (P evp ) and radiative (P rad ) cooling and the power difference (P net ) between incident and output energies of IPC-2, which indicates the near static state and accuracy of the thermal dynamic model. (h) Experimental cooling power test of IPC-2 in a typical sunny day in Singapore (tropical climate). (i) Cooling power comparison among IPC-2 and state-of-art multifeature passive coolers under different RHs.

Figure 3 .

 3 Figure 3. Cooling performances of IPC-2 under various weather conditions. (a) Outdoor cooling performance of IPC-2 in a typical sunny daytime in Singapore. (b) Outdoor cooling performance of IPC-2 compared with the state-of-the-art PRC at different faces. (c) Day-long temperature recording of ambient air and passive cooler, with no water replenishment during the test. The tested IPC-2 sample has dimensions of 25 cm (L) × 20 cm (W) × 4 mm (T). (d) Outdoor subambient cooling achieved by colored IPC-2 samples (lower panel inset), indicating effectiveness of the integrated cooling structure. (e) Outdoor cooling performance test of IPC-2 at University of Illinois Urbana�Champaign in Illinois, United States on May 25th, 2023. (f) Outdoor cooling performance test of IPC-2 at Sorbonne University Pierre and Marie Curie Campus in Paris, France, on May 25th, 2023.

  , resulting in an ∼28 °C skin temperature despite the high ambient temperature (∼33 °C) under 1 sun irradiation. Moreover, IPC could be easily composited into a commercial fabric through a facile one-step cross-linking process (Figure 4b-i,ii).

Figure 4 .

 4 Figure 4. Demonstration of promising applications and engineering potential of IPC. (a) IR image showing outdoor subambient cooling of IPC on human skin (the same sample as that in Figure S22c). (b) Demonstration of direct composition between IPC and commercial fabric (10 cmL × 7 cmW × 4 mmT). (c) Schematic illustration of IPC as a composite cooling layer in a commercial food box. (d) Temperature recording of IPC cooling for a commercial food box. (e) IR image of concrete slabs inside two food boxes at 2:50 pm, August 4th, 2022, in Singapore (Figure 4d).

Nano Letters pubs.acs.org/NanoLett Letter

	Supplementary Note 9), suggesting less intermediate water
	content that leads to a higher evaporative enthalpy in the IPC
	(Figure 2c). The distinct water states fundamentally differ-
	entiate the evaporation behavior of water in IPC from that of
	normal water, which leads to fine control of evaporation
	behavior dynamically, resulting in a drastically increased
	evaporative cooling efficiency. Specifically, with the same
	amount of evaporated water, IPC-2 saves 92% ACD on a sunny day (solar intensity, 1000 W/m 2 ; RH, 60%; wind speed,
	∼1.5 m/s). Moreover, IPC-2 shows a 50% longer cooling cycle
	than normal water, indicating a much more sustainable cooling
	effect (Figure

https://doi.org/10.1021/acs.nanolett.3c03694 Nano Lett. XXXX, XXX, XXX-XXX C higher evaporative enthalpy is more favorable to cooling since it can also prolong the cooling time. It is worth noting that both mass and thermal transfer are considered in the theoretical model, which is further validated by experimental results, proving its accuracy. The Raman OH stretching mode shows a lower intermediate/free water (IW/FW) ratio in IPC-2 than in Pure-Gel (Figure 2d and Figure S11 and

Nano Letters pubs.acs.org/NanoLett Letter

The intersection of cooling efficiency and reciprocal of ACD indicates 10 mm to be an appropriate thickness for outdoor https://doi.org/10.1021/acs.nanolett.3c03694

Nano Lett. XXXX, XXX, XXX-XXX D 295 cooling (under ambient conditions for Figure

Nano Letters pubs.acs.org/NanoLett Letter

  °C was recorded under strong solar illumination outdoors, 416 where the IR images of concrete slabs (taken when the food 417 box is open) further prove the cooled atmosphere even in the 418 absence of inner convection (Figure4d,e). Additionally, the 419 IPC shows excellent fire retardant performance (FigureS24420 and Supplementary Note 13) thanks to the contained water 421 and added BSPs.[START_REF] Xue | Pyrolysis products and thermal degradation mechanism of intrinsically flameretardant carrageenan fiber[END_REF] Videos 1 and 2 clearly show that the inner 422 structure of IPC was maintained after long-term burning from 423 both vertical and parallel directions; only a thin char layer 424 formed at the surface. Altogether, IPC opens a route for 425 adaptive passive cooling for various potential applications in a 426 tropical climate and beyond (Supplementary Note 14).427We have successfully engineered a hydrogel-based IPC 428 composite that provides a solution to the subambient passive 429 cooling challenge in harsh tropical climates. Incorporating a 430 variety of mechanisms such as high solar reflection, heat 431 insulation, long-wave infrared (LWIR) radiation, and evapo-432 rative cooling, the composite delivers a remarkable cooling 433 power of 350 W/m 2 and achieves 4-8 °C subambient passive 434 cooling in tropical climates. Demonstrating excellent adapt-435 ability to fluctuating ambient conditions, our IPC cooler 436 maintains stable subambient temperatures despite changing 437 weather conditions. It exhibits superior cooling performances 438 in climates beyond the tropics, recording up to 10 °C 439 subambient temperatures in America and Europe. Additionally, 440 the IPC can work as a colored passive cooler which could 441 achieve subambient cooling under strong sunlight, which is 442 challenging for a pure radiative cooler. Given its cost-443 effectiveness, adaptability, abundant choice of spectral 444 regulation agents, and high engineering flexibility, the IPC 445 cooler holds great promise for wide application across various 446 climates.

https://doi.org/10.1021/acs.nanolett.3c03694 Nano Lett. XXXX, XXX, XXX-XXX F 415 to 7 447 ■

ASSOCIATED CONTENT 448 * sı Supporting Information 449

  The Supporting Information is available free of charge at 450 https://pubs.acs.org/doi/10.1021/acs.nanolett.3c03694. School of Mechanical and Aerospace Engineering 458 and School of Electric and Electronic Engineering, Nanyang

	451	Char test of the IPC under a vertical condition (MP4)
	452	Char test of the IPC under a parallel condition (MP4)
	453	Details of sample preparation, theoretical simulations,
	454	and supporting experimental results (PDF)
	455 ■ AUTHOR INFORMATION 456 Corresponding Author
	459	Technological University, Singapore 639798, Singapore;
	460	CINTRA CNRS/NTU/THALES, UMI 3288, Singapore
	461	637553, Singapore; orcid.org/0000-0002-6975-7787;
	462	Email: ehongli@ntu.edu.sg
	463 Authors
	464	Jipeng Fei -School of Mechanical and Aerospace Engineering,

457 Hong Li -465 Nanyang Technological University, Singapore 639798, 466 Singapore 467 Di Han -School of Mechanical and Aerospace Engineering, 468 Nanyang Technological University, Singapore 639798, 469 Singapore 470 Xuan Zhang -Department of Energy and Power Engineering, 471 School of Mechanical Engineering, Beijing Institute of

472 Technology, Beijing 100081, People's Republic of China Ke Li -Institute of Materials Research and Engineering, Agency for Science, Technology and Research, Singapore 138634, Singapore; School of Materials Science and Engineering, Nanyang Technological University, Singapore 639798, Singapore; orcid.org/0000-0002-3140-3983 Nicolas Lavielle -Univ Gustave Eiffel, CNRS, ESYCOM, Marne la Vallée F77454, France; orcid.org/0000-0001-5842-2992 Kai Zhou -Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, United States; orcid.org/0000-0003-4653-6436 Xingli Wang -CINTRA CNRS/NTU/THALES, UMI 3288, Singapore 637553, Singapore Jun Yan Tan -School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore 639798, Singapore Jianwei Zhong -Pillar of Engineering Product Development, Singapore University of Technology and Design, Singapore 487372, Singapore; orcid.org/0000-0001-6198-868X Man Pun Wan -School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore 639798, Singapore Elyes Nefzaoui -Univ Gustave Eiffel, CNRS, ESYCOM, Marne la Vallée F77454, France Tarik Bourouina -Univ Gustave Eiffel, CNRS, ESYCOM, Marne la Vallée F77454, France; CINTRA CNRS/NTU/ THALES, UMI 3288, Singapore 637553, Singapore; orcid.org/0000-0003-2342-7149 Shuzhou Li -School of Materials Science and Engineering, Nanyang Technological University, Singapore 639798, Singapore; orcid.org/0000-0002-2159-2602 Bing Feng Ng -School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore 639798, Singapore; orcid.org/0000-0001-8112-1151 Lili Cai -Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, United States; orcid.org/0000-0003-1222-248X Complete contact information is available at: https://pubs.acs.org/10.1021/acs.nanolett.3c03694


Hicham Benniou

Abdallah Accary

Yann Malecot

Matthieu Briffaut

Laurent Daudeville
email: laurent.daudeville@univ-grenoble-alpes.fr

Discrete element modelling of concrete under high-stress level: influence of saturation ratio

Keywords: DEM, Discrete element model, Concrete, Saturation ratio, Confined compression

The discrete element model proposed in this paper addresses the macroscopic behavior of concrete taking into account the presence of free water in pores, thanks to a new interaction law between spherical discrete elements (DE). When concrete structures are subjected to a severe loading, e.g. an impact, material exhibits high triaxial compressive stresses which are highly influenced by the saturation ratio. In this new constitutive model, cracking and compaction are modeled at the interaction level between DEs and free water effects are taken into account by introducing a dependency between the water saturation ratio and the inelastic deformation due to the pore closure. The present numerical model has been implemented within the YADE (Yet Another Dynamic Engine) code in order to deal with extreme loading situations leading to stress states characterized by a high mean stress level.

Introduction

When a concrete structure is subjected to an intense loading, e.g. an impact, the material in the vicinity of the loading zone undergoes high levels of stress leading to irreversible compaction, whereas farther from this location, compression with a moderate triaxial stress level occurs [START_REF] Gran | In-target radial stress measurements from penetration experiments into concrete by ogive-nose steel projects[END_REF], [START_REF] Daudeville | Concrete structures under impact[END_REF], [START_REF] Smith | Numerical analysis of projectile penetration and perforation of plain and fiber reinforced concrete slabs[END_REF].

The quasi-static constitutive behavior of an ordinary concrete was extensively studied at 3SR laboratory thanks to triaxial compression tests performed with a large capacity press named Giga [START_REF] Vu | Strain measurements on porous concrete samples for triaxial compression and extension tests under very high confinement[END_REF]. Under high confinement (some hundreds MPa), the water saturation of concrete plays a major role [START_REF] Vu | Experimental analysis of concrete behavior under high confinement: effect of the saturation ratio[END_REF]. On the one hand, the hydrostatic behavior of wet or saturated concrete clearly becomes stiffer than that of dry concrete. On the other hand, the shear strength of wet or saturated concrete seems limited to a maximum value independent on the confining pressure, while the shear strength of dry concrete increases almost linearly with confining pressure. The limit shear strength value is directly correlated with the water saturation ratio of concrete [START_REF] Vu | Experimental analysis of concrete behavior under high confinement: effect of the saturation ratio[END_REF]. It is also worth noting that all these tests are performed under quasistatic undrained conditions similarly to impact condition; as the water has not enough time to migrate during impact.

Because protective concrete structures are generally massive, they may have a core partially or fully water saturated whereas the skin is dry. Accounting the effect of water saturation ratio is then particularly relevant for massive concrete structures submitted to very high stress level.

The PRM coupled model [START_REF] Pontiroli | Prediction concrete behaviour from quasi-static loading to hypervelocity impact An overview of the PRM model[END_REF], [START_REF] Vu | Influence of the Saturation Ratio on Concrete Behavior under Triaxial Compressive Loading[END_REF] is a finite element (FE) model that uses the effective stress concept developed by Mariotti et al. [START_REF] Mariotti | A numerical approach for partially saturated geomaterials under shock[END_REF] for wet geomaterials under high stress level. Thus, such a model takes into account the saturation ratio, but finite element models are not well adapted for large discontinuity and perforation simulations. The discrete element method (DEM) has several advantages. It can easily represent discontinuities caused by cracking or fragmentation and reproduce the macroscopic behavior of concrete. Several DE models were developed to reproduce the behavior of granular and cohesive materials [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], [START_REF] Cusatis | a) Lattice Discrete Particle Model (LDPM) for failure behavior of concrete[END_REF], [START_REF] Cusatis | Lattice Discrete Particle Model (LDPM) for failure behavior of concrete. II: Calibration and validation[END_REF], [START_REF] Donze | Numerical study of compressive behavior of concrete at high strain rates[END_REF], [START_REF] Hentz | Identification and Validation of a Discrete Element Model for Concrete[END_REF], [START_REF] Poinard | Compression triaxial behavior of concrete: the role of the mesostructure by analysis of X-ray tomographic images[END_REF], [START_REF] Shiu | Penetration prediction of missiles with different nose shapes by the discrete element numerical approach[END_REF], [START_REF] Shiu | Influence of the reinforcement on penetration and perforation of concrete targets: a discrete element analysis[END_REF], [START_REF] Suchorzewski | Discrete element method simulations of fracture in concrete under uniaxial compression based on its real internal structure[END_REF], but none of them accounts for the effects of saturation ratio.

This paper presents a new local constitutive behavior for discrete element modelling of concrete structures that attempts to take into account the influence of free water into the porous concrete material under very high quasistatic stress level. The proposed model is based on the use of spherical DEs, it was developed for concrete that is a cohesive material. Thus, contrarily to models developed for granular materials [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], cohesive interactions are considered between elements that are not in contact. As long as the material is cohesive, the proposed approach is similar to lattice models [START_REF] Schlangen | Fracture simulations of concrete using lattice models : Computational aspects[END_REF], [START_REF] Lilliu | 3D lattice type fracture model for concrete[END_REF], [START_REF] Yip | Irregular lattice models of fracture of multiphase particulate materials[END_REF]. Once fracture and fragmentation occurs, contact interactions may be created and easily handled with spherical DEs. The DE model is implanted is the open source code YADE [START_REF] Kozicki | YADE OPEN DEM: an open source software using a discrete element method to simulate granular material[END_REF] with the aim of modelling the behavior of concrete at the macroscopic scale to deal with impact problems on structures [START_REF] Yankelevsky | Local response of concrete slabs to low velocity missile impact[END_REF], [START_REF] Yankelevsky | Uniaxial compressive strength effect on high velocity penetration into thick NSC and HSC targets[END_REF]. Thus the internal structure of concrete (aggregates, cement paste, pores) is not described like in [START_REF] Kozicki | Discrete simulations of a triaxial compression test for sand by DEM[END_REF], [START_REF] Nitka | Modelling of concrete behaviour in uniaxial compression and tension with DEM[END_REF] and [START_REF] Suchorzewski | Discrete element method simulations of fracture in concrete under uniaxial compression based on its real internal structure[END_REF]. In this study, the three dimensional DEM constitutive model implemented in YADE will be described. The model calibration and validation is based on triaxial compression tests performed with the large capacity press Giga of 3SR laboratory at the university Grenoble Alpes. Numerical simulations of tests performed on ordinary concrete samples at different saturation ratios including triaxial and oedometric tests will be discussed.

Discrete element background

The DEs are rigid spheres of different radii. Each of them has a mass and a rotational inertia. Note that DEs do not represent aggregates, interaction laws between DEs are chosen to represent the macroscopic behavior of concrete as well as discontinuities (cracks) that may appear in the medium. Cohesive interactions are spring-like connections in extension, shear and rotation. The interaction force F represents the action between two elements a and b is given in eq.( 1). The interaction stiffness K and displacement U are decomposed to K n , K s and U n , U s acting in the normal and shear direction respectively.

F = KU (1)
Using the constitutive model of each interaction, the numerical model solves the equation of motion of the spheres assembly. The new displacement, velocity and acceleration for each discrete element are calculated by solving the dynamic equilibrium equations based on an explicit time integration scheme conditionally stable. The condition of stability applied for motion equations is inspired from [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] and [START_REF] Potyondy | Stiffness matrix at a contact between two clumps[END_REF].

Constitutive model description

At macroscopic scale, the concrete can be considered as a homogeneous isotropic material. Under confined compression and increasing mean stress levels, concrete first follows a linear behavior, and then its tangent stiffness decreases due to the damage of cement matrix, which leads to porosity closure. During the porosity closure, the stiffness of the material increases and tends to the elastic stiffness of the fully consolidated material [START_REF] Malecot | Strength and damage of concrete under high triaxial loading[END_REF]. This behavior at the macroscopic scale can be reproduced through simple interaction laws between the discrete elements. Shiu et al. [START_REF] Shiu | Compaction process in concrete during missile impact: a DEM analysis[END_REF] developed an interaction law featuring an elastic quasi-brittle behavior in tension, and a three linear behavior in compression, but this model does not take into account the effect of the concrete saturation ratio. In this new model, the effect of free water is taken into account by introducing a dependency between water saturation ratio and inelastic deformations due to pore closure. The porosity closure is defined as a local variable at the interaction level using the decrease of the distance between DEs. A maximum shear stress criterion is also introduced to take into account the correlation between limit shear strength and saturation ratio. The proposed model is calibrated by means of test results performed on a reference ordinary concrete named R30A7. Gabet et al. [START_REF] Gabet | Triaxial behavior of concrete under high stresses: Influence of the loading path on compaction and limit states[END_REF], Vu et al. [START_REF] Vu | Effect of coarse aggregate size and cement paste volume on concrete behavior under high triaxial compression loading[END_REF] and Piotrowska et al. [START_REF] Piotrowska | Experimental investigation of the effect of coarse aggregate shape and composition on concrete triaxial behavior[END_REF], studied the triaxial behavior of dry, wet or fully saturated R30A7 concrete under a confinement stress up to 600 MPa. Experimental results obtained previously emphasized on the fact that the presence of free water in pores has an effect only when the free porosity (not occupied with water) is closed. According to our assumptions, if the boundary allows the water outflow, then the effect of water vanishes and the constitutive behavior of concrete will be that of dry concrete. Note that differed phenomena due to water migration are not taken into account in the model since, due to the low loading time, no water migration is assumed.

Interaction laws

In this model, interactions between DE neighbors that are not in contact are allowed by introducing an interaction radius coefficient (I r ≥ 1). Two types of interactions are defined: initial links and residual contacts. The value of I r ,which concerns only initial links, influences directly the number of interactions per sphere at the beginning of the calculation (eq.( 2)). Rousseau et al. [28] showed that in order to achieve a realistic isotropic elastic behavior of concrete, this number should be around 12. The interaction radius coefficient I r will be adjusted accordingly. As calculation is proceeded and due to external loading, changes in DEs arrangement occur and new contact interactions may be created while others vanish.

Cohesive interactions

Cohesive interactions are created between DEs within the same range at the beginning of the simulation. For instance, a link is created between elements a, and b of radius R a and R b respectively if:

I r (R a + R b ) ≥ D 0 ab (2)
Where D 0 ab is the initial distance between the centroids of elements a and b. D 0 ab is also the reference distance used later to convert displacement to dimensionless strain. If C 0 a and C 0 b are initial centroids position of elements a and b respectively, then:

D 0 ab = |C 0 a -C 0 b | (3)
Tensile stress appears when the current distance between elements a and b:

D ab = |C a -C b | is larger than the initial distance D 0
ab because the interaction is cohesive.

Contact interaction

Contact interactions can be created between DEs after the beginning of the simulation, either when a cohesive Link interaction is broken or when two DEs that were not in interaction are in contact (R a + R b ) ≥ D 0 ab .

Elastic interaction law

Interaction force F can be decomposed in a normal force F n and a shear force F s :

F n = K n U n (4)
F s = K s U s (5)
U n and U s are respectively normal and shear displacements. K n and K s are respectively normal and shear stiffnesses, both expressed using E an equivalent Young modulus and G an equivalent shear modulus for the link:

K n = EA eq D 0 ab (6)
K s = GA eq D 0 ab (7)
The macroscopic elastic properties, here Young modulus E (GPa) and shear modulus G (GPa), are thus considered to be an input parameters. A eq is the interaction cross-section, and is defined as follows:

A eq = π × min(R a , R b ) 2 (8)
Thus, the forces could be defined from equivalent stresses as follows:

F n = σ n A eq (9)
F s = σ s A eq [START_REF] Frangin | On the use of combined finite/discrete element method for impacted concrete structures[END_REF] with,

σ n = Eε n (11) σ s = Gε s (12)
where ε n and ε s are the equivalent normal and shear strain respectively

ε = ε n ε s (13)
A fictitious interaction point P c is defined at a mid-distance from the two elements and a vector n which is the normal to the contact plane that passes through P c (always perpendicular to the contact plane). U n is aligned with n. The shear displacement U s must be expressed in the global coordinates while satisfying the condition U s ⊥ n. This is due to the lack of an interaction-local coordinate system. Indeed, each spherical contact only defines the contact normal and the contact plane (Figure 1). As large strains are expected, logarithmic strains, also equivalent to the sum of incremental strains are used (eq. ( 14)). The strain tends to -∞ if two centers of spheres approach from one another, which avoid them to penetrate through each other.

ε = ε n ε s = log( U n D 0 ) ∑ δU s D 0 (14)
For torque computation, the following equations are used:

M b = ∑ δ θ b K b (15) M t = ∑ δ θ t K t (16)
δ θ b and δ θ t are the incremental bending and twisting rotation angle respectively. K b and K t bending and twisting stiffnesses that are estimated by considering that a circular beam, with a section A eq , connects the two centroids of spheres a and b (eq. ( 8)). To limit the torque, two dimensionless parameters α b , α t ∈ [0, 1] are used such that:

K b = α b K s A eq (17)
K t = α t K s A eq (18)

Nonlinear interaction behavior

As shown in equations (eq. ( 11) and eq. ( 12)), σ n and σ s are related to ε n and ε s respectively. This relation between stress and strain at the link scale characterizes the interaction behavior in both tension and compression.

Tension

Beyond the elastic limit ε 0 , a damage behavior in tension is applied. Normal stress σ n is formulated as follows:

σ n = [1 -ω(κ)H(κ -ε 0 )]Eε n ( 19
)
κ is the maximum normal strain (κ = max(ε n )), ω(κ) is the damage evolution function, and H is a Heaviside function that deactivates damage effect if κ < ε 0 . The damage evolution function ω is described as follows:

ω(κ) = 1 -ε 0 κ 1 -ε 0 ε f if ε 0 ≤ κ ≤ ε f ( 20
)
ε f is the maximum strain corresponding to the maximum damage (ω = 1 when ε n = ε f ). Beyond this point, the link between the spheres is deleted and a new contact interaction is created only if these two spheres touch again. Figure 2 shows the evolution of the damage and the tensile stress in the normal direction.

Compression

The interaction model between DEs is a beam-like model whose constitutive behavior is phenomenological. It means that the local constitutive behavior is inspired from observations at the macro scale, including compaction (pore closure). Thus, under normal compression the stiffness varies between the initial elastic stiffness K n (directly linked to the Young modulus) and the consolidated material stiffness k h x k n (k h ¿1 is an non-dimensional parameter). k p and k h values are calibrated by means of simulation of confined compression tests. Thus, in compression the behavior phases are defined as follow (see Fig. 3):

Phase [A B]: linear zone of the link, characterized by the initial elastic stiffness of the material. ε el is the elastic limit strain in compression. Phase [B C]: compaction zone, the elastic stiffness of the link varies linearly between the initial elastic stiffness and the consolidated material stiffness. Let us note that at this point porosity is inherently defined at the link scale since ε pl is the maximum deformation corresponding to a complete porosity closure. k p is the strain hardening modulus introduced as:

σ n = Eε n for ε n < ε el (21)
σ n = E k p (ε n -ε el ) + σ el max for ε el < ε n < ε pl (22)
Phase [C D]: consolidated material zone. Beyond this point, the link is fully compacted and the response is incrementally elastic again as described in (eq. 23). Figure 3 also shows the cyclic behavior of a link in the normal direction (interaction level).

σ n = Ek h (ε n -ε pl ) + σ pl max for ε pl < ε n (23)

Shear

As explained in (sec. 3.2) the elastic shear stress σ s can be defined as shown in (eq. ( 12)). The shear response is then limited by the maximum shear stress surface (Figure 4). The shear stress criterion σ s max (see eq. ( 24)) is characterized by

σ s max = (1 -ω)C 0 + λC 0 σ n tan(Φ) σ n tan(Φ) + (λ + 1 -ω)F t (24)
ω is the damage function, F t the undamaged tensile strength (F t =Eε 0 ) and λ a dimensionless parameter that allows controlling the maximum shear stress asymptote.

Effect of saturation ratio (Sr) on the compressive behavior

The effect of free water on concrete mechanical behavior has been widely studied under high or moderate strain rates [START_REF] Rossi | The dynamic behaviour of concrete: influence of free water[END_REF], [START_REF] Zhou | Influence of free water content on the compressive mechanical behaviour of cement mortar under high strain rate[END_REF], [START_REF] Wu | Influence of strain rate and water content on mechanical behavior of dam concrete[END_REF], [START_REF] Sun | Effect of pore liquid viscosity on the dynamic compressive properties of concrete[END_REF]. Due to a lack of experimental data, there is no available discrete element model dealing with free water influence on the mechanical behaviour of concrete under high stress level. The proposed model is the first attempt aiming at filling this gap. Saturation ratio is introduced as a local variable at the interaction level between spheres and its evolution is related to the porosity closure during the compaction phase described by (eq. ( 22)). In this paper, an initial homogeneous, saturation ratio distribution is assumed in the sample however the structural heterogeneity, or a heterogeneity based on the DE size or position, could easily be taken into account. Note that, the effect of saturation ratio is taken into account for compressive stress only. As it was explained in (sec. 3.3.2), for a dry link (S r = 0%), the consolidation point is defined accordingly to normal strain, and it corresponds to a complete porosity closure (ε n = ε pl ). Then, the general idea of the model is to assure that a smaller deformation is needed to close a partially saturated porosity as a fraction of it is already occupied by water. The higher the saturation ratio is, the earlier the consolidation point will be reached and the effect of water will begin to appear on the mechanical behavior. Thus, the new consolidation point ε n,plC is defined as a linear function of S r which is the saturation ratio of the link:

ε n,plC = ε el + (ε pl -ε el )(1 -S r ) (25)
From equation (eq. ( 25)), it is clear that if the link is dry (S r = 0), ε n,plC is equal to ε pl while if the link is completely saturated (S r = 1) and then ε n,plC is equal to ε el . A kind of effective stress concept is then introduced to take into account the water contribution in the total stress at the link scale using Biot coefficient assumed to be, at the first order, equal to the porosity ϕ. The variation of the porosity is defined accordingly to the volumetric strain ε v which is defined from the normal strain of the unidimensional link as ε v = 3ε n and to ϕ i which is the initial porosity of the sample. Thus:

σ n total = σ n + 3ϕσ n water (26)
One can remark from (eq. ( 26)) that the water contribution will increase with the porosity and with the saturation ratio. Note that a Mie -Gruneisen like equation of state is used to take into account water compressibility and to compute the water contribution in the total stress [START_REF] Davis | Thermodynamic response of Mie-Grüneisen materials at high pressures[END_REF] (eq. ( 27)):

σ water = ρ 0 C 2 w0 (ε v -ε v,ps ) (1 -s(ε v -ε v,ps )) 2 [1 - Γ 0 (ε v -ε v,ps ) 2 ] ( 27
)
Where ε v,ps is the volumetric strain at the consolidation point (ε v,ps = 3ε n,plc ), C w0 is the sound-wave velocity, ρ 0 is the mass density, s and Γ 0 are two Mie-Gruneisen coefficients. The interstitial water pore pressure was recently investigated and measured experimentally by [START_REF] Accary | Design and evaluation of a deformable sensor for interstitial pore pressure measurement in concrete under very high stress level[END_REF] showing that it might reach several hundreds of megapascal under high confinement pressure. Figure 5 shows the normal compressive behavior of a dry link compared to that of a wet link (S r = 44%) and the corresponding water pressure. The consolidation point (CP) is reached earlier in the wet link hence the water contribution is added. In the dry link, reaching the consolidation point does not add any water contribution.

Fig. 5 Effect of the water contribution on the normal compressive stress-strain law at the link scale for S r = 44% compare to S r = 0%

Effect of saturation ratio on the maximum shear stress

In order to control the plastic sliding threshold of interactions, a coefficient λ was introduced in eq. ( 24). This coefficient control the maximal limit shear strength of a link (λC 0 ) once the material is fully consolidated as observed experimentally by [START_REF] Vu | Experimental analysis of concrete behavior under high confinement: effect of the saturation ratio[END_REF]. Two parameters λ 0 and λ 100 , are calibrated for a dry and saturated links respectively and introduced to express λ as a linear function of the saturation ratio:

λ = (λ 0 -λ 100 )(1 -S r ) + λ 100 (28)
Thus, when λ = λ 0 then ε n,plC = ε pl and when λ = λ 100 then ε n,plC = ε el . Figure 6 shows the dependency of the criterion on the saturation ratio. We can observe that the higher the saturation ratio is, the lower is the shear stress limit.

Numerical tests preparation and monitoring

The calibration is done by running simulations on sufficiently large numerical specimens to give continuum-like behavior. Samples, made of spheres, have to be isotropic to ensure a homogeneous interaction distribution in the sample and to prevent the forces from being in privileged directions [28], [START_REF] Frangin | On the use of combined finite/discrete element method for impacted concrete structures[END_REF]. The isotropy of the numerical medium is verified by plotting the cumulative orientation distribution of interactions (Figure 7).

The histogram of the size distribution obtained for a specimen is shown in Figure 8(a). The size distribution goes from 2 to 4 mm and is homogeneously distributed between these extremum. The numerical compacity of the sample is about 0.6 and the interaction radius I r is chosen in order to have an average of 12 interactions per element (I r = 1.5). Figure 8(b) shows the numerical sample. It is a cuboid-shaped packing (70mm×70mm×140mm) of approximately 4000 spheres. It was checked that the shape of the numerical specimen does not affect the results since the stress state is homogeneous across the specimen. ence concrete (unconfined compressive strength of 34 MPa, 7 cm slump, w/c = 0.64, 12% porosity accessible to water, see [START_REF] Gabet | Comportement triaxial du béton sous fortes contraintes: Influence du trajet de chargement[END_REF] for more details).

Model calibration and comparison with experiments

Uniaxial compression and tensile tests

Uniaxial tension or compression tests allow determining elastic parameters E, ν, initial cohesion C 0 , limit elastic strain in tension ε 0 and the limit failure strain ε f . The displacement is applied on boundary particles symmetrically on both ends of the specimen while restraining their other degrees of freedom. Average axial stress is obtained by averaging forces on both boundary particles divided by the specimen crosssection (assumed constant during the test). Figure 9 shows the numerical stress-strain curve in uniaxial compression, it is compared a cyclic test result [START_REF] Poinard | Damage of concrete in a very high stress state: Experimental investigation[END_REF]. The stress-strain curve in tension is not available, only the tensile strength (3.2 MPa) could be obtained experimentally. Figure 10 shows a comparison between numerical and experimental failure patterns in the sample for uniaxial tension and compression tests respectively. In tension, a single crack appears; the stress state is homogeneous far from the ends of the specimen so its location depends on the statistical distribution of spheres into the sample. In compression, the ends are fixed then a typical cone failure pattern can be observed in the central part of the specimen. Benniou [START_REF] Benniou | Modélisation par éléments discrets du comportement des matériaux ci-mentaires sous impact sévère prise en compte du taux de saturation[END_REF] has shown in his PhD that there is no mesh dependency for a moderate variation of DE sizes; this result was also shown in [START_REF] Potapov | Mixed DEM/FEM modeling of advanced damage in reinforced concrete structures[END_REF]. Note that Benniou [START_REF] Benniou | Modélisation par éléments discrets du comportement des matériaux ci-mentaires sous impact sévère prise en compte du taux de saturation[END_REF] has also verified that uniaxial compresssion and tensile test results do not depend on the initial saturation ratio of concrete samples.

Triaxial compression tests at high level of confinement (600 MPa)

Tests at 600 MPa of confinement are used to calibrate parameters given in Table 2. It has been checked experimentally that the stress state during confined compression tests Poisson ratio 0.2 ε 0 limit elastic strain in tension 1e -4 ε f limit failure strain 20ε 0 Φ c (radians) contact friction angle 0.8 C 0 (MPa) shear cohesion 4

is homogeneous in the circular section specimens [START_REF] Vu | Strain measurements on porous concrete samples for triaxial compression and extension tests under very high confinement[END_REF]. In a sake of simplicity, the DE specimen has a parallelepipedic shape with 6 contacts surfaces. This choice allows prescribing simple boundary conditions by means of six rigid walls; it has no influence since the obtained stress state is homogeneous. Triaxial calibration tests are run on dry (S r = 0%) and completely saturated (S r = 100%) numerical sample. These tests are conducted in two steps. The first phase is the hydrostatic compression, for which all the wall displacements are the same with a controlled velocity until the target confinement pressure P c is reached σ x = σ y = σ z = σ m = -P c . This step allows the calibration of non-linear parameters k p , k h , ε el and ε pl , (see Table 2, left part). The second phase is the deviatoric compression, for which the axial displacement is controlled on the top and bottom wall whereas all four other walls are monitored by their normal stress so that constant σ y = σ z = -P c . This step allows the calibration of the two last parameters λ 100 and λ 0 , (see Table 2). [START_REF] Malecot | Influence of free water on concrete triaxial behavior: The effect of porosity[END_REF] and of numerical simulations obtained for an hydrostatic and deviatoric compression phase respectively for (S r = 0%) and (S r = 100%). Numerical simulations show a good agreement with experimental tests and the capability of the model to reproduce results at the macroscopic scale. The effect of saturation ratio is pronounced for hydrostatic phase as shown in Figure 11, a relative difference of about 25% between the volumetric strains of dry and saturated samples at a mean stress of 600 MPa is observed. For deviatoric phase, the saturation ratio notably affects the shear limit as shown in Figure 12. At 600 MPa of confinement pressure, the shear resistance drops from 800 MPa for a dry sample to 250 MP for a saturated sample which is only due to the presence of free water in the saturated sample. Note on Figure 12 that the comparison between the experimental and numerical results is not perfect for dry concrete and an axial strain greater than 8%; it is a minor drawback of the DE model that the authors did not try to fix because the maximum deviatoric stress varies slightly for different tests performed on dry concrete [START_REF] Gabet | Triaxial behavior of concrete under high stresses: Influence of the loading path on compaction and limit states[END_REF]. Figure 14 shows simulation results of the deviatoric phase of triaxial tests at 100 MPa and 200 MPa compared to experimental results obtained by [START_REF] Malecot | Influence of free water on concrete triaxial behavior: The effect of porosity[END_REF] for different saturation ratios. Triaxial tests at moderate confining pressures are difficult to simulate because these mean stress levels correspond to the transition from brittle to ductile behavior, which means, to have in the same time a significant effect of both damage and plasticity. However, Figure 14 shows that the simulation results are in good agreement with the experimental ones.

From the one hand, Fig. 12 and 14 indicate that the deviatoric behavior of dry concrete is highly influenced by the confining pressure since as the confining pressure increases from 100 to 600 MPa the deviatoric stress increases from 250 to 800 MPa. From the other hand, for a given confinement pressure and for intermediate saturation ratios (Fig. 14(a) and 14(b)), the peak deviatoric stress lies between the dried and saturated limit cases and seems to be associated directly to the saturation ratio. Thus the limit stress decreases when the saturation ratio increases. Fig. 14(b)) shows nonsatisfying agreements for the 200 MPa confinement triaxial tests with initial saturation ratio Sr equal to 78% and 100%; all other results are very good. The authors explain these results by a possible overestimation of S r [START_REF] Vu | Experimental analysis of concrete behavior under high confinement: effect of the saturation ratio[END_REF], the two samples might have slightly dried just before the tests.

Oedometric tests

The static quasi oedometric compression (QOC) test allows testing the mechanical response of the material under a quasiuniaxial strain loading path while applying a passive confinement. The difference between hydrostatic and quasi oe-dometric tests lies mainly in the volumetric strains reached for the same mean stress. During this test, a cylindrical specimen, tightly enclosed in a confinement vessel, is axially compressed by means of high strength compression plugs as stated by [START_REF] Burlion | Experimental analysis of compaction of concrete and mortar[END_REF]. To simulate the quasi oedometric test, the sample is placed between six rigid walls. The axial displacement of both top and bottom walls are controlled while blocking all four other walls to simulate a uniaxial confined compression test under uniaxial strain. Figure 15 shows predictive mean stress evolution obtained for different saturation ratios. An experimental result obtained on a dry sample is also shown. The behavior of the numerical sample for S r = 0% is very close to the experimental behavior. The effect of free water is also shown by means of numerical simulations. As expected, the results reveal that when the sample saturation increases, its oedometric stiffness also increases.

Conclusion

The DE model presented in this paper is dedicated to the modelling of concrete behavior at very high stress level. To take into account the influence of free water during porosity closure that is observed in experiments, a new behavior law was developed at the scale of links between DEs. This constitutive behavior is based on a dependency of inelastic deformations of concrete and maximum shear stress criterion with the saturation ratio. The results show that the model is able to well reproduce the behavior of concrete samples on a wide range of stress levels (triaxial compression with a confining pressure up to 600 MPa) and on a wide range of saturation ratios (S r = 0 to 100%). The identification of parameters is carried out thanks to tests performed on dry samples. The simulation results obtained for various loading paths and saturation ratios demonstrate the efficiency of the approach to take into account the influence of free water at very high stress level. As a perspective, simulations at a larger scale on concrete structures are feasible and will allow the evaluation of saturation ratio effects at a structural scale. The effect of saturation ratio on the behavior of concrete structures submitted to penetration or perforation tests will be evaluated thanks to this model.

Fig. 1

 1 Fig. 1 Interaction between two elements. Force and displacement components

Fig. 2 a

 2 Fig. 2 a Damage evolution ω vs. normal strain ε n in a cohesive inter- action. b Corresponding tensile stress σ n vs. ε n

Fig. 3

 3 Fig. 3 Cyclic loading of a single cohesive link. a In the normal direction. b Zoom of A region

Fig. 4

 4 Fig. 4 Maximum shear stress surface versus uniaxial stress at the link scale

Fig. 6

 6 Fig. 6 Maximum shear stress surface for different saturation ratios (S r = 0%, S r = 44%, S r = 78% and S r = 100%) at the link scale

Fig. 7 Fig. 8

 78 Fig. 7 Initial distribution of the interaction links' orientations in the numerical sample

Fig. 9

 9 Fig.9Numerical stress-strain curve in uniaxial compression and comparison with experimental curve of a cyclic test[START_REF] Poinard | Damage of concrete in a very high stress state: Experimental investigation[END_REF]

Fig. 10

 10 Fig. 10 Comparison between numerical and experimental damage after uniaxial test: blue elements links are undamaged, red elements links are completely damaged. a,b Tension [11]; c,b Compression test[START_REF] Poinard | Damage of concrete in a very high stress state: Experimental investigation[END_REF]

Fig. 11

 11 Fig. 11 Hydrostatic compression test; mean stress vs. volumetric strain: comparison between experiment and modelling for S r = 0% and S r = 100%

Fig. 13

 13 Fig.13Simulation of hydrostatic compression test; mean stress vs. volumetric strain for different saturation ratios (S r = 0%, S r = 44%, S r = 78%, S r = 93% and S r = 100%)

Fig. 14

 14 Fig. 14 Deviatoric phase of triaxial tests, deviatoric stress vs. axial strain for different saturation ratios (S r = 0%, S r = 44%, S r = 78% and S r = 100%): comparison between experiment and modelling. a 100 MPa confinement. b 200 MPa of confinement

Fig. 15

 15 Fig. 15 Simulation of oedometric compression tests: mean stress vs. volumetric strain for different saturation ratios (Experiment Sr=0%; Numerical simulation S r = 0%, S r = 44%, S r = 78% and S r = 100%)

Table 1

 1 Model parameters calibration values from uniaxial tests

	Parameters	Physical meaning	Values
	E (GPa)	Young modulus	30
	ν

Table 2

 2 Model parameters calibration values from triaxial compression tests at 600 MPa

	Parameters Physical meaning	Values
	k p		strain hardening modulus	2
	k h		consolidated material parameter	1
	ε el		limit elastic strain in compression	20ε 0
	ε pl	maximum compaction strain	200ε 0
	λ 100	saturated sliding threshold coefficient 1
	λ 0	dry sliding threshold coefficient Figures 11 and 12 display the results of experimental
	tests carried out by

Acknowledgements This research were financially supported by CEA DAM, the Direction of military applications of the French Alternative Energies and Atomic Energy Commission. The authors would also like to thank Dr. Christophe Pontiroli for scientific advice.

Conflict of interest

The authors declare that they have no conflict of interest.