Tridimensional Long-Term Finite Element Analysis of Reinforced Concrete Structures with Rate-Type Creep Approach

: This paper presents a general procedure for a rate-type creep analysis (based on the use of the continuous retardation spectrum) which avoids the need of recalculating the Kelvin chain stiffness elements at each time step. In this procedure are incorporated three different creep constitutive relations, two recommended by national codes such as the ACI (North-American) and EC2 (European) building codes and one by the RILEM research association. The approximate expressions of the different creep functions with the corresponding Dirichlet series are generated using the continuous retardation spectrum approach based on the Post–Widder formula. The proposed rate-type formulation is implemented into a 3D ﬁnite element code and applied to study the long-term deﬂections of a prestressed concrete bridge built in Romania, which crosses a wide artiﬁcial channel that connects the Danube river to the port of Constanta in the Black Sea.


Introduction
An accurate simulation of creep and shrinkage behavior is necessary for certain types of structures such as long-span prestressed box girders; cable-stayed or arch bridges; large bridges built sequentially in stages by joining parts; nuclear containments and vessels; large gravity, arch or buttress dams; cooling towers; large roof shells; and very tall buildings [1,2]. In this paper reference will be made to prestressed concrete beams, but the model developed is of general application.
Concrete creep can be modeled for service stress state levels (less than half of compressive strength) using the framework of linear visco-elasticity with aging [1]. As a consequence, the principle of superposition holds and the material behavior is uniquely described by a compliance or by a relaxation function. However, in practice, all the reinforced concrete design codes use the compliance (or creep) function to characterize the visco-elastic behavior, because the creep tests are much more common and easy to do than the relaxation tests. They usually describe the compliance function by a suitable formula with several parameters that can be calibrated by fitting experimental data or estimated using empirical formulae that take into account the concrete mix composition, curing conditions and time, member size and shape, and external relative humidity and temperature.
The compliance function expresses the evolution of uniaxial strain over time in a creep test for a unit uniaxial stress. This stress-strain relation can be written as where J(t, τ) is the aging creep function and the times t and τ correspond to the age of concrete starting from the initial setting time. However, the analytical evaluation of the integral in Equation (1) is possible only for simple models and simple stress histories. For general applications in finite element codes, the numerical integration of Equation (1) requires the storage of the entire stress history at each integration point of each finite element and, as a consequence, the evaluation of the integral requires an extensive memory allocation and an increasing number of calculations for each time step as the observation time t progresses. Several methods have been presented in the literature to simplify the calculation of creep strain under time history loading as in Equation (1), such as the effective modulus method [3], rate of creep method [4], the ageing coefficient method (AAEM method) [5] utilized in many applications, among others [6,7], the approach based on the aging linear viscoelastic theory [8], and, recently, the parallel creep method [9] that is extended in [10,11] for a general age-dependent constitutive law. The computational cost of the integration of Equation (1) can be substantially reduced by replacing the integral stress-strain relation with a differential one (rate-type formulation, first proposed in [4]). This approach is based on the approximation of the compliance function by a Dirichlet series (i.e., a sum of exponentials) corresponding to a Kelvin rheological chain whose units are described by a differential equation that can be easily integrated in a step-by-step manner. This procedure transforms the original integral approach into a differential (rate-type) approach in which it is not necessary to store the entire load history but only a limited and fixed number of internal (history) variables. It also gives rise to a number of numerical calculations constant for each time step, independently of the length of the time interval considered.
In the rate-type approach the crucial point is the procedure adopted to convert the continuous creep function into its Dirichlet series approximation, due to the difficulties which arise in the identification of the series coefficients for aging formulations (typically used in codes and recommendations). One method was proposed in [12,13] for the study of the deflection behavior of the Koror-Babeldaob Bridge in Palau through a 3D finite element procedure. In their approach, to deal with the aging properties of the compliance function in each time step and at each integration point, use was made of the Widder's formula to convert the aging compliance function into a continuous retardation spectrum for the current age of concrete. A discretization of the spectrum would then yield the current elastic moduli of the Kelvin units.
The purpose of this paper is to formulate a general procedure, also based on the use of the continuous retardation spectrum, capable of readily converting the integral creep problem into a rate-type one without the recalculation of the Kelvin chain stiffness elements at each time step. This procedure is then particularized with reference to three different creep constitutive relations, two important because recommended by professional associations in USA and Europe (ACI and CEB-Euro-Code) and one important for its diffusion in the research community, the RILEM B3 model.
The procedure has then been implemented in a three-dimensional finite element model for bridge design, capable of taking into account both the different shrinkage and drying creep properties of the various parts of the bridge cross-section and the shear lag effect, both of which cannot be captured by the classical beam theory commonly used in bridge design. In particular the shear lag effect influences drastically the accuracy of the calculation of the prestress loss and of the bridge deflection. As reported in [12,13] the shear lag occurs in different ways: In the transmission of vertical shear force due to vertical reaction at the pier and in the transmission of the concentrated forces of tendon anchors into the horizontal slabs and the vertical walls of the box.
In order to validate the proposed finite element formulation, two examples are considered, a prestressed simply supported beam and a prestressed cantilever box girder, for which the results of the 3D finite element analyses are compared with the one-dimensional calculations. Then the long-term behavior of a real bridge of "balanced cantilever girder" type, characterized by a central 155 m span, with two side symmetric 77.5 m spans, has been analyzed (and has been the motivation of this work). The bridge belongs to the new Medgidia-Constanta Motorway crossing a wide artificial channel that connects the Danube River to the port of Constanta in the Black Sea. It holds the record span-length in Romania for prestressed concrete deck type and has required for its construction concrete volumes of 14,700 m 3 for foundations and 7700 m 3 for piers and decks, total reinforcing steel weight of 2600 t, and total post-tensioning tendons weight of 390 t. For this bridge, the predictions of the various creep formulations in terms of deflections, prestress loss and stress state in the upper and lower slabs of the cross section are compared.
It must be remarked that an exact numerical analysis of this type of structures should be performed through a multi-physics time-dependent approach based on a hygro-thermo-chemical-mechanical model, as in [14][15][16], that could be useful in predicting the concrete time-dependent response. However, the multi-physics modeling is out of the scope of the manuscript and also the dimensions of this type of structures prevent its use for now. For general constitutive formulations that starts from early-ages, an interested reader can also refer to the formulations presented in [17][18][19][20][21].

