A practical creep model for concrete elements under eccentric compression

Many prestressed concrete bridges are reported to suffer from excessive vertical deflections and cracking during their service life. Creep softens the structure significantly, and therefore an accurate prediction of creep is necessary to determine long-term deflections in elements under eccentric axial compression such as prestressed concrete girders. This study proposes a modification to the creep damage model of Model Code 2010 to account for the effect of load eccentricity. The modified creep model considers damage due to differential drying shrinkage. Initially, the creep behaviour of small scale concrete specimens under eccentric compression load is investigated experimentally. Twelve small-scale concrete prisms were subjected to eccentric axial loading to assess their shrinkage and creep behaviour. The main parameters investigated include the load eccentricity and exposure conditions. Based on the experimental results, an inverse analysis is conducted to determine the main parameters of the modified creep model. Subsequently, a numerical hygro-mechanical simulation is carried out to examine the effect of load eccentricity on the development of shrinkage and creep, and on the interaction between drying, damage and creep. The results indicate that eccentric loading leads to different tensile and compressive creep through the cross section, which contradicts the current design approach that assumes that tensile and compressive creep are identical. The proposed model also predicts accurately the long-term behaviour of tests on reinforced concrete elements available in the literature. This study contributes towards further understanding of the long-term behaviour of concrete structures, and towards the development of advanced creep models for the design/assessment of concrete structures.

were subjected to eccentric axial loading to assess their shrinkage and creep behaviour. The main parameters investigated include the load eccentricity and exposure conditions. Based on the experimental results, an inverse analysis is conducted to determine the main parameters of the modified creep model. Subsequently, a numerical hygro-mechanical simulation is carried out to examine the effect of load eccentricity on the development of shrinkage and creep, and on the interaction between drying, damage and creep. The results indicate that eccentric loading leads to different tensile and compressive creep through the cross section, which contradicts the current design approach that assumes that tensile and compressive creep are identical. The proposed model also predicts accurately the long-term behaviour of tests on reinforced concrete elements available in the literature. This study contributes towards further understanding of the long-term behaviour of concrete structures, and towards the development of advanced creep models for the design/assessment of concrete structures.

Introduction
Due to their structural efficiency and cost-effectiveness, prestressed concrete (PC) box girders are extensively used in the construction of long-span bridges (100-300 m). The cross section of box girders is typically subjected to bending and eccentric compression forces due to the opposite action of dead loads and prestressing forces. Eccentric compression forces can generate second order effects that affect the long-term behaviour of the structure. For instance, in recent years, many concrete box girder bridges have been reported to suffer from excessive mid-span deflections and cracking, which affect their safety and serviceability [1][2][3]. Whilst time-dependent vertical deflections of bridges can be calculated using existing design guidelines, actual measured deflections usually exceed the calculations due to the inaccuracy of existing creep and shrinkage models [4][5][6].
Due to the technical and economic constraints in testing and analysing large PC box girders, creep behaviour is generally examined using simple beam or column elements. For instance, previous research has investigated the influence of second order effects [7,8] and material nonlinearity [9][10][11] on the creep behaviour of such elements. In these studies, the magnitude of tensile and compressive creep strains was assumed to be equal. Existing state-of-the-art design guidelines (e.g. Model Code 2010 (MC2010) [12]) also assume that the compressive and tensile creep are identical. However, recent experimental evidence indicates that the tensile creep strains can be between 2.4 and 8.9 larger than the compressive creep strains [13,14]. The difference between tensile and compressive creep has been attributed to microcracking and drying-induced microprestressing [15]. The force applied on a concrete element affects significantly the development of drying-induced microcracking. For instance, compressive forces tend to close microcracks and inhibit further microcracking, whereas tensile forces promote the growth and widening of microcracks during the concrete drying process. The microprestress-solidification model [16] attributes the drying creep to microcracking and microprestressing effects. This model has been widely used in recent studies on creep [17][18][19][20][21] and expresses the rate of the stress-induced shrinkage as a function of the rate of internal humidity loss. An alternative creepdamage model was proposed based on concrete damage [22,23], where nonlinear creep is accounted for using the interaction between the growth of creep strain and damage over time. However, such models neglect the influence of shrinkage and drying-induced microcracking, which always occur during the service life of structures.
To account directly for the effect of microcracking due to the non-uniform distribution of internal humidity, recent studies [24][25][26][27][28] have adopted a hygro-mechanical approach to assess creep in concrete elements. This approach calculates the nonuniform development of stresses and the dryinginduced cracks using a shrinkage coefficient that correlates the internal humidity and shrinkage deformations. Drying-induced creep is accounted for through a Kelvin-Voigt chain associated with an ageing dashpot [29][30][31]. Although this approach has proven effective at estimating creep in concrete elements, tests are necessary to obtain the input parameters. Moreover, the analysis is computationally demanding and thus difficult to implement in practice. On the other hand, the creep and shrinkage models included in current design codes (e.g. MC2010 [12] and ACI 209 [32]) do not take into account the interaction between drying, creep and damage explicitly, which can lead to inaccurate and unsafe calculations of long-term deformations of structures with drying-induced cracks. Consequently, a more robust and practical model is needed to assess the timedependent behaviour of RC structures that considers the combined effect of (1) eccentric compression, and (2) ageing, drying and cracking effects.
This study proposes a modification to the creep damage model of Model Code 2010 to account for the effect of load eccentricity in concrete. The proposed model considers damage occurring during the creep process. Section 2 of this article presents an experimental programme on eccentrically loaded smallscale specimens (prisms) that assesses their tensile and compressive creep, as well as the interaction between drying-induced damage and creep strains. Subsequently, Sect. 3 proposes a new modelling approach that simulates the interaction between drying, creep and damage using a one-way hydro-mechanical coupling model. In Sect. 4, the model proposed in Sect. 3 is used to determine the main parameters of the new creep model using the test results from the prisms. The proposed model is subsequently validated using results from long-term behaviour tests on reinforced concrete beams available in the literature. Concluding remarks of this study are given in Sect. 5. The results of this study contribute towards developing a further understanding of the long-term behaviour of concrete structures, as well as towards the development of practical models to evaluate the serviceability of RC structures.
2 Experimental programme

