Three-dimensional simulations of ageing concrete structures using a multiphase model formulation

The durability of concrete structures is in no small degree determined by the quality and integrity of the concrete, where structural damages such as cracks negatively affect many of the functions of the structure. Often cracks are formed due to restrained thermal and hygral deformations, where the risk is exceptionally high during the early stages after casting. This study presents a hygro-thermo-chemo-mechanical model that accounts for phenomena such as hydration, external and internal drying, self-heating, creep, shrinkage and fracture. The model is derived as a modified version of a fully-coupled multiphase model recently proposed by Gasch et al. (Cem Concrete Res 116:202–216, 2019. https://doi.org/10.1016/j.cemconres.2018.09.009) and implemented in the Finite Element Method. Here the governing equations are simplified, and a more efficient solution method is proposed. These modifications are made with the intention to obtain a model more suitable for structural scale simulations. To validate the model, one of the end-restrained beams tested within the French research project CEOS.fr is analyzed. Laboratory data on the concrete is used to calibrate to model and recordings of ambient conditions makes it possible to define accurate boundary conditions. Results from the simulation are compared to measured temperatures and deformations from the first 60 days after casting and are found to generally be in good agreement. Compared to the fully-coupled model by Gasch et al. (2019), the modifications proposed in this study reduce the computational cost by a factor five; without any noticeable differences to the structural results.

Laboratory data on the concrete is used to calibrate to model and recordings of ambient conditions makes it possible to define accurate boundary conditions. Results from the simulation are compared to measured temperatures and deformations from the first 60 days after casting and are found to generally be in good agreement. Compared to the fully-coupled model by Gasch et al. (2019), the modifications proposed in this study reduce the computational cost by a factor five; without any noticeable differences to the structural results. For concrete structures in general, and infrastructure applications such as nuclear power, hydropower and bridges in particular, the quality of the concrete is essential in order to ensure a safe and durable structure. Early-age cracking due to restrained thermal and hygral deformations as well as poor hardening conditions not only affect the load carrying capacity of a structure greatly but also many other functions that the structure is required to fulfil. Such functions can include stiffness requirements for bridges and leaktightness for both nuclear and hydropower structures.
In addition, pre-existing cracks cause an increased risk and rate for other degradation processes such as frost action, leaching and alkali-aggregate reactions (AAR) when in contact with free water, thereby significantly reducing the expected service life of the structure. Although much progress has been made on this topic during the last decades, there is still a need to increase further the knowledge on how to prevent these, often unnecessary, structural damages due to restrained deformations, both in connection to casting but also during later stages due to for example environmental loads. The need exists for both experimental and theoretical studies on the material behaviour of concrete, and moreover to realize the knowledge gained from such studies into design aids and predictive models that are suitable for analyses on the structural scale, the latter being the focus of this study. The aim is to present and verify a numerical model useful for structural scale studies of early-age cracking, including effects such as cracking, creep, drying and thermal shrinkage.
Many predictive models have been presented and applied to the early-age behaviour of concrete structures in the literature. However, often these only consider the thermo-mechanical effects [1][2][3]. It is mainly during the last decade that models additionally accounting for hygral effects have been considered and also more refined models for the hardening process devised and coupled to other physical fields [4][5][6][7][8]. The development of models to also consider hygral effects for early-age concrete is in no small degree motivated by the increased usage of highperformance concrete, where self-desiccation and autogenous shrinkage during hydration is more pronounced than in conventional concrete. However, these effects are also noticeable and significant for conventional concrete. Hygral effects can be vital to consider when studying concrete structures subjected to realistic environmental conditions when for example considering drying shrinkage. A class of models that have proven fruitful in describing many of the processes essential for early-age and mature concrete is the adoption of multiphase porous media theory to early-age concrete by Gawin et al. [9,10], together with later developments on the same framework, e.g. [11][12][13][14]. In this class of models, concrete is described as a three-phase porous media and the governing equations employed to describe its hygro-thermochemo-mechanical (HTCM) properties are formulated based on a mechanistic approach through the averaging theorems outlined by, for example, the Thermodynamically Constrained Averaging Theory (TCAT) [15]. While the original model by Gawin et al. [9,10] presents a fully-coupled approach to the HTCM behaviour and is mainly validated using laboratory scale examples, it has later been modified to be more suitable for structural applications by, for example, Scium'e et al. [14] through the decoupling of the hygro-thermo-chemical problem and the mechanical problem. Sciumé et al. applied and verified this modified model to repair problems, and also extended the original model to consider concrete cracking through damage mechanics. Similar modelling approaches have, furthermore, been successfully applied to study other applications related to concrete structures, for example, exposure to fire [16,17], frost action [18,19] and AAR [20].
This paper presents a HTCM model based on a multiphase porous media approach intended for predictive analyses of concrete structures from casting and during the later stages of its service life. It is based on the model recently proposed by Gasch et al. [11], in which a fully-coupled approach is adopted to the HCTM problem. As done by Sciumé et al. [14], the current study presents a decoupling of the governing equations and the suppression of several weak physical couplings, but still maintains the same general format of the governing equations and other properties of the model. While most constitutive equations employed also remains the same as in [11], some new and modified ones are also provided. Furthermore, a new solution algorithm is presented to replace the fully-coupled solution method employed in [11], in order to arrive at a model that is more computationally efficient and thus more suitable for large-scale applications. To analyse the applicability of the proposed model to describe cracking due to restrained deformations in a structural scale application, one of the end-restrained beams tested within the French research project CEOS.fr [21] is modelled. This example includes both the early-age behaviour with cement hydration, self-desiccation and self-heating, as well as actual environmental conditions including hourly variations in temperature, relative humidity and solar radiation.
First, the general mathematical framework of the multiphase model is reviewed, and its numerical implementation discussed. It is followed by the application of the model to the CEOS.fr beam, including calibration of concrete material properties, model set-up and a discussion on the simulation results in comparison to the measured data. Finally, some general conclusions are given.
2 Multiphase model for early-age concrete Mechanics of porous media is today a major branch of continuum mechanics and of interest for many engineering applications within for example geomechanics and biomechanics. When formulating the governing equations that describe the behaviour of the porous media, either a phenomenological or a mechanistic approach can be adopted [22]. In the phenomenological approach, the behaviour of the multi-phase system is derived directly on the macroscopic scale through observation of experiments and by fitting differential equations to these; not necessarily considering all processes involved. In contrast, a mechanistic model aims at a true physical description of the phenomena causing the observed behaviour of the material. For a porous media such as concrete, this means working with multi-scale considerations and up-scaling procedures. In this work, the latter approach is used, with the Thermodynamically Constrained Averaging Theory (TCAT) developed by Gray and Miller [15,23] as a starting point for deriving the conservation equations for early-age and mature concrete. The used model builds on the recently presented model by Gasch et al. [11] in which concrete is considered as a three-phase material consisting of a solid, a wetting and a non-wetting phase. Compared to other similar models in the literature, such as the models by Gawin et al [9,10], Jefferson et al. [12], Sciumé et al. [14], the solid phase is by Gasch et al [11] decomposed into several species, including the structural water held in gel pores (\3 nm), aggregates, unhydrated cement and hydration products. The same assumption is used in this study. However, the proposed model is modified to be more suitable for large-scale structural applications. As made by for example Sciumé et al. [14], this is done by introducing the following hypotheses: 1. The velocity of the solid phase is negligible compared to those of the two fluid phases. 2. The effects of gravity on fluid flow is negligible.
3. All species are considered incompressible. 4. The hygro-thermo-chemical phenomena does not depend on mechanical part of the problem.
These hypotheses allow us to decouple the hygrothermo-chemical part from the mechanical part of the model, which significantly reduces the computational cost and complexity of the model. It should be mentioned that Sciumé et al. [14] also assumes that the advective of heat transport is negligible, which is not done in the current study. Although similar to the equations presented by Gasch et al. [11], the final form of the governing equations of the multiphase model is presented in the following to highlight the differences. However, most constitutive and physical relationships are for brevity only presented in Appendix A. For more background on the model equations, the reader is referred to reference [11], as well as [9,10,14,24,25] for more general descriptions of multiphase modelling of concrete.

Hygro-thermo-chemical behaviour
The hygro-thermo-chemical (HTC) behaviour of the concrete is given by evolution equations to describe the extent of cement hydration together with the conservation of mass and energy of all species in all phases. Compared to the more detailed model by Gasch et al. [11] it can be noticed that all terms related to the movement of the solid skeleton are omitted from the presented equations. The dependent variables of the HTC part of the model is the solid volume fraction s , the capillary pressure p c , the gas pressure p g and the temperature T. Furthermore, two state variables are used to track the extent of cement hydration through the hydration degree n C and an ageing variable j that is used to model the evolution of several material properties.

Cement hydration and ageing
The chemical model to describe the hydration process employed herein follows the work of Ulm and Coussy [26,27], using a normalized description of the reaction extent through n C . An Arrhenius type evolution law can be formulated for the reaction rate following Ulm and Coussy [27] and including later developments [4,28], such that The normalized affinity A C controls the reaction rate and is here defined using the formula proposed by Sciumé et al. [14], see Table 3. The first term on the right-hand side of Eq. (1) is introduced to account for the influence of the moisture state on the cement reaction in a simplified manner, where s w is the capillary saturation of the porous media [11]. The evolution of many material properties of concrete, such as strength, depends on the characteristics of the microstructure and how it evolves with time, and can thus be difficult to relate directly to the extent of hydration. To consider also how the curing temperature effects this, the model suggested by Cervera et al. [28] is adopted in which the evolution of an ageing degree j is introduced on the form where the first term accounts for the effect of different curing temperatures. It is done by prescribing a maximum temperature T T at which hydration can take place (% 100 C) [28]. The ageing degree can then be used to define the evolution of mechanical properties such as strength f AE and fracture energy G AE f , and also transport properties such as the intrinsic permeability k. For example, the mechanical strengths are given as where f AE 1 is the strength at complete hydration for a concrete cured at a constant temperature T ref . Similarly, G AE f is assumed proportional to ffiffiffi j p , see [11,28].

Conservation of mass
Using the characterization of concrete as a porous media by Gasch et al. [11] means that a total of seven mass balances needs to be considered, one for each species in each phase. However, to reduce the number of model equations these are grouped such that only three mass balances are included in the final model: one for the solid species A; C; H ð Þ , one for water W and one for dry air D. First, the mass balance of all solid species is utilized to obtain an evolution equation for the solid phase volume fraction which is defined by the rate of chemically bound water r Ws and the composition of the solid phase expressed by Notice that all volume occupied by each species is temperature dependent and that ia evolves with the extent of hydration, see Table 3. Furthermore, since the capillary porosity ¼ 1 À s , it follows that o s =ot ¼ Ào=ot, which is used in several of the following equations.
Water is in the framework proposed by Gasch et al. [11] present in all three phases of the porous media in contrast to most other similar models. This fact is of course reflected in the mass balance formulated for the total water content of the material, which takes the form o ot where the second term accounts for the water held in the gel pores; this term is normally not present and often lumped together with the capillary water. A consequence of this is that the gel water does not contribute to the fluid flow in Eq. (6), which is in line with the argumentation by Ulm et al. [29]. Also, the term s r Ws originates from the mass balance of gel water, meaning that it describes how the cement only reacts with water in the gel pores. All other terms are standard and present in most multiphase models of concrete, see for example Lewis and Schrefler [24]. Lastly, the mass balance of dry air is given on a standard form as o ot

Conservation of energy
The conservation of energy is formulated for the concrete as a whole by summation over all species and phases. Furthermore, as usually done for concrete, it is assumed that all phases are in local thermodynamic equilibrium and that velocities are considered small and the time scale of interest large (days). From the latter assumption, it follows that any contributions from viscous dissipation and mechanical work can be omitted [11,22]. Then by expressing the energy balance in terms of enthalpy, we obtain a model equation on a similar form to the classical heat equation:

Constitutive relationships of the hygro-thermochemical problem
Most constitutive and physical relationships that are necessary are summarized in Table 3 and described in detail by Gasch et al. [11]. To define the capillary pressure that is used as a dependent variable of the model, the Kelvin equation is used where also p Wg =p Wg sat is identified as the relative humidity u. Moreover, the capillary pressure can also be defined as the pressure difference between the nonwetting and wetting phases, i.e.
which is derived from the second law of thermodynamics given a state of local equilibrium, as shown by for example Pesavento et al. [30].

Initial and boundary conditions of the hygrothermo-chemical problem
To complete the HTC problem formulation, the governing equations have to be complemented by initial and boundary conditions. Initial conditions are specified for each dependent variable over the entire model domain X such that For the three physical fields, p c 0 , p g 0 and T 0 are typically identified from ambient data or measurements, while s 0 can be calculated from the properties of the concrete mixture, e.g. through the water-to-cement ratio w=c.
Additionally, boundary conditions are needed for the governing Eqs. (6)-(8) that contain spatial derivatives, either by directly prescribing the value of their respective dependent variable on surface C 1 , such that or the boundary flux. Since the proposed model is intended for concrete structures located in an outdoor environment, a mixed boundary condition is often appropriate to define the boundary flux. Hence, on the surface C 2 , the fluxes are given by When exposed to realistic environmental conditions, the surface flux of water through C c 2 is in Eq. (14) controlled by the ambient concentration of water vapour q Wg amb and the surface transfer coefficient f W . Similarly, for the heat flux through the surface C T 2 , both convective and radiative contributions are considered in such a case. The convective heat transfer is controlled by the ambient temperature T amb and the surface coefficient f T , which is calculated as a function of wind speed and the ambient surface film temperature [31]. For the radiative heat flux, the total irradiation G rad on the surface is calculated from the sum of the ambient and the solar radiation.

Mechanical behaviour
The mechanical behaviour of the concrete is governed by the conservation of linear momentum of all species and phases. As will be shown, it also depends on the current state of the HTC properties, which thus serve as an input to the mechanical part of the model. However, as already pointed out, no coupling in the opposite direction is considered. The dependent variable of the mechanical part of the model is the displacements d. In addition, to account for the complex behaviour of early-age concrete a number of state variables are utilized for the computation of damage and inelastic strains, see Sect. 4.

Conservation of linear momentum
The conservation of linear momentum is formulated for the concrete as a whole by summation over all species and phases. As for the HTC part, it is assumed that velocities are small and the time scale of interest large (days), which means that inertial forces and forces due to mass exchanges can be neglected [22]. Hence, the linear momentum balance takes the quasistatic form wheret is the damaged total stress tensor introduced to account for the fracture of the solid skeleton. In continuum damage mechanics literature,t is often referred to as the nominal stress tensor or simply the macroscopic stress tensor. The total density q is calculated as

Constitutive relationships of the mechanical problem
Most constitutive and physical relationships that are necessary are summarized in Table 4 and described in detail by Gasch et al. [11]. To clarify the phenomena accounted for by the model, the definition oft is, however, given. Fracture of the solid skeleton is accounted for using the d?/d-damage framework by Cervera et al. [32], defined in a normalized stress space. In this model, the total undamaged stress t is decomposed into its positive and negative parts In continuum damage mechanics literature, t is often referred to as the effective stress tensor. This terminology is not used here in order to avoid confusion with the effective stress concept in the sense of poromechanics. Two damage variables x AE are introduced such that Furthermore, using the effective stress concept and the definition by Gray et al. [33], the total stress tensor is where the second term on the right-hand side acts as an internal load that causes volumetric changes, i.e. shrinkage or swelling, due to changes in for example the moisture condition. The resistance of the solid matrix to withstand external and internal loads is then given by the effective stress where the elasticity tensor D e is considered to be ageindependent and given byÊ and m. The direct effect of a change in temperature enters the model through the thermal strains e th . However, it also affects the creep strains e cr , which are calculated using the Microprestress-Solidification (MPS) theory [34][35][36][37][38][39]. The creep strains thus not only depend on the current state of deformation but also on the temperature, moisture state and the extent of cement hydration, see Table 4. Apart from the inelastic strains contributions, the MPS model introduces the microprestress viscosity g f as an auxiliary state variable controlling the long-term viscous creep. Following the modified MPS theory presented by Gasch et al. [11], its evolution is given by In contrast to the original MPS theory, it is here assumed that the effective stress causes creep deformation, following Gawin et al. [40]. Hence, the used variant of the MPS theory describes not only the creep deformation but also a significant part of what usually is considered stress-free shrinkage. The latter is here considered in part as a creep deformation caused by P s . Lastly, the total strain tensor e relates the momentum balance to the dependent variable d and is defined using the standard small strain definition

Initial and boundary conditions of the mechanical problem
As for the HTC problem, initial and boundary conditions are needed to complete the mechanical part of the model. These are given as where apart from the displacements, an initial value is also needed for the microprestress viscosity which is identified as g f0 ¼ t 0 =q 4 by Jirásek and Havlásek [39].
For the conservation of linear momentum in Eq. (15), boundary conditions are additionally required, either of Dirichlet type or of a mixed type prescribing the surface tractioñ directly throught or by using a spring-type boundary condition with a stiffness f d .

Numerical implementation
The coupled system of governing equations comprised of Eqs. (1)-(4), (6)-(8), (15) and (21) is solved using the Finite Element Method (FEM) as implemented in COMSOL Multiphysics [41]. This implementation is done using standard procedures of the FEM (see for example [24,42]), starting with the conversion of the point-wise governing equations to their weak-form and application of the Galerkin method. The respective dependent variable of each governing equation is expressed in terms of its nodal quantity, such that by using the FE shape functions contained in matrices N i . Throughout this study N p , N T and N d are defined using ordinary linear Lagrangian shape functions, while the shape functions for remaining variables are chosen as a discontinuous Lagrangian function, also of first order, since their corresponding equations contain no spatial derivatives. Following the decoupling of the hygro-thermo-chemical equations from the mechanical problem, one arrives at two set of equations on the form with u 1 ¼n;ĵ; s ;p c ;p g ;T h i | and u 2 ¼d Â Ã . Notice that the auxiliary dependent variableĝ f of the mechanical problem only depends on quantities of the HTC problem and is thus not assembled into Eq. (28). To obtain the coefficients of the system matrices C 1 , K 1 , f 1 , C 2 , K 2 and f 2 , all unknowns of the weak-form governing equations are expressed in terms of their respective dependent variable using the physical and constitutive relationships described in Sects. 2.1.4 and 2.2.2 as well as in Tables 3 and 4, followed by the introduction of the FE discretization in Eq. (26).
The decoupling of the HTC and the mechanical problems allows for Eqs. (27) and (28) to be solved in series. Time integration of each equation follows a fully implicit backward differentiation formula (BDF), using an adaptive order. The solution of the partially coupled non-linear system of equations at each time increment schematically proceeds as follows: 1. Set all dependent variables in Eq. (26) to their initial value or their previously converged value. 2. Iteratively solve Eq. (27) using a full Newton algorithm to update u 1 .
-If convergence, continue to step 3 -Else, reduce time step and return to 1 3. Update the auxiliary variable g f . 4. Iteratively solve Eq. (28) using a full Newton algorithm to update u 2 . As part of the iterative solution, the inelastic strain components are updated using an exponential algorithm for the creep strains e cr and a mid-point integration rule for the thermal strains e th [11].
-If convergence, continue to step 5 -Else, reduce time step and return to 1 5. Update internal state variables for inelastic strain components (e cr and e th ) and damage (r AE ). 6. Continue to next time increment.
Compared to the fully-coupled solution employed by Gasch et al. [11], the current solution scheme is much more computationally efficient and thus better suited for large-scale structural applications.

End-restrained beam: CEOS.fr
To validate the proposed HTC-M model intended for structural scale applications, the well-documented tests of end-restrained beams performed as s part of the French research project CEOS.fr are used. These tests represent massive reinforced beams that are subjected to restrained deformation due to both selfheating and self-dessication, as well as being exposed to the outdoor environment. Measurements are available both from laboratory experiments on the concrete that can be used to calibrate the material model and for the actual test where environmental conditions, as well as mechanical and thermal data from the beam, were recorded. Buffo-Lacarriére et al. [21] summarizes the entire experimental set-up, procedure and results from the research project. Furthermore, the experiment has been investigated theoretically to validate different models for early-age concrete, both as part of the CEOS.fr project and presented in the literature, see for example [43][44][45][46]

Model calibration
The model and material parameters of the proposed HCT-M model for early-age concrete are calibrated using laboratory data from the CEOS.fr project. The concrete used for the end-restrained beams is of grade C50/60 and cast from a CEM I 52.5 N cement. The concrete mixture contains 400 kg/m 3 cement, 1765 kg/ m 3 aggregates of two different fractions and 185 kg/ m 3 water. Additionally, a small amount of superplasticizer (5.4 kg/m 3 ) was used. This yields a water-tocement ratio w=c of 0.46 and a aggregate-to-cement ratio a=c of 4.4. Data on the CEOS.fr concrete from two different laboratories is utilized in the calibration procedure, the first presented in the report by Julliot and Mary-Dippe [47] and the second by Kolani [48].

Cement hydration and ageing
The characteristics of the cement hydration and heat development were tested under adiabatic conditions by Julliot and Mary-Dippe [47]. The results of this test are presented in Fig. 1a, together with the calibrated As seen, the agreement is excellent, both with regards to the kinetics of the reaction controlled by the affinity function proposed by Sciumé et al. [14] and the final temperature. The latter is given both by the latent heat of hydration, which is here set to 415 kJ/kg cement (see Table 2), and the thermal properties of the different constituents of the concrete as summarized in Table 1.
Kolani [48] tested the evolution of several mechanical properties with the hydration including Young's modulus E, compressive strength f À and tensile strength f þ , and tensile fracture energy G þ f under isothermal conditions. The calibrated model response from the hydration model together with the ageing model in Eq. (2) for f À , f þ and G þ f are presented in Fig. 1b and compared to the measured data. Additionally, the material values predicted by the model under adiabatic conditions are presented as dashed lines to exemplify the consequence of considering the curing temperature when evaluating the material properties. As can be seen, the higher curing temperature leads to a reduction of 8 to 15 % at complete hydration compared to isothermal conditions. However, it should be pointed out that the chemical reaction is significantly faster under adiabatic conditions. Hence, for a given time the adiabatic sample can still exhibit higher strengths than the isothermal sample. The actual value obtained in-situ often lies in between these two extremes. However, for massive structures such as the CEOS.fr beam, the core can exhibit near adiabatic conditions while the concrete at the surface is significantly cooler.

Delayed deformations
Several measurements of delayed strains such as autogenous and drying shrinkage as well as creep on sealed and unsealed specimens were performed by both laboratories used to calibrate the creep and shrinkage parts of the model. The measurements reported by Julliot and Mary-Dippe [47] were carried out on prismatic specimens with dimensions 70 Â 70 Â 280 mm 3 , while those reported by Kolani [48] were carried out on cylindrical specimens with a height equal to 225 mm and a diameter equal to 120 mm. Hence the FE model of the specimens is set up using either plane stress conditions or axial symmetry depending on which specimen is modelled. In both cases, all specimens were kept at 20 C. Recalling that shrinkage in the proposed model is in part accounted for as a creep deformation caused by internal loads to the solid skeleton, it must be pointed out that the different measurements of delayed strain cannot be looked at in isolation. Consequently, the parameters related to the creep and shrinkage in the model are chosen by looking simultaneously at all reported measurements. Also, all simulations are performed with the calibrated chemical model. It is therefore mainly the parameters of the moisture transport model and those controlling creep and shrinkage that is determined using this data.
Specimens under sealed and unloaded conditions were reported by both laboratories to characterize the autogenous shrinkage of the concrete. The measurements reported by Julliot and Mary-Dippe [47] started one day after casting which is thus set as the zero level in Fig. 2a. In the figure, it can be observed that the model slightly overestimates the measured  [48] started 2 days after casting, which is why the presented autogenous shrinkage in Fig. 2b is smaller than in Fig. 2a. Still, the model prediction of the autogenous shrinkage is good also for this set of measurements, although slightly underestimated. Measuring the autogenous shrinkage at early-ages is difficult, so with this, in mind, the model performance is good especially since data from two laboratories can be reproduced. Both laboratories also reported measurements on specimens subjected to unsealed and unloaded conditions to characterize the drying shrinkage when exposed to an ambient relative humidity of 0.5. Julliot and Mary-Dippe [47] also reported the moisture loss of the specimen, which is valuable for the calibration of the moisture transport model parameters. Simulations of the moisture loss compared with the measured data are shown in Fig. 2a, and the two are in excellent agreement. After 100 days of drying the model predicts that the specimen weight has decreased by 0.8 kg. As for the corresponding deformations, the model somewhat overestimates the total shrinkage although the overall behaviour is captured. The measurements, however, only show the early shrinkage behaviour of the concrete so no conclusions can be made on the final drying shrinkage, but the model predicts a value close to 1100 microstrain for the given conditions. As for the measurements of drying shrinkage by Kolani [48], the prediction from the calibrated model is in excellent agreement with the reported measurements, but again only the early stages of drying shrinkage are reported and compared.
Creep deformations were reported by Kolani [48], both for sealed and unsealed specimens loaded in compression. The sealed specimens were loaded at either 2 or 7 days from casting, while for the unsealed specimens measurements were only reported for loading and start of drying at 2 days from casting. The applied pressure was set equal to 40% of the compressive strength at the age of loading, corresponding to 12.2 MPa at 2 days and 19.6 MPa at 7 days from casting. The results of the simulation for the three specimens are presented and compared to the measured data in Fig. 2b. The simulations of both sealed specimens are in good agreement with the measurements. However, the measured deformation on the unsealed specimen is severely underestimated, although the drying creep effect is present also in the simulation. This underestimation could likely be improved by changing the model parameters of the Micropresstress model, most notably l c . However, since the drying creep effect is also active for the sealed specimens due to the self-desiccation, this would also increase the simulated deformations of these specimens. Another explanation is that the deviation is related to the excessive delay and reverse size effect of drying creep reported by Havlásek [49]. A remedy to this was recently proposed by Rahimi-Aghdam et al. [50], but it is not considered in the current study.
A summary of all physical parameters are given in Table 1 and of all calibrated model parameters in Table 2 for the calibrated model. In total 17 of the 31 parameters in Table 2 were determined through the calibration. For the remaining parameters, typical values are assumed either due to a lack of available data or a low expected variability of the parameter.  Fig. 2 Comparison of measured (markers) and simulated (lines) results for delayed strains using specimens tested by a Julliot and Mary-Dippe [47] and b Kolani [48]. In a, also the weight loss of the specimen due to drying is reported

RG8 beam model set up
The calibrated model can now be applied to the actual end-restrained beam from the CEOS.fr project. A total of three similar beams were tested with for example different reinforcement contents and cast at different times. Here the beam entitled RG8 is analysed as it was considered a reference configuration during the research project. The geometry of the beam is depicted in Fig. 3, which includes the concrete beam and the two restraining steel struts. The central part of the concrete beam is 5.1 m long, 0.5 m wide and 0.8 m high with 2% longitudinal reinforcement content and stirrups of 16 mm bars with a spacing of 200 mm; more details are given by Buffo-Lacarriére et al. [21]. The casting of the RG8 beam started at 10:30 h on the 7th of April in 2010, and the test ended 60 days later with a static bending test. However, the bending test is not considered here as the study focuses on the early-age behaviour of concrete and the effects of restrained deformations. The beam was cast outside and thus exposed to the ambient environment; hence, quantities such as external temperature, relative humidity, wind speed and solar radiation were monitored during the entire test program. This data means the ambient conditions for the hygro-thermal part of the problem are welldefined. For all surfaces, no external loss of moisture is assumed to occur for the first 2 days after casing when an insulated formwork encloses the entire beam. After its removal, the concrete is exposed to the ambient air with the measured relative humidity and temperature used to calculate q Wg amb in Eq. (14), and a constant moisture surface coefficient f W ¼ 0:25 mm/s is assumed. The ambient gas pressure is assumed  (14) is used for all boundaries of the heat balance equation by both considering convective and radiative heat transfer. The effect of the insulated formwork is modelled by modifying the heat transfer coefficient f T accordingly. According to Buffo-Lacarriére et al. [21], two different zones regarding thermal properties of the formwork can be identified. Here it is assumed that for all lateral faces the formwork corresponds to surface coefficient equal to 1 W/m 2 K and for the upper and lower faces equal to 4 W/m 2 K, which approximately corresponds to average values calculated from the formwork drawings. The final value of f T used during these 2 days is then calculated by combining the values of the formwork with the convective heat transfer calculated for natural and forced convection due to wind speed using in a series model for the surface resistance. After the removal of the formwork, f T is only given by the convective heat transfer. For the solar radiation included in G rad , the solar position is also required and is calculated for the site of the test (approximately 48 33 0 N, 2 17 0 E) at each time instance [41]. Additionally, according to Buffo-Lacarriére et al.
[21] a 6 m high building is located 7 m from the western face of the beam causing shadowing effects during the afternoon. This effect is also accounted for when calculating G rad , by including an isothermal domain with the above dimensions and location. The isothermal domain only affects the calculation of G rad and does not add any extra degrees-of-freedom to the model. The boundary conditions of the mechanical problem are less well defined, but can have a significant affect on the results. Following drawings provided by the CEOS.fr project, the concrete beam is assumed to be supported at four locations; two supporting the central part and one support for each head. All four supports are assumed to be distributed over two finite elements and modelled as springs resisting vertical deformation. Also, a small but non-zero spring stiffness is assumed for the in-plane directions of the supports of both heads in order to cancel any rigid body motion. The reinforcement bars are modelled using a linear elastic truss formulation since no plastic deformations are expected. All rebars are connected to the concrete using a spring-type interface, where the stiffness in the normal direction degrades with the tensile damage x þ of the concrete in order to account for the loss in bond between the reinforcement and the concrete during cracking in a simplified manner. The two restraining steel struts are also modelled using a linear elastic truss formulation, where a constraint is added to suppress any out-of-plane deformation. Their connection to the concrete heads is assumed to be without any rotational stiffness.
The average element size is set to 0.1 m resulting in a mesh that in total consists of approximately 10,000 hexahedrons for the concrete, and 1600 edge elements for the reinforcement and the restraining struts, see Fig. 3. The time-stepping is set up with a maximum allowed step size of 0.1 days, although the actual step size determined by the adaptive algorithm often becomes smaller in order to resolve the continuous variation of the ambient conditions and non-linear effects of the model, e.g., damage growth with strain localization.

Results
During the test program, an extensive amount of sensors were installed in the concrete beam including temperature and strain sensors at various location. Fibre optical cables were installed in the middle section of the beam to measure its deformation. Also, the force in the reinforcement and the restraining struts are monitored throughout the test program. Buffo-Lacarriére et al. [21] gives a more detailed description of the installed sensors.

Temperature results
One of the most critical quantities measured during this type of test is the temperature evolution at different locations in the beam; especially during the early-age where the self-heating is very pronounced, but also at later times where the variation in temperature is mainly due to the ambient conditions. The measured temperature at four different locations in the mid-section of the beam is shown as dashed lines in Fig. 4 for the first 10 days of the test program. The corresponding results from the simulation are shown as solid lines, which are overall in good agreement with the measurements. The peak temperature during the hardening phase is almost in perfect agreement at the three central locations, where it is clear that the bottom part of the cross-section reaches a slightly higher peak temperature (compare blue and yellow lines in Fig. 4). One reason for this might be that the top surface is more exposed to the wind causing a higher surface heat transfer coefficient; such an assumption is used in the simulations. For the temperature sensor closest to the surface, the simulations slightly underestimate the measurement data, where it seems that too much heat is lost to the environment.
After the formwork removal, the internal temperature of the concrete beam is mainly controlled by the ambient conditions, both the temperature (black line in Fig. 4) and also the solar radiation (black area curve in Fig. 4). The latter causes the internal temperature of the concrete beam to be several degrees above the measured ambient temperature in some parts of the cross-section, as seen from both the measurements and simulation results in Fig. 4. While the measurements and simulation results are in good agreement for the location closest to the surface, the simulation seems to underestimate the temperature for the inner parts of the cross-section. The cause of this underestimation could, for example, be related to poorly assumed thermal conductivities or other material properties. Another possibility could be the mutual radiation between the concrete and the ground or other effects the ground may have on the boundary heat flux, which is not considered in the simulation. However, the overall trend observed from the measurements is recreated by the simulation.
The temperature variation can be further studied in Fig. 5 showing the temperature distribution at the midsection of the beam during the third day after casting.
The figure, also, shows the in-plane deformation of the cross-section with respect to its undeformed configuration at casting. In the first part of the figure (t ¼ 2:5 days), one can see the remaining elevated temperatures from the self-heating, where the central parts of the cross-section have a significantly higher temperature compared to its outer parts. As time moves on, the temperature decreases until t % 2:8 days when the east face of the beam is exposed to the sun, as shown by the thick border of the cross-section. This exposure significantly raises the temperature in the crosssection close to surfaces radiated by the sun, and also causes the beam to deform towards the east. When the sun is eventually shadowed by the building at t % 3:2 days at which only a small portion of the top surface is exposed to the sun, the internal temperature of the concrete starts to drop, and the cross-section moves closer to its undeformed state. As seen for example in t ¼ 3:2 days, the solar irradiation causes the thermal field of the beam to become highly irregular. Lastly, during the night the temperature continues to decrease so that the outer parts of the cross-section again becomes colder than its inner parts.

Mechanical results
Several sensors were installed during the test program to monitor the deformation of the beam. In Fig. 6 measurements of the strain-like relative displacement obtained from three cast in 2.5 m fibre optical cables placed along the mid portion of the beam at different vertical positions are shown as dashed lines. Results from the simulations evaluated as the difference in displacement over the length of the cables are shown as solid lines. As can be noticed in the figure two sets of solid lines are included. The set with thick lines corresponds to a situation where the tensile strength properties of the concrete (i.e. f þ 1 ) are homogenous. Such an assumption causes the localization of cracks to be less well-defined and cracks may not appear at the same location as observed during the test program, making it difficult to directly compare it to the relative displacement over the optical fibres. For example, during the simulation two localized cracks form over this distance (see Fig. 8), while during the test program only one single crack was observed. Moreover, these two cracks appear 7 days after casting in the simulation as seen by the jump in displacements in Fig. 6, whereas the single crack observed in the test program appears already after 3 days. For reference, it can be mentioned that two or more crack planes was observed for the other two beams of the CEOS.fr project [43]. However, these have different reinforcement content and were cast under different ambient conditions, so the results are not directly comparable.
To investigate the observed deviation, a second simulation indicated by the thin solid lines in Fig. 6 was performed. In this simulation, a weak-zone is introduced in the mid-section as suggested by Buffo-Lacarriére et al. [43] motivated by the high concentration of instrumentation at this location. The weakzone is introduced by reducing f þ 1 with 5 % over a length corresponding to a single mesh element. This modification forces a single crack to appear in the weakened section, and it does so at approximately the same time as observed during the test program, i.e. at t % 3 days after casting. After removal of the formwork, the simulations seem to underestimate the contraction of the beam slightly as seen for the displacement magnitude when the first crack appears (t % 3 days). This might point to an issue with for example the assumed mechanical boundary conditions or the connection between the concrete heads and the two steel struts which have been simplified, causing the restraint to be overestimated. It could also be related to the material response underestimating, for example, the thermal contraction during the cooling phase after the removal of the formwork. However, after cracking it should be emphasized that the overall trend and magnitude of the deformation is recreated by the simulations, where the top of the beam is exposed to a larger contraction than the bottom due to the increased curvature caused by the cracking. In order to fully verify the mechanical response of the beam, a more thorough analysis of the mechanical boundary conditions is definitely required. Moreover, it might be necessary to include the random spatial distribution of especially mechanical strength properties and possibly other properties as well.
Lastly, the simulated crack pattern is shown in Figs. 7, 8 and 9 for the east and west faces of the beam where e 1 is the maximum principal component of the total strain tensor e. The crack width is estimated in the order of 0.25 mm after 60.5 days, which can be compared to an estimated crack width in the order of 0.1 mm obtained from the fibre optical cables at 3.5 days [21]. Looking at Fig. 7 that shows the obtained crack pattern at 6.5 days after casting, it can be observed that larger cracks have mainly formed at the two ends of the beam. This fact is likely a consequence of the geometric discontinuities at these locations.
Continuing to the end of the analysis in Fig. 8 (t ¼ 60:5 days) the number of cracks have increased significantly. The simulated crack pattern is now dominated by two new cracks in the middle part of the beam in addition to those at the ends that were also present in Fig. 7. Moreover, several surface cracks are now visible, especially at the east face, which is likely caused by a combination of the continuous drying of the cross-section and the thermal dilation. Interesting to notice is the effect of the solar irradiation on the cracking, where the east face shows both larger cracks and the extensive surface cracking compared to the west face that is shadowed by the building seven meters from the beam during the afternoon. The same general conclusions can be made from Fig. 9, showing the crack pattern from the simulation where a weakzone was introduced.

Comment on the numerical solution
A significant goal of the present study is to show how the model formulation and solution strategy of a fullycoupled HTCM for early-age concrete can be adapted for structural scale applications. For comparison, a simulation of the RG8 beam was also done with the fully-coupled model presented by Gasch et al. [11]. The overall computational time increased by a factor of five when using that model, while only insignificant differences were observed for structural scale quantities as those presented in previous sections.

Conclusions
This paper presents a multiphase model for the hygrothermo-chemo-mechanical behaviour of early-age concrete intended as a predictive analysis tool for use in structural applications. The formulation of the model is a modified version of the fully-coupled HTCM model by Gasch et al. [11] in which the concrete is considered as a three-phase material consisting of a single solid phase and two fluid phases. Through a number of simplifications of the model by Gasch et al, mainly the decoupling of the hygrothermo-chemical problem from the mechanical problem, it is shown how it is possible to adopt a more efficient solution strategy to reduce the computational cost of the model. For the studied structural scale application, this reduced computational cost is obtained without any notable loss in accuracy. This model property can be vital for three-dimensional cases, where the number of degrees of freedom can be significant. Using this modified model, essential effects for the early-age behaviour concrete, e.g. cement hydration, self-desiccation, self-heating and strength growth, as well as mechanical effects also vital for mature concrete such as cracking, drying shrinkage and creep can be considered while maintaining a reasonable computational effort. With this type of model, questions central to improve the understanding of early-age cracking of concrete structures can be investigated. Such questions among others include the importance of creep for relaxing the restrained deformations, and whether one should place most care on the thermal or moisture treatment of a structure.
To show the applicability of the model, one of the end-restrained beams from the French research project CEOS.fr was simulated also considering the effects of realistic environmental conditions. The latter was done using measured data of the ambient relative humidity, temperature, wind speed and solar radiation. This data was used to define the boundary fluxes of the hygro-thermal problem accurately. Furthermore, the behaviour of the concrete was calibrated using laboratory data on the used concrete mixture. The simulations showed how the proposed model was able to recreate the behaviour of the end-restrained beam while there was still some discrepancies when qualitatively comparing the measurements and simulations. These discrepancies were mainly observed for the mechanical response, where the boundary conditions are less well-defined and the behaviour more non-linear due to cracking. While not a truly predictive analysis, the simulations showed how the model is useful in modelling and understanding structural damages due to restrained deformations.
The constitutive equations that are required to complete the hygro-thermo-chemo-mechanical model are summarized in Tables 3 and 4. All equations are derived or taken directly from Gasch et al. [11] if no other references are given in the tables. Table 3 Constitutive equations for the hygro-thermo-chemical part of the model Affinity function [14] A C ¼ Capillary degree of saturation s w ¼ exp g s À 10n C À Á u Â Ã À 1 exp g s À 10n C À Á À 1 , and s g ¼ 1 À s w