Integral Formulation
Using the principle of superposition and knowing the compliance function J(t, t 0 ) it is possible to determine the uniaxial strain evolution from any general uniaxial stress history, σ(t), in a integral-type form where t is the current time, t 0 is the historic time and t − t 0 is the elapsed time and the integrals are understood in the Stieltjes sense, so that they can be evaluated even for discontinuous stress and strain evolutions. Equation (2) where matrix G is constant over time The value ν = 0.18 will be used in all numerical simulations presented in this paper. Let us express Equation (2) in the incremental discrete form for finite element applications. In the time interval ∆t = t i+1 − t i we can define the stress and strain increments: Using Equation (2), the strain increment in the time step ∆t is and rearranged as Assuming a linear stress variation in the time interval, i.e.,σ = ∆σ/∆t, we can rewrite Equation (7) as To calculate the second expression in Equation (9) one needs to know the entire load history in each material point and in each time interval. This requires to store during the analysis a huge amount of information that increases drastically with the length of time interval considered. This makes the integral approach for creep prohibitive in a finite element program for long-term analysis, while it might be utilized for short-term processes such as early-age analyses, although this is not the most efficient procedure.

Rate-Type Creep Law
As a remedy, it was showed [1,22] that through the approximation of the creep function as a sum of negative exponentials (i.e., Dirichlet or Prony series) it is possible to convert the integral expression (in Equation (2)) into a set of linear differential equations, which turn out to be the governing equations of the well-known Kelvin chain rheological model with aging spring moduli E j (t) and dashpot viscosities η j (t). These functions can be also identified as the coefficients of the Dirichlet series expansion of J(t, t 0 ). Using this approach, a convenient approximation of the J(t, t 0 ) function can be written as where E j (t 0 ) are aging material parameters function of time t 0 being identified from the fitting of the given function J(t, t 0 ) at any given fixed t 0 . A simple and convenient choice for the other material parameters τ j , called the retardation times, is to assume them constant. However, they must be chosen suitably with a not too large spacing in the log scale in order to approximate adequately the J(t, t 0 ) function in all time duration of interest during a the structural analysis. The function J(t, t 0 ) is the design compliance function being defined by some formula arising from the adopted code or standard. The term for instantaneous deformation 1/E 0 (t 0 ) can be considered included as the first term of the summation with a retardation time extremely small (τ 0 ∼ 0). Substituting the approximation of the compliance function as in Equation (10) into the expressions in Equation (9) one gets That can be rewritten as where the integrals have been calculated using the mid point rule with t * = 0.5(t i + t i+1 ) and * j (t i ) can be considered as internal variables that must be stored for the integration points of each finite element. Those internal variables can be updated according the following recurrence formula (obtained by calculating σ(t i+i ) assuming a linear stress variation in the time step and constant time step increments ∆t) * It must be noted that the number of the internal variables for the rate-type formulation (Equations (8), (11), (13), and (14)) is now fixed and limited to N, in contrast to the integral formulation (Equations (8) and (9)) that requires the storage of the entire load history for each material point, which increases drastically with the final time of the analysis (for long-term behavior is typically many decades).

Aging Kelvin Chain
It is worth showing that the quasi-elastic constitutive law in Equations (8), (11), and (13) that has been just derived is fully equivalent to the solution of a rheological model consisting of a aging Kelvin chain with springs of stiffness E j (t) in parallel with a dashpot with viscosities τ j E j (t). To simplify the mathematical derivation, let us admit that in each time interval ∆t = t i+1 − t i of the loading history the spring stiffnesses is constant and equal to E j (t * ) where t i < t * < t i+1 . Under this assumption the constitutive law of the aging Kelvin element in a time interval is given by the following first order differential equation An effective numerical integration of Equation (15) can be done by virtue of the so-called exponential algorithm, which makes possible the use of increasing time steps with the same accuracy and numerical stability [1]. The exponential algorithm assumes that the stress varies linearly in the time interval ∆t = t i+1 − t i so that the differential equation in (15) can be integrated exactly. The linear stress variation can be assumed as where ∆σ is the stress increment over the time step ∆t. With this stress variation the general solution of the differential equation in (15) can written as The integration constants A and B of the particular integral can be obtained from the substitution of Equation (17) into (15) while C is calculated from the initial condition in each time step, = (t i ). This yields the strain, (t i+1 ), at the end the time step (t = t i+1 = t i + ∆t) as So the strain increment is Considering now a rheological model chain of N Kelvin elements in series with a spring with stiffness E 0 (t * ), the total strain increment is given by Substituting in this expression ∆ j from Equation (19) and ∆ 0 = ∆σ/E 0 one gets the stress increment in the time step from t i to t i + ∆t as where (22) in which * j (t i ) are the internal variables that can be updated according the recurrence formula in Equation (14).
Since E * = E and σ * = σ the two approaches are fully equivalent. This means that the approximation of an aging compliance function through a Dirichlet (or Prony) series transforms the classical Volterra integral equation of creep (Equation (2)) into a rate-type formulation governed by an aging Kelvin chain with spring moduli E j (t * ) and viscosities τ j E j (t * ) obtained directly from the coefficient of the Dirichlet series approximation of the compliance function. At this point, to generalized the application of the rate-type approach, there is the need for a robust and reliable procedure capable of determining the coefficients of the Dirichlet series that approximate a given compliance function.
As pointed out in [23] the retardation times must be chosen "a priori" because their calculation from experimental data can give an ill-conditioned equation system. A suitable choice is τ n = τ 1 10 n−1 with n = 1, 2, 3...N (23) which means that the retardation times are equally spaced in a logarithmic scale and this gives smooth enough creep curves. Each of these times is representative of one order of magnitude, covering the interval from τ n / √ 10 to τ n √ 10. For a general analysis, 10 Kelvin units, i.e., N = 10, are enough to consider a wide spectrum time, i.e., from 10 −4 days to 10 +5 days.