Test specimens and material properties
To examine the effect of load eccentricity and drying of concrete over time, twelve 100 9 100 9 450 mm concrete specimens (prisms) were tested. Specimens D1, D2, D3 and D4 were tested in concentric compression to obtain creep strains, whereas specimens W1, W2, W3 and W4 were tested in eccentric compression. Likewise, two specimens (DS1 and DS2) were used to measure drying shrinkage strains, whereas two specimens (BS1 and BS2) were epoxy coated on all faces to examine basic shrinkage without concrete drying effects.
The specimens were cast using one batch of normal-strength concrete typical of bridge construction in China. A w/c ratio of 0.40 was used according to the following mix proportions (kg/m 3 ): water 163, cement 410, sand 545, and gravel 1275. The 28-day mean concrete compressive strength f c was 54.5 MPa, as obtained from concrete cubes tested according to the Chinese bridge design code JTG D62 [33]. All the specimens were moist-cured for 3 days after casting. After this, epoxy coating was applied onto the faces of specimens BS1 and BS2 to prevent drying, whereas DS1 and DS2 were fully exposed to the laboratory environment.

Test setup and instrumentation
The specimens were tested in compression using the standardised test rig in ASTM C512 [34], but with some minor modifications. The self-reacting test rig consisted of three horizontal steel plates (1-3, as shown in Fig. 1) and three vertical steel bars. The concrete specimens were placed between plates 1 and 2. The compressive load was applied using a hydraulic actuator of 500 kN capacity, which in turn reacted against plate 3 during the tests.
To investigate the linear and nonlinear creep under direct compression, the following axial loads were applied: 100 kN on specimens D1 and D2, and 250 kN on specimens D3 and D4. These loads (equivalent to normalised axial load ratios of 0.18 and 0.45, respectively) simulated typical moderate and heavy loading conditions. Eccentric compression creep tests were carried out on specimens W4, W3, W2 and W1 under a sustained load of 100 kN, and using eccentricities of 5, 10, 15 or 20 mm, respectively. The eccentric load was applied using steel plates of width 60, 70, 80 or 90 mm (see Fig. 1) fixed to the specimen ends. These plates were rigidly fixed using four 8 mm bars (welded to 100 9 100 mm steel plates, see Fig. 1) cast inside each specimen for a length of 50 mm. To ensure a uniform application of the load, the 100 9 100 mm steel plates were used as part of the mould during casting and vibration of the specimens. Accordingly, no gaps existed between the plates and the specimens, as shown in Fig. 1. However, in a few cases where small voids (\ 5 mm in size) existed between the 100 9 100 mm steel plates and the specimens ends, epoxy padding material was added to fill in such voids. Spherical steel bearings at the top and bottom of the specimens (plates 1 and 2) allowed free end rotations of the samples. This ensured the application of vertical load with no bending moments during the tests. The loading on specimens D1-D4 and W1-W4 started on day 9 after casting (casting day = day 0). The applied vertical force was monitored and controlled by a pressure sensor in the hydraulic jack. All the specimens were kept in standard laboratory conditions at ambient temperatures of 20 ± 2°C and RH = 70% ± 10% during testing. It should be noted that the authors tried to test a fifth specimen with an eccentricity of 25 mm, but the specimen failed suddenly after a few hours of testing due to excessive bending.
The instantaneous strains due to the applied force on specimens D1-D4 and W1-W4 were recorded during the tests. Before applying the load, instantaneous deformations were monitored using four vibrating wire strain gauges, one fixed on every face of each specimen. Two temporary strain gauges installed on the direction perpendicular to the specimen's flexion plane ensured that the loading was applied uniformly in this direction. During the initial loading process, all the measured strains were compared with the numerical results to adjust the applied bending moments and the position of the specimens in the perpendicular direction. The development of strain over time was measured with two vibrating wire strain gauges installed on two opposite faces of each specimen: P1 on the compression face, and P2 on the tension face. The gauges had a strain range of ± 1500 le and accuracy of 1 le, which were sufficient for the expected strains. Measurements started on day 4 after casting, whereas creep data from D1-D4 and W1-W4 were recorded from day 10 onwards. The shrinkage and creep data were recorded for 321 days and 331 days, respectively.  Figure 2b shows the development of the total strain (i.e. the sum of instantaneous, shrinkage and creep strains) over time. The results in Fig. 2b show that the mean values of the instantaneous strains are -305 le and -705 le for axial loads P of 100 and 250 kN, respectively. Figure 2b also shows that, at the end of the tests, the shrinkage and creep strains for P = 250 kN is about 2.4 times that for P = 100 kN. Figure 3 compares the instantaneous strains recorded during the tests on specimens W1 to W4 and the results given by MC2010. The instantaneous strains on the compression face of the specimens ranged from -610 le (W1) to -305 le (W4), while the instantaneous strains ranged from 60 le (W1) to -200 le (W4) on the tension face. The figure also shows that, as the load eccentricity increases, the differential strain between the compression and tension faces increases from 516 to 2172 le at day 331, thus confirming that the bending moment induces additional curvature deformations. A significant nonlinear relationship between curvature deformation and eccentricity is also observed. For instance, as the eccentricity increases from 5 mm (specimen W4) to 20 mm (specimen W1), the differential strain (i.e. the curvature deformation of the specimen) increases 4.21 times, which is not proportional to the increase in applied bending moment (four times). This may be attributed to different creep behaviour under compression and tension, and to additional eccentricity induced by lateral deflections. The results in Fig. 3 also show that the MC2010 results match well the test results on the tension face but not so well on the compression face. Hence, there is a need to address this issue by coupling damage and creep on the longterm deformation potential of concrete.