Non-Aging Kelvin Chain
For non-aging material the expressions previously derived are simplified, since the equivalence is imposed with a non-aging Kelvin element (with spring stiffnesses constant). Moreover, the approximation of a non-aging compliance function, J(t − t 0 ), through a Dirichlet (or Prony) series, is obtained using constant moduli E j as The values of the coefficients E j can be determined by the best fitting a given non-aging compliance function J(t − t 0 ) either using a minimization algorithm (Least Squares or Levenberg-Marquardt algorithm) or passing through the continuous retardation spectrum [24], which is a more general approach that is summarized in the following.
As shown in Equation (24), a Kelvin chain model with N units gives a non-aging creep compliance function C(t − t 0 ) which can be approximated by a Dirichlet series as [22] where A n = 1/E n , t = age of concrete, and t 0 = age of concrete at the moment of loading. To deal with a general procedure and to avoid some weak points of other approaches [22,25], the continuous Kelvin chain model with infinite units (continuous spectrum), in which the retardation times are infinitely close, is used [24]. Passing through the continuous spectrum, the discrete spectrum can be obtained by discretizing the continuous one. The creep compliance function C(t − t 0 ) may be approximated in a continuous form as where ξ = t − t 0 and A(τ) is the continuous retardation spectrum, which has the same meaning in the logarithmic time scale as A n in the real time scale. Following the method developed by [26] and setting τ = 1/χ we have and The previous Equation (28) can be rewritten as where f (ξ) if the Laplace transform of the function A(χ −1 )χ −1 and f (0) is a constant. Using the inversion formula of Widder [27] the function A(χ −1 )χ −1 can be obtained asymptotically as where f (k) is the kth derivative of the function f . Remembering that f (ξ) = f (0) − C(ξ) and χ = 1/τ, we obtain the continuous retardation spectrum The approximate spectrum of order k is obtained by assuming a finite value of k in Equation (31). For practical purpose an approximate spectrum of third order (k = 3) may be used with enough accuracy [24,28], i.e., Decomposing the integral by a finite sum over finite time intervals given by Equation (23), the Equation (26) can be rewritten in the time of interest, i.e., 0 < t < τ N √ 10, as The integrals can be evaluated using the n-point Gaussian quadrature rule. Using the one-point quadrature rule in the intervals ∆ (lnτ n ), which with Equation (23) are given by ln10∆ (logτ n ) = ln10log10 = ln10, the Equation (26) can be written as (34) and with the coefficients of the Dirichlet series given by Comparing Equation (34) with Equation (24) we observe that the constant term A 0 must be added to 1/E 0 and that 1/E n = A(τ n )ln10.
However, as shown by Jirásek and Havlásek [28], the accuracy of the Dirichlet series which approximates a compliance function expressed by Equations (34) and (35) is not always good and it can often be substantially increased by appropriate modifications of the discrete retardation times adopted in the Dirichlet series. Jirásek and Havlásek [28] demonstrated that in many cases the accuracy of the Dirichlet series approximation cannot be increased by increasing the derivative order in Equation (32), but rather by the adjustment of the discrete retardation times, applied after the evaluation of the compliance coefficients. With this adjustment the Dirichlet series approximation is expressed as The expressions of such adjustments, ψ j , are presented in the next section for each creep function considered in this study.
The presented approach for non-aging compliance function can be extended for aging creep function of the different codes and standards as it is presented in the next section.

Numerical Algorithm
As already reported in the literature [17,22,29,30], the finite element analysis of long term behavior with creep is much more efficient if a rate-type approach is used instead of an integral-type form. A Kelvin chain, as well as a Maxwell chain, arrangements of springs and dashpots can described the most general creep behavior [29,30]. Since the material constitutive law of concrete is typically based on the assumption of the strain additivity, a Kelvin chain is more convenient [16,31] than the Maxwell chain. Following the original idea of Bažant [32], the structural creep problem can be reduced to a sequence of elastic finite element analyses using an elastic stress-strain relation with inelastic strain, i.e, step-by-step linear elastic analysis for each time step. This means that in each time step the rheologic model can be considered non-aging and, consequently, its spring moduli and viscosity are constant and updated only at the beginning of the time step.
As shown in the previous section, the incremental quasi-elastic stress-strain relation suitable for a general finite element program can be written as (Equation (8) or Equation (21)) where E * and σ * (t i ) have been derived in the previous section (see Equations (11) and (13) or Equation (22)) and ∆ 0 is the inelastic strain increment in each time step, such as the shrinkage strain or thermal strain [16]. This incremental stress-strain relation represent the quasi-elastic response of a non-aging kelvin chain that approximates a non-aging creep function J(t, t 0 ) = C(t − t 0 ). The formula of the compliance functions generally adopted in standards or codes can be put in the form where d 0 (t 0 ) indicates the instantaneous elastic strain caused by a unit applied stress at the time t 0 and c 0 C 1 (t 0 )C(t − t 0 ) represents the creep deformation, which is expressed as the product of a constant c 0 , an aging term, C 1 (t 0 ), and a non-aging term, C(t − t 0 ). The compliance function in Equation (38) can be approximated by the following Dirichlet series where the coefficients E j are identified from the Dirichelt approximation of C(t − t 0 ) using the expression in Equation (35). In this case, the general incremental quasi-elastic stress-strain relation is still given by Equation (37) with the following expression for E * and σ * (t i ) where t can be taken as the time at the middle of the time step, t i + ∆t i /2. When a specific creep formulation is considered, from the expression of C(t − t 0 ) the coefficients E j must be calculated first and then, using the specified expressions for d 0 , C 1 , and c 0 , one can calculate the value of E * and σ * (t i ) to be utilized in the constitutive law of Equation (37). The formulation proposed for the rate-type creep analysis is based on the continuous retardation spectrum of the adopted constitutive relation for the compliance function. This spectrum is derived below for the most significant creep models.

EuroCode 2 Model
The EuroCode 2 (EC2) expresses the creep behavior through the creep coefficient, φ(t, t 0 ), so that the compliance function may be expressed as where E(t 0 ) is the modulus of elasticity at loading age t 0 . Comparing the compliance function in Equation (38) to the expression in Equation (42), we have Considering the EuroCode 2 formulation reported in Appendix A, the previous expressions in Equation (43) can be rewritten as The approximate continuum retardation spectrum of third order (k = 3 in Equation (31)) is given by with A 0 = 3.96 × 10 −8 0 for β H = 600. The optimum adjustment factors, ψ j , in Equation (36) for the best approximation of the EC2 creep function are given by The quality of the approximation of the EC2 compliance function by Dirichlet series is shown in Figure 1a for different t 0 and with β H = 600, with negligible relative errors (≤ 1%, Figure 1b). The inelastic strain incremental for the EuroCode 2 model due to shrinkage in each time step is given by

Time [days]
where the expression of cd (drying shrinkage) and ca (autogenous shrinkage) are reported in Equations (A16) and (A18), respectively.

ACI Model
According to the ACI-209R-92 code provisions, the compliance function is expressed as where E(t 0 ) is the modulus of elasticity at loading age t 0 and the creep coefficient φ(t, t 0 ) is given in Equation (A25). Comparing the compliance functions in Equation (38) to the one in Equation (48), we have Recalling the ACI-209R-82 formulation, reported in Appendix B, the previous expressions in Equation (49) can be written as and The approximate continuum retardation spectrum of third order (k = 3 in Equation (31)) is given by with A 0 = 1.28 × 10 −9 0 for ψ = 0.6 and d = 10. The optimum adjustment factors, ψ j , in Equation (36) for the best approximation of the ACI creep function are given by The inelastic strain incremental for the ACI model due to shrinkage in each time step is given by where the formula of sh (t i+1 , t c ) can be found in Equation (A20).