Proposed modified creep model
The proposed model considers the combined effect of ageing of concrete, coupling effects of damage and nonlinear creep strain. These effects are accounted for using a damage model (described in Sect. 3.1), a drying and shrinkage model (Sect. 3.2), and a creep model (Sect. 3.3), which are integrated and solved numerically in a commercial software (Sect. 3.4). A full list of symbols used in this section is included in Table A of the Supplementary Material.

Damage model
For computational efficiency, this study adopts the isotropic elastic damage model proposed by Mazars [35] to analyse the creep and drying-induced stresses [22,23,31]. The material properties used in Mazars' model depend on the concrete age, especially at the early stages when concrete is most vulnerable to damage. In this article, some parameters of the original Mazars' model are modified to consider the changes in the strength and stress-strain relationship over time.
The time-dependant constitutive law is defined as: where r and e e are the elastic stress and strain, respectively; D(t) is the damage variable ð0 Dt ðÞ 1Þ at a concrete age t; and E 0 is the initial stiffness matrix of the material. To consider the behaviour under tension and compression stresses, the damage variable D(t)is defined as the linear superposition of a tensile damage variable D þ t ðÞ , and a compressive damage variable D À t ðÞ :  Total strain(µε) where A þ ,B þ ,A À and B À are empirical parameters obtained from uniaxial material testing; a is a ratio between tension and compression stresses; k 0 t ðÞis the damage threshold at a concrete age t; and e eq is an equivalent strain. The value a is calculated as: where hr i i is the ith principal effective stress which is controlled by the Macaulay function. The Macaulay brackets in Eq. (3) indicate that only positive stress is considered.
For uniaxial tension, k 0 t ðÞcan be calculated as: where f ctm t ðÞis the mean tensile strength of concrete; and E 0 t ðÞis the initial modulus of elasticity. The damage criterion function, which is linked to D(t), can be used to determine the damage level in a multi-axial stress state. The strain-based damage criterion function x is expressed as: where all the variables are as defined before.
To consider the contribution of creep strains to the development of damage over time, Mazzoti and Savoia [23] proposed the use of a weighting coefficient b c to calculate the equivalent strain e eq , which is defined as: where e i is the ith principal elastic strain; and e ic is the creep strain in the ith principal direction. The Macaulay brackets in Eq. (6) indicate that only positive strain is considered. Typical values of b c range from 0.05 to 0.40 [22,23,31].
In real situations, the nonlinear creep strain is coupled with the plastic strain, which increases microcracking over time and can lead to structural failures at high stress levels. The affinity hypothesis [36] has been used to describe the nonlinear effect of a high stress/strength ratio on the linear creep coefficient (b ; , as described later). In this study, the affinity hypothesis is adopted to represent the influence of the stress/strength ratio on b c : where a c is a modification coefficient for creepinduced damage which needs to be calibrated using test results;r i is the ith principal stress; f c is the concrete compressive strength; HF is a heaviside function defined as HF= 1 for e eq À k 0 t ðÞ!0 or HF= 0 for e eq À k 0 t ðÞ\0 (i.e. a discontinuous function to control the effect of creep strain). Therefore, b c is taken into account only after the damage criterion [i.e. Eq. (5)) has been reached.
Equations (4)to (7) show that the effect of concrete ageing on damage is controlled by the tensile strength f ctm t ðÞand the initial elastic modulus of concrete E 0 t ðÞ . In this article, the strength development model in MC2010 [12] is adopted to calculate f ctm t ðÞand E 0 t ðÞ using the test results from Sect. 2.
It should be mentioned that, in real concrete structures, the evolution of damage and the development of material properties are mutually linked. For instance, as the strength and elastic modulus of concrete increase over time, the threshold of damage and the damage criterion also change over time. However, limited experimental data exist on the interaction between damage and strength increase. Consequently, the calculations in this study assume that the strength and elastic modulus do not change once the damage criterion is reached.

Drying and shrinkage model
The time-dependent strain of concrete consists of both shrinkage and creep strains. The hygro-mechanical simulations in this study consider the effect of both drying and autogenous shrinkage. According to Fick's Second Law, the diffusion equation with the effect of self-desiccation [37] can be expressed as: where D(H) is the moisture diffusion coefficient; H is the relative pore humidity; f is the surface factor; H en is the environmental humidity; H s is the relative humidity on the exposed face S and oH a ot is the rate of internal humidity change due to self-desiccation (full details of the self-desiccation model can be found in Ref. [6]). To determine the parameters needed for the hydro-mechanical model, a moisture diffusion analysis is first performed. This study adopts the approach suggested by MC2010 [12], where D(H) is calculated as: ; and H c is the relative pore humidity at D(H) = 0.5D 1 . MC2010 suggests using b = 0.05, H c = 0.80 and n = 15 in practical calculations. The hydro-shrinkage coefficient b sh couples the free shrinkage strain and the varying moisture. Hence, the incremental shrinkage strain can be calculated as: where DH is the variation of pore humidity. The result of the shrinkage strain calculated with the above equation is substituted later in Eq. (17) to consider the interaction between drying, damage and creep (see Sects. 3.4.1 and 3.4.2).

Creep model
In this study, the Dirichlet series is used to calculate the creep coefficient ;: where g j t 0 ðÞis the jth age function; t 0 is the initial loading time (in days); t is the age of concrete (in days); and s j is the jth discrete retardation time.
The creep strain at the end of the ith and ith?1 time steps considering the effect of concrete ageing can be expressed as: where t i and t iþ1 are the loading ages; g j t i ðÞis the jth age function; and s j is the jth discrete retardation time. From Eq. (12), the creep strain increments from t i to t iþ1 are given by: where Note that during the explicit iteration process, the stress is updated at the beginning of each time step and remains unchanged within a time step. Dr t i ðÞis taken as the size of each time step.
From Eqs. (12)to (14), the incremental relationship can be defined by: To carry out the above creep incremental analysis, the creep coefficient ; [Eq. (11)] needs to be calibrated/calculated so that the parameters g j t i ðÞand s j can be determined. In this study, the creep model in MC2010 [12] is adopted for such calculation.
For simplicity, the loading time function and concrete age function are represented by curves. The age function g j t i ðÞ is assumed as constant at t 0 ¼ 1, 10, 100, 1000 and 10,000 days. At t 0 ¼ 10 days, Eq. (11) can be expressed as: where ; 10; t ðÞ is the creep coefficient at a loading age of 10 days.
By approximating the creep curve in MC2010, g 1 10 ðÞ to g m 10 ðÞ can be determined at t 0 ¼ 10. Subsequent values of g j can be calculated at any loading age t i . The jth age function g j t i ðÞ is established by fitting the values of the discrete points from g j 1 ðÞ to g j 10000 ðÞ . The results indicate that the exponential expression (for m = 5) reproduces accurately the creep model included in MC2010.

Coupled hygro-mechanical numerical modelling
The above mentioned shrinkage, creep and damage models are implemented into the subroutine CUSER2 for 2D elements of ADINA Ò [38]. The time-dependent stress integration and numerical algorithm are described in the following sections.

Humidity and shrinkage
The development of shrinkage deformation is simulated using the drying process described previously in Sect. 2.2. The consumption of internal humidity due to self-desiccation is modelled using a time-dependent uniform drying process. In each time step from t i to t iþ1 , the incremental shrinkage strain De sh is calculated as: where b sh is the hydro-shrinkage coefficient; and H t iþ1 and H t i are the humidities at t iþ1 and t i , respectively.