B3 Model (RILEM)
According to the B3 model developed by Bažant and coworkers at Northwestern University [33] and recommended by RILEM, the compliance function (Equation (A34) of Appendix C) is expressed as Comparing the compliance functions in Equation (38) to the previous one in Equation (54), one gets Since the term C 0 (t, t 0 ) needs a specific different calculation respect to the previous formulations, it does not appear in Equation (55). Recalling the B3 formulation, reported in Appendix C, the previous expressions in Equations (55) and (54) can be written as The function C 0 (t, t 0 ) describes basic creep by a log-power law with aging incorporated through the solidification theory [22,34] with an additional logarithmic term that reflects viscous flow. This term can be described by a dashpot with age-dependent viscosity and can be treated directly in the rate form, without the need to construct a Dirichlet series approximating its compliance function. Therefore, to describe the basic creep of B3 model, only the spectrum of the log-power Its third order approximation (k = 3 in Equation (31)) with the continuous retardation spectrum (Equation (32)) is given by with A 0 = 0.29209 (from Equation (35)) for τ 1 = 10 −4 days. The quality of the approximation of the B3 compliance function by Dirichlet series with the coefficient obtained from Equation (58) is shown in Figure 3a for different t 0 and with q 3 = 1.5 * 10 −5 . The approximated curves exhibit a maximum error, plotted in Figure 3b, below 1% with the adjustment factors ψ j = 1.2 for j = 1, 2, 3...9 and ψ 10 = 0. For the numerical calculation of the basic creep function of the B3 model, C 0 (t, t 0 ) in Equation (56), reference must be made to the solidification theory [22,34] for which the visco-elastic strain rate is given by˙ Substituting in the second expression of Equation (59) the Dirichlet approximation of the compliance function, Φ, expressed through the coefficient A(τ) given by Equation (58) we have where the γ n are now the internal variables that can be calculated with the following recursive formula, obtained by evaluating ∆γ n = γ n (t i+1 ) − γ n (t i ) with the second expression in Equation (60) and assuming a linear variation of the stress increment in the time step from t i to t i+1 = t i + ∆t Rewriting the first expression of Equation (59) in discrete form, ∆ = (q 2 √ 1/t + q 3 )∆γ, and substituting the ∆γ obtained from Equations (61) and (60), one gets the contribution of the basic creep, ∆σ = E * b (∆ − * b ), which must be added to the stress-strain relations in Section 2.5, in the following form where t = t i + ∆t i /2, γ n (t i ) are the internal variable obtained from Equation (61), and the last term represents an additional logarithmic term that reflects the purely viscous flow. Those two terms in Equations (62) and (63) must be added to the expression in Equations (40) and (41), respectively. In addition to the basic creep, in the model B3 there is a separate drying creep term C d (t, t , t 0 ) (see Equation (A37)), which is capable of reproducing the Pickett effect. Like the shrinkage, the drying creep term is bounded and depends on the humidity and the cross section thickness. Without losing generality we assume t = t 0 (otherwise a constant term should be added to the compliance function) and comparing the compliance function in Equation (38) to the expression of B3 drying creep compliance we have

Time [days]
The approximate third order continuum retardation spectrum of this drying creep function can be expressed as (35) with τ 1 = 10 −4 days, τ sh = 3600, and h = 0.6). The reliability of the approximation of the B3 drying compliance function by Dirichlet series with the coefficient obtained from Equation (65) is shown in Figure 3c for different t 0 and with q 5 = 10 −4 , τ sh = 3600, and h = 0.6. The approximated curves in the figure exhibit a maximum error below 2% as showed in Figure 3d. To have that lever of error (certainly admissible for practical applications) the terms of the Dirichlet series are increased up to 13, with a denser set of retardation times with retardation time interval of 10 0.75 . The Dirichlet approximation of the B3 drying creep function reads where N 1 = 13, τ n = τ n−1 10 0.75 with τ 1 = 10 −4 days, A d (τ n ) from Equation (65), and with t = t 0 , i.e., drying and loading act simultaneously at time t 0 , and ψ n = 1.3. The hygral strain incremental for B3 model in each time step is given by where the equations for sh (t i+1 , t c ) can be found in Equations (A28)-(A31).

Numerical Validation of the Finite Element Model
In this section two numerical applications of the proposed approach are presented. The previously illustrated constitutive relations have been implemented in three-dimensional finite element program, because only this approach can easily analyze cases with complicated loading history, with different distribution in the structure of the material properties, with different construction phases, with cross-section in which the inhomogeneous effects of the diffusion processes are not negligible, with changes of the structural configuration, and with the shear-lag effect. This effect is characterized by out-of-plane warping of cross-sections where there is high shear force and by a nonlinear stress distribution, that are neglected by the classical beam-type analysis. The reason for taking it into account is that the shear lag due to the self weight is stronger than the one generated by the prestress and since the total deflection is a small difference of the downward deflection due to self-weight and the the upward deflection due to prestress, a small error in only one contribution produces a larger error in the total deflection. The finite element program employed for all the numerical simulations presented in this work is a Fortran code written by the first author in which the implicit time (real time) integration is performed using the Newton-Raphson method with a constant global stiffness matrix obtained with the value of the Young modulus at 28 days.
All the calculations and the numerical simulations are done assuming concrete in the uncracked stage, which can be justified by the limited tensile stresses in concrete. However, the proposed formulation can be easily extended to include non-linear behavior and cracking, for instance see [13,16,30,35] for mesoscale formulation. Moreover, the relaxation of the steel prestress bars or strands is not considered in the analyses which follow.

Numerical Simulation of a Prestressed Beam with I-Shaped Solid Cross-Section
The first example concerns a simply supported prestressed concrete beam with a constant I-shaped cross-section and with a single equivalent prestressing strand with steel diameter of 36 mm. The cross-section geometry of the considered beam is showed in Figure 4a, the boundary and loading conditions in Figure 4b. In order to validate the numerical formulation of the rate type creep model and its implementation into a three-dimensional finite element code, the mid-span deflection in time obtained from the finite element analysis is compared with its analytical calculation resulting from the Effective Modulus Method (EMM) [36] applied to the type of load history considered (see Figure 4c). Using the EMM the mid-span deflection can be calculated with the Principle of Virtual Work as where E e f (t, t 0 ) = 1/J(t, t 0 ) is the effective modulus calculated from the creep function J(t, t 0 ), G e f (t, t 0 ) = E e f (t, t 0 )/2(1 + ν) is the effective shear modulus, I(x, t, t 0 ) is the cross-section moment of inertia, A(x, t, t 0 ) is the cross-section area, and κ(x, t, t 0 ) is the shear factor [37]. The geometric properties depend on time through the homogenization coefficient n = E s /E e f (t, t 0 ). Assuming a constant cross-section and substituting the bending moment and the shear corresponding to the external uniform distributed load of 30 kN/m applied at time t 0 (sketched in Figure 4c) the mid-span deflection, f q (t, t 0 ), is given by The prestressing tendon is modeled using beam finite elements connected rigidly to the nodes of the three-dimensional mesh (no slip). No regular reinforcing steel bars are considered in this application. The prestress force, N p , is applied in the tendon by assigning an initial equivalent thermal deformation which generates a stress which accounts for the initial elastic loss. The initial prestress force is of 1297 kN that is equal to prestressing force of N p = 1250 kN after the initial elastic loss. Substituting the bending moment (M p = N p e with e = the eccentricity of the tendon) the mid-span deflection, f p (t, t 0 ), is given by It has been considered no relaxation of the steel and no shrinkage of concrete in this first example. The creep model adopted is the EC2 creep model assuming h = 0.5, h 0 = 20.445 m 2 /5.397 m = 164.9 mm, f ck = 55 MPa, f cm = 63 MPa, t 0 = 14 days. The comparison between the mid-span deflection of the beam as calculated using the beam theory with the EMM (Equations (69) and (70)) and the deflection obtained with the 3D finite element analysis with rate type formulation is presented in Figure 5 for different loading conditions showing a coincidence of the calculated deflections with the two methods. The excellent agreement is also conformed by the evolution of the normal stresses and by the force in the steel tendon, as shown in the following.  The normal stresses in the concrete cross-section can be calculated as in which y is the distance from the centroid of the homogenized section and M is the bending moment. In the following, tensile stresses are considered positive. The force in the prestress steel tendon is also compared in Figure 6a with the analytical value that can be calculated using the age-adjusted effective modulus approach, given by where E aae f (t, t 0 ) = E c (t 0 )/(1 + χ(t, t 0 )ϕ(t, t 0 )) is the age-adjusted effective modulus with χ(t, t 0 ) 0.8, σ c is the stress in concrete at the level of the tendon and ∆σ c = σ c (t) − σ c (t 0 ) is its variation in time, calculated using Equation (71) with the appropriate y. Usually, in the design procedure and also in the code recommendations, only the first term in Equation (72)  In addition also the capability of the proposed approach for complex load histories is considered. Complex load histories that are reported on the left of Figure 7 are adopted and the corresponding numerical solution is compared with the exact integral formula In addition in Figure 7 is also reported the approximate solution give by the AAEM

Numerical Simulation of a Prestressed Box Girder
The second example deals with a cantilever prestressed concrete box girder with variable cross-section and with two equivalent prestressing bars with steel diameter of 36 mm. In Figure 8a,b are shown the cross-section geometries at the fixed and at the free end, respectively, of the considered girder. The boundary and loading conditions are presented in Figure 8c. For this type of geometry, the drying process, which drives shrinkage and drying creep, causes nonhomogeneity of shrinkage and creep properties throughout the cross-section that are nonsymmetric with respect to neutral axis. This phenomenon can not be captured by a 1D beam analysis and it has often been one major cause of gross mispredictions of long-time deflections of structures [38,39]. Therefore, only using a 3D analysis one can capture the different properties in the cross-section.
Moreover, for this application, the 3D finite element numerical solution of the rate type creep model is compared with the analytical calculation using the beam theory (with the above-mentioned limitations) and the Effective Modulus Method (EMM) using the load history plotted in Figure 8c). Using the EMM the free end displacement can be easily calculated with the Principle of Virtual Work as in Equation (68) where the effective modulus and effective shear modulus are calculated as in previous example. The moment of inertia, I(x, t, t 0 ), the area, A(x, t, t 0 ), and the shear factor, χ(x, t, t 0 ), are calculated according to the current the cross-section, which varies from the cross-section at the fixed end in Figure 8a to the cross-section at the free end. No relaxation of the steel and no shrinkage of concrete is adopted for this second example. The creep model adopted is the EC2 creep model assuming h = 0.6, f ck = 55 MPa, f cm = 63 MPa, t 0 = 14 days, and a cross-section value of h 0 which is assumed varying lineally from 218.6 mm at the fixed end to 155.56 mm at the free end. The comparison between the displacements of the beam as calculated using the beam theory with the EMM (Equations (69) and (70)) and the displacements obtained with the 3D finite element analysis with rate type formulation is presented in Figure 9 for different loading conditions showing a perfect coincidence of the calculated deflections with the two methods. The excellent agreement is also confirmed by the evolution of the normal stresses and by the force in the steel tendon. The normal stresses in the concrete cross-section can be calculated using the expression in Equation (71) with the correction factors proposed by [40] as k α1 = 1 + 4 tan(α) 2 and k α2 = 1 − 4 tan(α) 2 for the straight and tapered side, respectively. In Figure 10a the time evolution of the axial force in the steel tendon is shown and compared with the analytical solution obtained from Equation (72) presenting a very good agreement between them. Figure 10b-d report the evolution of the normal stresses obtained with the analytical (Equation (71) with the above correction factors) and the numerical solution which basically correspond with each other.

Numerical Simulation of the Long-Term Behavior of a Bridge
The bridge is of the "balanced cantilever girder" type, built with segments cast from piers with mobile equipment, and it is characterized by a central 155 m span (the longest span in Romania for prestressed concrete box girder bridges), with two side symmetric spans of 77.5 m (see Figure 11). The deck has a varying depth provided by a curved soffit of parabolic shape. The depth of cross section varies from maximum value of 10 m at the pier axis, to a minimum value of 2.4 m, at mid span and at the abutment supports. The upper slab is 14.75 m wide and transversally inclined of 2.5% (see Figure 12).  The layout of the post tensioned internal tendons (total length 10 km/single deck), follows the upper and lower slabs (see Figure 12). The upper tendons counterbalance the action of the free cantilever under gravity loads. The lower tendons are located along bottom slab and anchored to reinforced concrete internal blisters. The layout is symmetric with respect to the midlength of the central span both for the central and the side spans.
The bridge deck is supported by two piers and two abutments through two seismic bearings at each support. The bearings are of the "friction pendulum" type, which develop friction forces both in static conditions, due to static forces and small displacements, and dynamic condition providing dissipation. Under seismic loads (moderate for that construction site) the whole concrete structure develops only elastic behavior, because dissipation and lengthening of the natural period of the structure are provided by the seismic bearings.
Piers 1 and 2 are characterized by same shape but different height, respectively, 17.4 m and 16.15 m. The pier has a hollow rectangular cross section with 8 m transverse and 6 m longitudinal external dimensions and a massive head capital at the top. The piers base sections are connected to a reinforced concrete massive rectangular footing erected on a ring of diaphragm walls capable of transferring to a deep ground level the forces coming from the superstructure. Both abutments are spill-through type. The beam seat is supported by a number of shear walls aligned with longitudinal deck axis. Each shear wall is founded on deep diaphragm walls that transfer to the deep ground level load due to superstructure and earth pressures.
The design of the viaduct was carried out taking into consideration many advanced issues, including creep and shrinkage deformations and seismic behavior. A finite element numerical model of the prestressed concrete deck, capable of simulating more than thirty construction phases was implemented with specific software. The long-term behavior was studied accurately by means of an experimental/numerical procedure which is presented in this paper. Because of piers low ductility capacity, the the seismic protection was achieved by integral isolation through "friction pendulum" devices that develop dissipation by friction mechanism. These type of seismic bearings show many advantages in term of cost/performance in comparison with traditional high dissipation rubber bearings, mainly if they sustain high level of axial loads (as it is in this case, 40,000 kN/each pier bearing, 170 cm diameter). Dynamic tests of "friction pendulum" bearing devices have been performed at SRMD lab, University of San Diego (USA) [41]. The detail project of foundation design was developed in collaboration with Technodata (Naples). Design assistance in the construction field, was developed with constant support by Italrom (Bucharest, affiliate company of Lombardi-Reico).
The concrete used for the bridge construction has the following mix composition: Cement CEM III A 42.5 N-LH; polycarboxylic superplasticiser; natural calcareous aggregates with maximum aggregate size of 25 mm with a water/cement ratio of 0.37 and superplasticizer/cement ratio of 1% (by weight).