Creep and damage
In the elastic state (i.e. damage criterion function x\0), the modulus of elasticity only depends on the age of concrete. The incremental strain and stress can be calculated using the iterative procedure described previously [Eqs. (13)- (15)]. For x ! 0, damage starts to develop and thus the damage-induced nonlinear strain-stress relationship needs to be used for the calculation of stresses [Eq.
(1)]. The evolution of damage is controlled by Eq. (2). Accordingly, the elastic strain increases from e t i to e t i þ1 in each time step from t i to t iþ1 . This in turn increases the damage variable from D t i to D t iþ1 , thus reducing the effective elastic modulus. The elastic strain e t iþ1 can be expressed as: where De, De c and De sh are the increments of total strain, creep strain and shrinkage strain, respectively. The increment of creep strain De c is calculated using Eq. (13). The increment of shrinkage strain is obtained using the humidity variation from t i to t iþ1 . The elastic stress r t iþ1 is calculated as: The incremental stress Dr, which is used to update the creep model parameters, can be obtained using: The subroutine in ADINA Ò uses an explicit Euler forward stress integration approach and the abovementioned algorithm. In summary, for each time step from t i to t iþ1 , the calculation involves the following steps:  (15)] and save the data, then go to next the time step and repeat the process until the calculation time is reached.
The following section adopts the proposed approach to provide further insight into the long-term behaviour of the specimens tested in Sect. 2. The approach is further validated using long-term test results from concrete specimens available in the literature.

Numerical analysis and discussion of results
This section describes Finite Element Analyses (FEA) that incorporate the diffusion model (described in Sect. 3.2) to simulate the hydro-mechanical behaviour of the tested specimens, as well as the mechanical models to assess the time-dependent shrinkage and creep strains considering the interaction between drying, damage and creep. This is conducted through an ad hoc subroutine developed in ADINA Ò [38], as described in Sect. 3.4.

FE modelling
The FEA package ADINA Ò was used to solve the proposed coupled diffusion-mechanical model. The internal distribution of RH was calculated using ADINA-T, while ADINA Structure was employed for the analysis of the time-dependent structural behaviour. For consistency, the mesh of the diffusion model and mechanical model was the same. In the diffusion model, 2D diffusion elements were adopted to calculate the internal distribution of RH, and 1D convection elements were used to model the drying boundary of the specimens. Based on the RH measured in the laboratory, a constant ambient RH = 70% was selected for the exposed face in the diffusion model. 2D plane stress elements were used in the mechanical model.

Modelling parameters
In this study, the parameters needed for the new analytical creep model were determined through inverse analyses, calibrated against the experiments described in Sect. 2.

Drying and shrinkage
The surface factor f and the hydro-shrinkage coefficient b sh were calculated through inverse analyses using the results of sealed (BS) and unsealed specimens (DS). In the inverse analyses, the values f were first varied from f 1 to f n , while b sh remained constant. Subsequently, b sh was varied from b sh1 to b shm , which in turn led to different values of f i . This resulted in m Â n sets of results for different combinations of f i and b shi . The inverse analyses with the test results indicate that, as expected, there is a nonlinear relationship between the calculated shrinkage strain and the surface factor f . Through a parametric study (see Fig. S1 in Supplementary Material), the best match between the analysis and test results is obtained for f = 0.70 mm/day and f = 0.01 mm/day for unsealed and sealed specimens, respectively, and therefore such values are adopted in subsequent analysis. A value b sh = 0.00234 matched well the experimental results of sealed and unsealed specimens, and therefore this value is used in the proposed analytical creep model.

Creep and creep-induced damage
To determine the linear creep coefficient b ; and the creep-induced damage modification coefficient a c [defined in Eq. (7)], the test data from specimens D and the linear creep model in MC2010 [12] (see Sect. 3.3) were used. To consider the errors due to variability of material properties and environment, MC2010 suggests the use of a mean creep coefficient with a 95% confidence interval, which leads to lower and upper limit values of 0.59 and 1.41, respectively. In this study, the mean total strain of specimens D1 and D2 (subjected to P = 100 kN) is calculated by performing FEA assuming b ; = 1.5, 1.25 or 1.0. The results from a parametric analysis indicate that the use of b ; = 1.25 matches accurately the test results. Since the compressive stress produced by the applied load is below the linear limit value (0.4f cm ), the calculated results confirm that the specimens do not experience damage.
In this study, the creep model in MC2010 is adopted for creep incremental analysis. The function b RH ðÞ (defined in MC2010) is used to take into account the effects of relative ambient humidity. Since all the specimens were kept under standard laboratory conditions at ambient RH = 70% ± 10% during testing, the approximate constant RH = 70% is adopted for creep calculations.
Unlike specimens D1 and D2, specimens D3 and D4 were subjected to a 250 kN compressive load, which led to some compressive damage in the specimen core and nonlinear creep strain. To examine the interaction between damage and nonlinear creep, the modification coefficient for creep-induced damage a c was varied from 0.2 to 1.0 in the FEA of specimens D3 and D4. The results in Fig. 4a indicate that the best matching of total strain was obtained for a c = 0.6. It is also shown that, at day 398, the calculated total strain increased from -2595 to -3260 le,asa c increased from 0.2 to 1.0. Figure 4b shows that as a c increased, damage also increased from 0.13 at 10 days to 0.37 at 398 days. This confirms that the creep-induced damage modification coefficient a c has a significant influence on the development of compressive damage over time.
To study the relationship between tensile creep and compressive creep, tensile creep strains were calculated using the same material parameters used previously in the compressive creep FEA models. The applied load was determined using a stress/strength ratio of 50% for both tensile creep and compressive creep models. The specific creep (creep strain per unit stress) is adopted here to show more clearly the difference of creep strains under different applied stress. Basic creep was also calculated to examine the contribution of the drying process to the total creep strain over time. To avoid the influence of drying shrinkage and creep, an ambient RH = 100% was selected for basic creep analysis. Figure 5 compares the calculated specific creep in compression and tension for a stress/strength ratio = 50%. The total creep (basic creep plus drying creep) results show that the specific creep stain in tension is significantly higher than in compression, and that the ratio of tensile to compressive specific creep is 1.78. It is also shown that the calculated basic creep in tension is equal to the one in compression, and that the ratio of basic creep to total creep is 0.24 and 0.43 for tension and compression, respectively. These results indicate that the difference between tensile and compressive creep can be attributed to drying-induced damage.