Structural Effects of Long-Term Deformations
The analysis of the structural behavior has been performed through a three-dimensional finite program applied, because of symmetry, to only one half of the bridge. The box girder of this bridge is considered a thick shell which is discretized by brick (8-node three-dimensional) finite elements and the generated mesh is shown in Figure 13. The prestressing tendons (see Figure 13) are modeled through beam elements connected rigidly to the nodes of the three-dimensional mesh (slip is allowed along each tendons). Non-prestressed reinforcing steel bars are not considered in the present discretization. The fineness of the mesh has been validated by checking that a finer mesh with the double of hexahedral elements would yield only a negligible improvement of the computed elastic deflections. The integration algorithm described in Section 2 reduces the creep problem to a series of elastic structural analyses in each time step. A finite element code developed by the first author was used to carry out the step-by-step elastic analysis. The introduction of prestress in each tendon was done by assigning the initial stress to the tendons through an equivalent thermal deformation. To reproduce the time sequence of segmental construction, the finite elements corresponding to different segments were activated at the time of their casting, which also introduced a different age of concrete for each segment. Each segment was activated taking into consideration the correction of level adopted during the construction procedure as reported in Table 1. These corrections have been introduced in the construction process for taking into account the cumulative deflection at the free end of the cantilever element and for obtaining the demanded level of the road. In other words, those corrections are introduced in order to compensate for the initial creep deflection at the free end of the cantilever beam, due also to the weight of the construction apparatus (see the first picture in Figure 12). During the deck construction the different type of post-tension tendons/bars (see Figure 12) were activated as described in the following. The free cantilever is built from the pier segment by adding up to 15 side segments, each one stabilized by the activation of two upper tendons, for a total of 32 upper cables. At the end of free cantilever construction, abutment segments and key center segment are connected by lower cables and pre-stressed bars for a total of 88 post-tension tendons/bars. Table 2 reports all the features of post-tension tendons and bars. The tendons of each segment were prestressed 2 days after the segment casting. The friction caused by curvature was neglected while the wobble friction was considered. The initial prestress diminishes also because of the relaxation of steel. However, this effect is not considered in the 3D finite element analyses, although it can be implemented in the numerical algorithm just imposing a reduction of the initial prestress in accordance with some formula taken from a code (for instance the EuroCode2 formula). The creep phenomenon is strongly influenced by the temperature value. In the numerical simulations however the effect of the temperature increase due to solar heating of the top slab is neglected as well as any effect of the environmental temperature variation on the concrete properties. Tensile cracking has not been taken into account because no significant tensile stresses arise in the numerical simulations and, therefore, the mechanical behavior can be assumed as linear elastic. The nonlinear effects for creep have not been considered because the stress levels are always below 0.4 f cm .
In the numerical analysis the self-weight has been applied through a volume force (γ = 25 kN/m 3 ). During the construction phases the weight of the mobile launching wagon were considered by means of a vertical force of 750 kN acting at the head of the overhang, and therefore, with a bending moment dependent on the size of a segment. At the end of the construction, an additional permanent load due to the pavement and superstructure equal to 35.8 kN/m is applied on the top surface of the upper slab.
The input parameters required for creep and shrinkage analysis depend upon the model adopted, however, usually they are: The 28-day elastic modulus of concrete, E c , or the required design strength f c28 , from which E c can be estimated; the starting age of drying, taken as t = 2 days, which corresponds to the segmental erection cycle; the average environmental humidity h = 0.6; the effective thickness of cross-sections, D = 2V/S (V/S is the volume/surface ratio) see Tables 3 and 4; some concrete composition parameters (especially for model B3). The age at the start of drying is taken here as 2 days, which is the time of formwork removing in the segmental erection cycle of 7 days. For the EC2 model the following values have been utilized: f cm = 62.90 MPa, f ck = 54.90 MPa, and cement of class N so α = 0, α ds1 = 4, α ds2 = 0.12 (see the Appendix A for all the details). For the ACI model the 28-day mean compressive strength of f cm28 = 62.90 MPa is adopted with the coefficients reported in the Appendix B. For the B3 model using in Equation (A39) the 28-day mean cylinder compressive strength of f cm28 = 62.90 MPa and the concrete composition we have q 1 = 1.598 × 10 −5 , q 2 = 9.248 × 10 −5 , q 3 = 5.026 × 10 −7 , q 4 = 7.107 × 10 −6 , α 1 = 1.1, α 2 = 1.2, k s = 1.0, s∞ = 6.87 × 10 −4 .