Drying-induced stresses and damage
The results summarised in the previous paragraphs are used here to simulate the drying-induced stresses and damage, as well as the creep behaviour of the tested specimens. Figure 6a-b show, respectively, the FE results of time-dependent RH and vertical stress distribution across the 100 mm width (at mid-height) of the shrinkage specimen DS. Figure 6a indicates that, as expected, the calculated internal RH decreases over time from 92% at t = 10 days to 77.5% at t = 345 days. It is also shown that, as expected, the RH is not uniform across the section due to surface drying.  For instance, at t = 10 days, the maximum RH = 92% is located in the middle of the section, whereas the minimum internal RH = 77.7% is at the surface.
The results in Fig. 6b show that drying causes vertical stresses due to the differential drying shrinkage strains throughout the cross-section of the specimen. At an early age (t = 10 days), the stress is not uniform across the width of the specimen. Tensile stresses develop up to depths of 20-22 mm, whereas the core remains in compression. Figure 6b also shows that, at t = 10 days, the tensile stresses exceed the concrete tensile strength at the surface, thus leading to damage. This is confirmed by a drop in the maximum tensile stress from 2.0 to 1.1 MPa at a 5 mm depth from the specimen surface. The maximum compressive stress is -1.0 MPa in the centre of the specimens. Note that the magnitude of drying-induced stresses decreases over time as the internal RH tends to balance with the ambient RH. For instance, at t = 345 days, the maximum tensile stress drops to 0.17 MPa, which is only 15% of that calculated at t = 10 days.
For the free shrinkage specimens, shrinkage dominates the total deformation of the specimens. However, creep significantly affects the distribution of drying-induced stress during the drying process of the specimens. To examine the contribution of creep to the variation of stress and damage distribution over time, the specimens were modelled without creep. Figure 7a-b compare the FEA results of drying specimens at t = 10 days with and without creep, respectively. Figure 7a shows that, due to stress relaxation in concrete, the drying-induced damage occurs mainly near the surface up to depths of 2.5-5 mm. Conversely, Fig. 7b shows that if creep is neglected, the tensile stresses caused by non-uniform drying cannot be relieved, thus leading to damage (cracking) growing from the surface to a depth of up 30 mm at t = 10 days. As discussed in Sect. 3.3,n o cracking was observed during the tests on drying specimens DS1 and DS2. These observations imply that the inclusion of concrete creep in the proposed model is necessary to simulate and assess accurately the behaviour of drying concrete. The results also indicate that, as expected, concrete creep has no influence on the development of shrinkage strain in sealed or unsealed specimens. However, creep significantly affects the drying-induced stresses and the damage distribution.

Creep under eccentric compression
In this section, specimens W1-W4 are modelled to examine the time-dependent strains and stresses. For the diffusion and mechanical models, the following values were adopted: surface factor f = 0.7 mm/day, hydro-shrinkage coefficient b sh = 0.00234, creep modification coefficient b ; ¼ 1:25, and creep-induced damage modification coefficient a c = 0.6 (see more details in Sect. 4.2). The input parameters used in the   Table B and Table C in the Supplementary Material. Figure 8 compares the FE results and strain gauge readings from gauges P1 and P2. The total strain shown in Fig. 8 includes the transient elastic strain caused by the applied load, as well as the time-dependent shrinkage and creep strains. The results indicate that the FE models match well the experimental data and provide more accurate results than MC2010 for the compression face of the specimens. The minor discrepancy between calculations and test results can be attributed to minor variations of (i) environmental conditions (e.g. variable RH), and (ii) applied load during a test. To examine the impact of different load eccentricities, Fig. 9a-c compare the FE results of vertical strain, stresses and the damage variable D(t) at the mid-height at t = 387 days (i.e. at the end of the test), respectively. To assess the difference between the proposed model and MC2010, the simplified elastic beam theory is adopted in the calculations. Figure 9a shows that, as expected, the slope of the vertical strain calculated using the elastic beam model is constant and varies linearly for specimens W1 to W4. However, when non-uniform drying shrinkage and nonlinear creep are considered, the slope of the vertical strain is steeper than the results from the beam model, while the vertical strain on the compression face is higher. It is also shown that the ratio between the maximum compressive strain given by the proposed model and MC2010 decreases nonlinearly from W1 to W4. For example, such ratio is 1.13 for W1, but only 1.05 for W4. It is also evident that the proposed model Total strain(µε) Compression face calculates a roughly constant slope of vertical strains for W1, which indicates that there is no cracking at the cut-line of W1. These results match well the experimental observations shown later in Fig. 10. As shown in Fig. 9b, the slope of the vertical stresses (which include normal and bending stresses) is approximately constant in the core of the specimen. However, due to drying shrinkage, vertical stresses reduce close to the specimen surface, particularly within the first 5 mm of the compression face. For example, the vertical stress at the compression face of W1 drops by about 10 MPa.
The results in Fig. 9c also show that, due to drying and bending, the magnitude and distribution of damage are different on the tension and compression faces of the specimens. On the tensile face, damage extends up to a depth of 5-7.5 mm. It is also shown that the damage variable depends on the instantaneous strain produced by the load. The damage variable decreases from 0.86 (W1) to 0.22 (W4) as the eccentricity decreases. The results in Fig. 9c also show that the depth of damage on the compression face of W1 extends up to 50 mm (i.e. half of the specimen width), where the damage threshold has been reached. Overall, the results in Fig. 9a-c indicate that both concrete drying and the applied load can produce and affect the development of damage. To further study the effect of eccentricity on the development of curvature, Fig. 10 compares the curvature of specimens W1-W4 at t = 321 days, as well as curvature given by the FE models and by MC2010. The experimental results in Fig. 10 were calculated using strains measured at both faces of the specimens. The figure also shows FE results of three models with theoretical eccentricities of 25, 30 and 35 mm, which are used here to investigate the trend of results beyond the eccentricities used in the tests. The results in Fig. 10 show that whilst the curvature increases initially linearly with the eccentricity, it becomes nonlinear for eccentricities larger than 25 mm. The rapid increase in curvature can be attributed to extensive damage on the tensile face of the specimens, as indicated by the cracks shown in Fig. 10. It is shown that the proposed model and MC2010 estimates accurately the curvatures for eccentricities below 25 mm, but the results diverge for large eccentricities since the model in MC2010 neglects the combined effects of cracking and drying. For instance, the calculated curvature for an eccentricity of 25 mm is 34.5 9 10 -3 /m, which is 1.3 times larger than the MC2010 results.
The proposed model can also estimate the development of cracking damage in concrete. Figure 11 shows that after applying the load (t = 10 days) in the model with a theoretical eccentricity of 25 mm, damage occurs at the tensile face and across a depth of 5-10 mm. Tensile damage extends towards the centre of the specimen over 1 day (t = 11 days), and the number and length of the cracks increase significantly, which somehow explains why the fifth specimen mentioned in Sect. 2.2 failed prematurely during the test. Tensile damage progresses until the end of the analysis (t = 387 days) where the depth of tensile damage increases up to 25 mm.