Long-Term Variation of Stress and Deformation States
The All these responses have been computed with the same finite element program and the same step-by-step time integration algorithm as presented in Section 2. It must be remarked that the relaxation of the steel tendons is not considered in the present analysis to emphasize only the effect of creep and shrinkage on the long-term behavior of the bridge. Figure 14 presents the mid-span deflection, both in linear and logarithmic time scales, in which the time is measured from the end of construction stage. It must be remarked that the deflection is also evaluated with reference to the configuration existing at the end of the construction stage, i.e., this deflection is not the absolute value of the mid-span deflection. The bridge deflection is calculated as the difference between two large numbers (both affected by some uncertainties): The downward deflection due to the dead load and the upward deflection due to prestress. The numerical analysis takes into account the effect of the differences in thickness of the slabs and webs on their drying rates through a different effective thickness (D = 2V/S) in the cross section (see Table 3). If the shrinkage and the drying creep compliance are considered to be uniform over the cross section the overall effective thickness, D = 2V/S of the whole cross section as reported in Table 4, can be used. The use of non-uniform creep and shrinkage properties throughout the cross section allows to simulate the initial upward deflection due to differential shrinkage and differential drying creep, as shown by Figure 14 for all the creep models considered in the first years after the end of the bridge construction. The predicted deflections at about 80 years (∼30,000 days after span closing) are quite different for shape and values depending on the adopted model. The Model B3 predicts the greater value of the deflection, ≈40 mm, while using the ACI model the deflection is of few millimeters upward and for the EC2 model it reaches an asymptotic value of ≈20 mm. The ACI deflection curve presents a shape that is rather different from those of the other models. It gives a deflection growth during the first years which is not compensated at a later time by the creep deformation and, as a result, the deflection tends to a small positive value. The EC2 deflection curve presents a smaller influence of the shrinkage effect at the beginning and then predicts a deflection which tends to an asymptotic value. On the contrary the B3 model deflection curve shows a large influence of the shrinkage for the first 6/7 years with an upward deflection, but then presents an increasing downward deflection with a reduced rate and with no tendency to an asymptotic value.   Figure 15 presents the prestress loss in three different tendons at the main (central) pier predicted by the various models in logarithmic time scales. The 100-year prestress force in the top tendon of the last segment, termed "upper cable-16" in Figure 15, is 96% and 92% of the initial force (after the instantaneous losses) when the ACI and EC2 models are used, respectively, but approximately it is 84% when the model B3 is used. Whereas the 100-year prestress force in the lower tendons located along the bottom slab is 98%, 95%, and 87% of the initial value for the ACI, EC2, and B3 model, respectively. Figure 16 shows the normal stresses at the upper and the lower part of different cross section predicted by the various models at two different times: Just after the end of the construction and after 32,000 days (∼85 years). It can be seen clearly that the stresses reduce in time due to the viscosity of the material and that the amount of the stress reduction is very similar to the deflection and prestress losses provided by the different creep models.

Conclusions
The paper presents a general procedure for the modeling of typical creep compliance functions indicated in codes and standards and for its implementation in a three-dimensional finite element code. The procedure is based on the approximation of an aging compliance function through a Dirichlet (or Prony) series, which transforms the classical Volterra integral equation of creep into a rate-type formulation, governed by an aging Kelvin chain whose coefficients are obtained directly from the coefficients of the Dirichlet series approximation of the compliance function. As shown in the manuscript, the calculation of the coefficients can be done in two steps: (1) For the non-aging term of the creep function the coefficients can be efficiently calculated on the basis of the continuous retardation spectrum approach, obtained by the application of the Post-Widder formula for the inversion of Laplace transform (using a third order approximation which reproduces with sufficient accuracy the exact spectrum); (2) after this evaluation all the coefficients of the Dirichlet series must be adjusted step by step in the time integration procedure to take into account the aging term of the creep function. Therefore, according this procedure, the series coefficients have to be calculated at the beginning of the numerical analysis for the determination of the non-aging term of the creep function and then updated at each time increment by multiplying all of them for the same terms to take into account the aging effect.
The particular creep compliance functions considered are those provided by the European Euro-Code 2, the North-American ACI model, the basic creep compliance function and the drying creep compliance function of the B3 model. The reliability of the proposed approach is demonstrated by comparing the numerical solution with the analytical solution for two pre-stress concrete beams, a simple supported beam and a cantilever beam. The 3D numerical model is then been applied to the simulation of the long-term behavior of a real bridge under construction (which was the motivation of this work) because only the three-dimensional analysis can capture the shear lag effects in slabs and webs. Instead of the commonly used simplified beam-type analysis which in general underestimates the deflections and prestress loss, the performed full three-dimensional analysis highlights the effects of the different drying properties (different rates of shrinkage and drying creep) caused by the diverse thicknesses of the upper, lower and lateral slabs in the cross section. In particular, a reliable behavior of the real structure has been obtained with an upward deflection in the first years caused by the differential shrinkage and differential drying creep for all the creep models considered. However, after a few years, a downward deflection sets in, with rate and magnitude which depend on the type of creep model considered.