Validation of proposed approach
The proposed approach is used to simulate the longterm behaviour of two simple-supported reinforced concrete beams tested by Gilbert and Nejadi [39]i n four-point bending. Whilst the beams were not subjected to direct axial compression, their behaviour is representative of that observed in the specimens tested in this study. Beams B1-a and B1-b had a clear span of 350 cm long and a cross-section of 25 9 34 cm. The main flexural bottom reinforcement consisted of two 16 mm bars. The 28-day compressive strength, tensile strength and elastic modulus of concrete were 24.8 MPa, 2.8 MPa and 24.9 GPa, respectively. After moist-curing for 14 days, two sustained concentrated loads were applied at 117 cm from the supports. The loads were 18.6 kN and 11.8 kN for specimens B1-a and B1-b respectively, or approximately 50% and 39% of the ultimate flexural moment at midspan. The time-dependent deflections and cracking of the specimens were recorded for 400 days.
3D solid elements are employed to simulate the hydro-mechanical behavior of the beams in ADINA (see Fig. S2 in Supplementary Material). Due to symmetry, only a quarter of the beams is modelled. A 1D rebar element (truss element) is used to model the reinforcement. Bond-slip behavior between bars and concrete is considered through distributed nonlinear spring elements.
The mechanical models described in previous sections (considering the interaction between drying, damage and creep) are used to assess the timedependent deflections of the beams reported in [39], where standard shrinkage and creep tests were also performed. The shrinkage and creep specimens were exposed to the same environmental, curing and drying conditions as the tested beams. Accordingly, the surface factor f and the hydro-shrinkage coefficient b sh are calculated through inverse analyses using the results of drying shrinkage tests reported in [39]. The parameters for the creep model are defined according to the results of concrete properties and To check the influence of differential dry shrinkage and nonlinear creep, the original shrinkage and creep models in MC2010 are adopted and examined. The damage effect (Sect. 3.1) is also considered in the MC2010 model to simulate the initial cracks due to the first loading. Figure S3 in the Supplementary Material shows the estimated distribution of internal humidity in the core and the surface of the beams at t = 15, 50 and 300 days. The FEA results of time-dependent RH show that the drying depth of the exposed surface is only 1-3 cm in the first 300 days, while the internal RH at the core varies marginally. The non-uniform RH between the surface and the core of the cross section causes a differential drying shrinkage strain, as well as drying-induced stresses. This is inconsistent with the uniform shrinkage model adopted in MC2010, where the differential drying shrinkage strain and stress are neglected. Conversely, the one-way hydro-mechanical coupling analysis approach proposed in this study considers the drying shrinkage and drying-induced stress according to the distribution of the internal humidity on the beams. It is shown that the proposed approach predictions match well the experimental beam deflections for both the instantaneous deflection and time-dependent deflection. This confirms the importance of considering the drying-induced stress in the analysis of longterm deformations of concrete elements.
The proposed model also matches well the crack distribution observed during the tests performed by Gilbert and Nejadi, as shown in Fig. S4 in Supplementary Material. It is shown that MC2010 matches Test [39] This study MC2010 Chong et al. [40] Defection(mm) Time(day) Fig. 12 Comparison of calculated deflection and test results of beams B1-a and B1-b well the cracks in the constant moment region (due to the first loading), but not the cracks outside this region that developed during the period of sustained loads. For instance, cracks 1 and 2 formed at 20 and 51 days after the initial application of the load, respectively. The underestimation of concrete cracking outside the constant moment region explains the underestimation of the long-term deflection given by MC2010 in Fig. 12. These results further confirm that both the applied load and the drying-induced stress contribute to cracking in concrete elements. To further validate the proposed model, Fig. 14 also includes the predictions of the model proposed by Chong et al. [40], which is widely used for the time-dependent analysis of reinforced concrete structures. The figure only includes Beam 1a from Gilbert and Nejad's study because the results of beam B1-b are not available in Ref. [40]. It is shown that the model proposed in this study matches the experimental results better than Chong et al.'s model. Overall, the results in this study show that the proposed numerical approach estimates the behaviour of the concrete specimens with good accuracy. Therefore, it is proposed to use this model to assess time-dependent behaviour and crack development in specimens subjected to eccentric compression taking into account ageing, drying and cracking effects. However, due to the limited test data available for the model calibration and validation, future research is needed to further generalise the model. It should also be mentioned that whilst an isotropic elastic damage model was adopted in this study, the use of an elasticplastic damage criterion could provide better results on damage evolution. Moreover, a one-way coupling method was adopted to determine the interaction between drying, damage and creep. This implies that although the proposed model considers the contribution of concrete drying to surface damage (cracking), it ignores the additional drying due to surface damage, which in turn could speed up moisture evaporation. Further research is needed to clarify this issue.

Conclusions
This study proposed a new creep damage model based on modifications to the linear creep approach given in Model Code 2010. The parameters of the model were determined and calibrated against 12 small-scale concrete specimens subjected to long-term eccentric axial load. The proposed model simulates accurately the development of shrinkage and creep, as well as the interaction between drying, damage and creep. Based on the results of this study, the following conclusions are drawn: 1. The results of shrinkage specimens calculated by the hygro-mechanical simulation show that the creep of concrete has no obvious influence on the development of shrinkage strain for both the sealed or unsealed specimens tested in this study, but it significantly affects the drying-induced stress and the damage distribution. 2. The model indicates that eccentric loading leads to different tensile and compressive creep on the specimens' faces, which contradicts the current design approach that suggests that compressive loads produce similar tensile and compressive creep. 3. The relationship between lateral deflection and load eccentricity in short columns is nonlinear. As the eccentricity increases from 5 to 20 mm, the differential strain (which represents the deflection curvature of the specimen) increases by 5%. 4. The proposed approach matches the experimental results of the tested specimens with good accuracy. The accuracy of the approach is confirmed by further validation with test data available in the literature. Therefore, it is proposed to use this model to assess the time-dependent behaviour and the development of cracking of concrete specimens subjected to eccentric compression and taking into account ageing, drying and cracking effects. Since only a few parameters are needed during calculations, the proposed approach is suitable for practical assessment of long-term behaviour of concrete bridges and buildings.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Human and animal rights This research does not involve human participants and/or animals.
Informed consent Informed consent is not needed as the research does not involve human participants.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.