Acknowledgments:
The authors wish to thank the master student Luisa Aletti for her contribution to some numerical simulations presented in the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Eurocode 2
In Europe the code recommendation for the prediction of creep and shrinkage are provided by the EuroCode2 [42]. The compressive strength of concrete at an age t depends on the type of cement, temperature and curing conditions. For a mean temperature of 20 • C and curing in accordance with ISO 2736/2 [43], the relative compressive strength of concrete at various ages f cm (t) may be estimated from the following equations where f cm (t) is the mean concrete compressive strength at an age of t (expressed in days), f cm is the mean compressive strength after 28 days, β cc is a coefficient which depends on the age t of concrete, s is a coefficient which depends on the type of cement (s = 0.20 for rapid hardening high strength cements RS, 0.25 for normal and rapid hardening cements N and R, and 0.38 for slowly hardening cements SL).
The modulus of elasticity of concrete at an age t = 28 days may be estimated from Equation (A2): where E c (t) is the modulus of elasticity at an age of t days, E ci is the modulus of elasticity at an age of 28 days, β E (t) is a coefficient which depends on the age of concrete with t [days], β cc is a coefficient according to Equation (A1).
The total strain at time t, c (t), of a concrete member uniaxially loaded at time t 0 with a constant stress σ c (t 0 ) may be expressed as follows where ci (t 0 ) is the initial strain at loading, cc (t) is the creep strain at time t > t 0 , cs (t) is the total shrinkage strain, cT (t) is the thermal strain. The prediction model for creep and shrinkage given below predicts the mean behavior of a concrete cross-section. Unless special provisions are given the model is valid for ordinary structural concrete (12 MPa < f ck ≤ 80 MPa) subjected to a compressive stress |σ c | < 0.4 f cm (t 0 ) at an age of loading t 0 and exposed to mean relative humidities in the range of 40 to 100% at mean temperatures from 5 • C to 30 • C. It is accepted that the scope of the model also extends to concrete in tension, though the relations given in the following are directed towards the prediction of creep of concrete subjected to compressive stresses.
Within the range of service stresses |σ c | < 0.4 f cm (t 0 ), creep is assumed to be linearly related to stress. For a constant stress applied at time t 0 this leads to where φ(t, t 0 ) is the creep coefficient and E ci is the modulus of elasticity at the age of 28 days. The stress dependent strain cσ (t, t 0 ) may be expressed as where J(t, t 0 ) is the creep function or creep compliance and E c (t 0 ) is the modulus of elasticity at the time of loading t 0 . The creep coefficient may be calculated from where φ 0 is the notional creep coefficient, β c (t − t 0 ) is the coefficient to describe the development of creep with time after loading, t is the age of concrete (days) at the instant considered t 0 is the age of concrete at loading (days).
The notional creep coefficient may be estimated from where f cm is the mean compressive strength of concrete at the age of 28 days (MPa), f cm0 = 10 MPa, RH is the relative humidity of the ambient environment (%), RH 0 = 100 (%), h = 2A c /u is the notational size of member (mm), where A c is the cross-section and u is the perimeter of the member in contact with the atmosphere h 0 = 100 mm.
The development of creep with time is given by The coefficients α1/2/3 are calculated as function of the strength as The effect of the type of cement on the creep coefficient of concrete may be taken into account by modifying the age t 0 at loading in accordance with Equation (A12) where t 0,T is the age of concrete at loading (days) adjusted according to Equation (A14), α is the power which depends on the type of cement (α = -1 for cement Class S, α = 0 for cement Class N, α = 1 for cement Class R). The effect of elevated or reduced temperatures within the range 0-80 • C on the maturity of concrete may be taken into account by adjusting the concrete age according to the following expression: where t T is the temperature adjusted concrete age which replaces t in the corresponding equations, T(∆t i ) is the temperature in • C during the time period ∆t i , ∆t i is the number of days where a temperature T prevails. For stress levels in the range of 0.4 f cm (t 0 ) < |σ c | < 0.6 f cm (t 0 ) the nonlinearity of creep may be taken into account using Equations (A15) where φ 0,k is the non-linear notional creep coefficient, which replaces φ 0 in Equation (A6), k σ = σ c / f cm (t 0 ) which is the stress-strength ratio, α σ = 1.5. The total shrinkage strain, cs (t) = cd (t) + ca (t), is composed of two components, the drying shrinkage strain, cd (t), and the autogenous shrinkage strain, ca (t). The development of the drying shrinkage strain in time follows from cd (t) = β ds (t, t s )k h cd,0 with β ds (t, t s ) = t − t s t − t s + 0.04 t [days] is the age of the concrete at the moment considered; t s [days] is the age of the concrete at the beginning of drying shrinkage or swelling (normally at the end of curing); f cm [MPa] is the mean compressive strength; f cm0 = 10 MPa; α ds1 is a coefficient which depends on the type of cement: =3 for cement Class S, =4 for cement Class N, =6 for cement Class R; α ds2 is a coefficient which depends on the type of cement: =0, 13 for cement Class S, =0, 12 for cement Class N; =0, 11 for cement Class R. The autogenous shrinkage strain, ca (t), follows from: ca (t) = β as (t) ca (∞) with β as (t) = 1 − exp −0.2t 0.5 and ca (∞) = 2.5( f ck − 10)10 −6 (A18)

Appendix B. ACI Model
In 2008, the American Concrete Institute (ACI) recommended the procedure for the prediction of creep and shrinkage in its code provisions. The most recent version, labeled as 209R-92, was published in 1992 and reapproved in 2008 [44,45]. This procedure is applicable to normal weight and all the light weight concretes (using both moist and steam curing and Types I and III cement) under the standard conditions. Correction factors are applied for conditions which are other than standard. For this model E c (t) = 0.043 ρ 3 f cm28 t a + bt (A19) with ρ = 2500, a = 4, and b = 0.85. The ACI-209R-92 code recommends the following expressions for shrinkage: where t c in days is starting time of drying, f = 26.0exp[1.42 × 10 −2 (V/S)], α = 1, and shu = 780γ sh × 10 −6 with γ sh = γ sh,tc γ sh,RH γ sh,vs γ sh,s γ sh,ψ γ sh,c γ sh,α where γ sh represents the cumulative product of correction factors. In the present work only the following factors are considered. The ambient relative humidity coefficient γ sh,RH is γ sh,RH = where V is the specimen volume in mm 3 and S the specimen surface area in mm 2 . Whereas the other coefficients are not considered (γ sh,tc = γ sh,s = γ sh,ψ = γ sh,c = γ sh,α = 1). The ACI-209R-92 code recommends the following expressions for creep compliance function: where φ u = 2.35, ψ = 0.6, d = 10, γ t 0 , γ λ , and γ h are correction factors for different loading ages t 0 , ambient relative humidity h in decimals, slump of fresh concrete s (in mm), the ratio of fine aggregate to total aggregate by weight ψ expressed as percentage, air content α in percent, and the member size effects through the volume-surface ratio V/S, respectively. where s is the slump of fresh concrete (mm), ψ 1 is the ratio of fine aggregate to total aggregate by weight expressed as percentage, α is the air content in percent. In the numerical analysis reported in this work, the correction factors to allow for the composition of the concrete are not considered, i.e., γ c,s = 0.82, γ c,ψ = 1, and γ sh,α = 1.

Appendix C. Model B3
In 1995, RILEM TC-107-GCS recommended the B3 model which is based on the statistical analysis of creep and shrinkage data in a computerized data bank involving about 15,000 data points and about 100 test series. The model is an improved version of the earlier models namely BP model and BP-KX model [46,47]. The prediction of material parameters of B3 model is restricted to the Portland cement concretes, having a 28-day mean cylinder compressive strength varying from 17 to 70 MPa, w/c ratio 0.30-0.85, a/c ratio 2.5-13.5 and cement content 160-720 kg/m 3 . For this model where D = 2v/s = effective cross-section thickness (v/s = volume to surface ratio of the concrete member), α 1 = 1 for type I cement, =0.85 for type II cement, =1.1 for type III cement, α 2 = 0.75 for steam-curing, =1.2 for for sealed or normal curing in air with initial protection against drying, =1.0 for for curing in water or at 100% relative humidity. k t is a factor given by and k s is the cross-section shape factor as (A33) The average compliance function is expressed as in which q 1 = instantaneous strain due to unit stress, C 0 (t, t 0 ) = compliance function for basic creep, and C d (t, t 0 , t 1 ) = additional compliance function due to simultaneous drying. where t 0 = max(t , t 0 ) if t ≥ t 0 , otherwise C d (t, t , t 0 ) = 0; t 0 is the time at which drying and loading first act simultaneously. The creep coefficient is calculated by the following expression: The model parameters can be estimated from the concrete strength and composition according to the following formulae q 1 = 0.6 × 10 −6 4734 f cm , q 2 = 185.4c 0.5 f −0.9 cm , q 3 = 0.29(w/c) 4 q 2 q 4 = 20.3(a/c) −0.7 , q 5 = 7.57 × 10 5 f −1 cm | sh∞ | −0.6 (A39)