Concrete creep is defined as the time-dependent deformation of a concrete specimen under sustained load. Its magnitude is closely related to the stress applied, time, cement type, mix proportion of concrete, and environmental conditions (ACI Committee 209, 2005). In practice, concrete creep can lead to additional displacements, stress redistribution, and even cracking in concrete structures during their service life (Hubert Rüsch and Hilsdorf, 1983; Bažant et al., 1997). As a result, a prestress loss is often observed in many large-scale prestressed concrete structures, such as long-span bridges and nuclear containments, which could significantly affect their safety and durability.

Under low stress, concrete can be considered an aging viscoelastic material, with concrete creep following the Boltzmann superposition principle. Thus, its strain rate can be expressed as (Bažant and Jirásek, 2018) ε ˙ t = σ ˙ t E t + ∫ 0 t J ˙ t , t 0 d σ t 0 , (1) where t is the age of concrete, t 0 is the time when the load is applied, J ˙ t , t 0 is the first derivative of the compliance function J t , t 0 with respect to t , E t is the instantaneous elastic modulus, and ε ˙ t and σ ˙ t are the strain and stress rates, respectively. Equation (1) can be solved by the finite difference method. However, this method has some limitations. First, the entire stress history is required to obtain the strain increment at the current time step. As a result, in the process of solving creep by the finite element method, the entire stress history needs to be stored for the integration points of each element. For a large-scale problem with many time steps, the evaluation of these history variables is quite time-consuming. Second, the effects of some variable factors, such as temperature, humidity, and concrete cracks, cannot be considered.

It should be noted that the calculation method of J t , t 0 comes from different creep models, reflecting the ratio of strain value to stress of the material at time t (loaded at t 0 ), which is usually related to conditions such as material mix ratio, specimen shape, and environmental factors.

To overcome the difficulties of the finite difference method, Eq. (1) can be transformed into a differential rate-type equation, and an efficient exponential algorithm can be employed to calculate the concrete creep which exhibits a quadratic convergence rate and is unconditionally stable. Zienkiewicz et al. (1968) first applied this method to nonaging viscoelastic materials, and Bažant extended it to aging viscoelastic materials. This method can effectively improve computational efficiency with fixed internal variables. When using the efficient exponential algorithm, it is crucial to select a proper rheological model, such as the Kelvin chain model, to describe the viscoelastic behavior of the material. From a mathematical point of view, the creep compliance function of such a viscoelastic material can be approximated by the Dirichlet series. This can be achieved by curve fitting—the so-called retardation spectrum method.

The curve fitting method is usually based on the least squares method to obtain the coefficients of the Dirichlet series from test data. However, this method lacks actual physical meanings and does not follow the second law of thermodynamics, which sometimes leads to negative coefficients when the test data is not statistically ideal (Schapery, 1962). Many efforts have been made to solve the issue (Baumgaertel and Winter 1989; Elster and Honerkamp, 1991; Kaschta and Schwarzl, 1994; Mead, 1994; Emri and Tschoegl, 1995; Ramkumar et al., 1997; Park and Kim, 2001). Furthermore, for aging viscoelastic materials such as concrete, the coefficients of the Dirichlet series are time-dependent, and the curve fitting method becomes inefficient as the computational process needs additional optimization techniques.

With the continuous retardation spectrum method, the coefficients of the Dirichlet series can be determined by discretizing the continuous retardation spectrum, avoiding the issues encountered in the curve fitting method. Bažant and Xi (1995) studied the continuous retardation spectrum for concrete solidification theory and used the Post–Widder method to approximate the spectrum. In practice, however, a high-order Post–Widder formula is often needed to meet the precision requirement, which significantly increases the analytical complexity. Fortunately, Jirásek and Havlásek (2014) solved this issue using a low-order Post–Widder formula with time adjustment factors of retardation times. A high convergence speed and good accuracy are demonstrated for various concrete creep models. However, the method is heuristically based on empirical analyses, and the determination of time adjustment factors is highly dependent on personal experience and numerical experiments for different creep models.

The purpose of this paper is to develop an improved approach for efficiently approximating the continuous retardation spectra of various concrete creep models. The continuous retardation spectrum is first introduced, then the process of calculating the continuous retardation spectrum by the Post–Widder method and its corresponding shortcomings are analyzed, and an improved approach for solving the continuous retardation spectra based on the Weeks inverse Laplace transform method is proposed. By taking the CEB MC90 creep model as an example, the numerical solution of the continuous retardation spectra solved by the improved approach is analyzed. The proposed approach is then applied to the ACI 209R-92, JSCE, and GL2000 concrete creep models. Finally, the numerical results are compared with the other methods and some conclusions are drawn.

To describe viscoelastic materials, their constitutive properties can be represented by the Kelvin chain model. In the Kelvin chain model, the deformation of a material can be characterized by a number of Kelvin units and an additional spring unit assembled in series (Figure 1). Each Kelvin unit is composed of a spring and a dashpot assembled in parallel. All these units bear the same stress, and the total strain ε is equal to the sum of the deformations of each Kelvin unit and the spring unit. ε = ε 0 + ε 1 + ε 2 + … + ε M . (2)

For the nonaging Kelvin unit μ shown in Figure 2, the total stress is equal to the sum of the elastic stress and the viscous stress. σ t = E μ ε μ t + η μ ε ˙ μ t , (3) where both the spring elastic modulus E μ and the dashpot viscosity η μ do not change with age.

For aging viscoelastic materials such as concrete, the aging Kelvin chain model is needed. For the aging Kelvin unit μ , the stress–strain relationship can be expressed in rate form as σ ˙ t = D μ t ε ˙ μ t + η μ ε ¨ μ t , (4) where the modified age-dependent modulus D μ t is equal to D μ t = E μ t + η ˙ μ t . (5)

Since σ ˙ t is equal to zero for creep tests, the boundary condition of ε μ t 0 = 0 is satisfied and the strain of the aging Kelvin unit μ can be obtained as ε μ t = σ D μ t 0 1 − ⁡ exp − t − t 0 / τ μ , (6) where τ μ = η μ t / D μ t is the retardation time of the Kelvin unit μ and t 0 is the time when the load is applied.

When the aging Kelvin chain model is subjected to a unit stress, the strains of all the Kelvin units ( μ = 0 , 1 , 2 , … , M ) and the spring unit are superimposed; the compliance function is given by (Bažant, 1988) J t , t 0 = 1 E 0 t 0 + ∑ μ = 1 M 1 D μ t 0 1 − ⁡ exp − t − t 0 / τ μ , (7) which can be considered a finite Dirichlet series. In the practical application process, once the compliance function is constructed in the form of Eq. (7), the long-term creep value of materials can be numerically calculated by exponential algorithm. For a given time t 0 , when there are infinite number of Kelvin units in the chain and the retardation time distributes continuously, Eq. (7) can be converted into an integral form (Bažant and Xi, 1995): J t , t 0 = 1 E 0 t 0 + ∫ 0 ∞ 1 τ L τ 1 − ⁡ exp − t − t 0 / τ d τ , (8) where L τ is the continuous retardation spectrum. For a specific creep model, once the continuous retardation spectrum is given, the direct discrete method can be used to construct the compliance function in the form of finite Dirichlet series. For convenience sake, the compliance function J t , t 0 is rewritten as J t , t 0 = 1 E 0 t 0 + φ ξ , (9) where ξ = t − t 0 is the loading duration, and φ ξ is defined as φ ξ = ∫ 0 ∞ 1 τ L τ 1 − ⁡ exp   − ξ / τ d τ . (10)

By using the inverse Laplace transform method, Tschoegl (Nicholas, 1989) approximated L τ from the Post–Widder method as L k τ = − − k τ k k − 1 ! φ k k τ ,   k = 1 , 2 , … , (11) where φ k k τ is the k th order derivative of φ k τ . When k → ∞ , L k τ converges to the exact solution: L τ = lim k → ∞ L k τ . (12)

In this method, a sufficiently smooth function φ ξ and a higher order k are required for an acceptable approximation of the continuous retardation spectrum. This method is straightforward and efficient for some problems with simple compliance functions. In practice, however, the computation of the high-order derivatives is usually quite complicated, and the process of solving them will become complicated. To improve this method, Jirásek and Havlásek (2014) analyzed the difference between the low-order Post–Widder method and the exact solution and proposed an approach for adjusting the retardation time by multiplying correction factor α τ to improve the accuracy of calculating continuous retardation spectrum. Thus, the continuous retardation spectrum is defined as L α τ τ = L k τ . (13)

When the correction factor α τ is applied, the accuracy of the numerical solution of the compliance function is effectively improved, and the approximation order is well controlled within low ranges. In this method, however, the value of the correction factor is determined empirically and is different for different creep models.

In Section 3, in view of the shortcomings of the Post–Widder method, an improved method for solving continuous retardation spectra based on the Weeks inverse Laplace transform is proposed. Here, Eq. (10) is processed in advance.

The differentiation of Eq. (10) with respect to ξ yields φ ′ ξ = ∫ 0 ∞ 1 τ 2 L τ exp   − ξ / τ d τ . (14)

By setting δ = 1 / τ , we have φ ′ ξ = ∫ 0 ∞ L 1 / δ exp   − ξ δ d δ . (15)

It can be seen from Eq. (15) that φ ′ ξ is the Laplace transform of L 1 / δ . Therefore, once the inverse Laplace transform of φ ′ ξ is determined, L 1 / δ can be obtained.

In the method by Weeks (1966), the Laguerre polynomials are used to numerically calculate the inverse Laplace transform. The main advantage is that an explicit solution can be obtained. In applying this method, the following two conditions should be fulfilled: the Laplace space function is a smooth function with bounded exponential growth, and it can be expressed as a Laguerre series. The above two conditions are fulfilled for commonly used creep models, including the CEB MC90, ACI 209R-92, JSCE, and GL2000 models.

For Laplace space function F x , the analytical expression of the time-domain function f t can be obtained by the Weeks method as f t = ⁡ exp σ t ∑ n = 0 ∞ a n ⁡ exp − b t Lg n 2 b t . (16)

In Eq. (16), σ and b are two parameters that fulfill the conditions of b > 0 and σ > σ 0 , where σ 0 is the Laplace convergence abscissa. a n contains the information of the Laplace space function. It could be a scalar, vector, or matrix but does not change with time. With the Maclaurin series, the analytical expression of a n is given by 2 b 1 − ω F σ + 2 b 1 − ω − b = ∑ n = 0 ∞ a n ω n ,   ω < R , (17) where R is the radius of convergence of the Maclaurin series.

L g n x in Eq. (16) is the Laguerre polynomial of degree n defined as L g n x = exp   x n ! d n d x n exp   − x x n . (18)

For practical numerical calculation, Eq. (16) can approximately be expressed as f N t = ⁡ exp σ t ∑ n = 0 N − 1 a n ⁡ exp − b t Lg n 2 b t , (19) where N is an integer. The exact solution of f t can be obtained as N approaches infinity f t = lim N → ∞ f N t . (20)

By using the Bromwich integral and the fast Fourier transform, an approximate expression of a n in Eq. (17) can be obtained as a n ≈ exp − i n π 2 N 2 N ∑ m = − N N − 1 2 b ⁡ exp − i n θ m 1 − ⁡ exp i θ m + 1 / 2 F σ − b exp i θ m + 1 / 2 + 1 exp i θ m + 1 / 2 − 1 , (21) where θ m = m π / N .

In the Laplace transform formula of Eq. (15), the Laplace space function is φ ′ ξ , the complex variable is ξ , the time-domain function is L 1 / δ , and the time-domain variable is δ . From Eq. (19), L 1 / δ can be obtained through the inverse Laplace transform as L N 1 / δ = ⁡ exp σ δ ∑ n = 0 N − 1 a n ⁡ exp − b δ Lg n 2 b δ . (22)

Substitution of τ = 1 / δ into Eq. (22) leads to L N τ = ⁡ exp σ / τ ∑ n = 0 N − 1 a n ⁡ exp − b / τ L g n 2 b / τ , (23) a n ≈ exp − i n π 2 N 2 N ∑ m = − N N − 1 2 b ⁡ exp − i n θ m 1 − ⁡ exp i θ m + 1 / 2 φ ′ σ − b exp i θ m + 1 / 2 + 1 exp i θ m + 1 / 2 − 1 . (24)

In applying the Weeks method, it is obvious that the truncation error can be reduced by using a larger value of N . In addition, a proper choice of σ and b can lead to a higher convergence speed. Therefore, it is very important to choose the values of σ and b reasonably. For Eq. (19), according to the method provided by Weeks (1966), when solving f N t in the range of 0 < t < t max , the value can be simply determined as σ = σ 0 + 1 t max u σ 0 + 1 t max , b = N 2 t max   , (25) u x = 0   if   x < 0 , 1   if   x ≥ 0 . (26)

For Eq. (23), let t max = 1 / τ , and, considering that for all creep models σ 0 = 0 , then σ and b can be taken as σ = τ , b = N τ 2 . (27)

Equation (27) is simple and suitable for all creep models and can directly improve the computational accuracy by increasing N , which is recommended to take a value between 10 and 50. In order to fully improve the computational efficiency, N is taken as 10 in this paper, and, to further determine the values of σ and b for different models, the analysis process is shown as follows.

To determine the parameters σ and b , Weideman (1999) proposed a method by minimizing the theoretical error E E = ⁡ exp σ t ∑ n = N 2 N − 1 a n + ϵ ∑ n = 0 N − 1 a n , (28) where ϵ is the machine round-off error.

Further analysis according to the method of Weideman (Weeks, 1966), taking the CEB MC90 creep model (CEB-FIP, 1993) as an example, the compliance function is J t , t 0 = ρ t 0 + ϕ 0 β c t − t 0 E c m 28 , (29) where E c m 28 is the average elastic modulus of concrete at 28 days, ρ t 0 is the correction factor related to the loading time t 0 , ϕ 0 is the nominal creep coefficient related to the material strength, loading time, and the relative humidity of the environment, and the creep development coefficient β c t − t 0 is equal to β c t − t 0 = t − t 0 β h + t − t 0 0.3 , (30) with β h being a coefficient related to the volume/surface ratio, the relative humidity, and the material strength, generally ranging from 250 to 1,500. By setting φ ξ = β c t − t 0 , the first-order derivative of φ ξ is given by φ ′ ξ = 0.3 β h ξ 0.7 β h + ξ 1.3 . (31)

When β h , N , and τ are given, the error E as a function of σ and b can be obtained by substituting Eqs (24) and (31) into Eq. (28). N should be selected by considering computational accuracy and efficiency. When a larger N is used, the choice of σ and b is more flexible.

By setting β h = 500 and N = 10 , log   E is obtained as a function of σ and b for different retardation times τ (Figure 3). The figure shows that the error E can be significantly reduced by a proper choice of the two parameters. In some zones where σ and b follow a certain relationship, the error E is kept within relatively low levels as the dark part shown in Figure 3.

Combine the three graphs in Figure 3 and represent the horizontal and vertical coordinates in exponential growth. Select several points with the smallest error E within different parameter ranges and use a few crosses to represent them, forming Figure 4.

If the relationship between σ and b can be determined under the condition of minimum error E , the method of parameter selection can be further determined. It was found that the relationship is related to the poles and branch points of the Laplace function (Weideman, 1999). If s 1 = α 1 + β 1 i and s 2 = α 2 + β 2 i are two farthest points to the origin of the Laplace function on the complex plane, we have b 2 − σ 2 + s 2 2 − s 1 2 α 2 − α 1 σ + α 2 s 1 2 − α 1 s 2 2 α 2 − α 1 = 0 . (32)

Equation (31) has two branch points: ξ = − β h and ξ = 0 . By using Eq. (32), the relationship between σ and b for the CEB MC90 model can be obtained as b = β h σ + σ 2 . (33)

By taking β h = 500 , Eq. (33) is plotted as the solid line shown in Figure 4. It is apparent from Figure 4 that the numerical results from error minimization are very close to Eq. (33) for the CEB MC90 model. Therefore, for the CEB MC90 creep model, Eq. (33) can be determined as the relationship between σ and b .

As seen from Eq. (23), the relationship between σ (or b ) and τ is also required to perform the inverse Laplace transform and will be determined empirically through error analysis. Thus, a reference solution L d τ is obtained by a combination of the Durbin inverse Laplace transformation (Durbin, 1974) and an adaptive numerical integration. It should be pointed out that, although the reference solution L d τ has higher accuracy, it is very time-consuming and needs to be calculated by long-term iterative calculation and is thus not suitable for practical applications. The relative error E r is defined as E r = L d τ − L N τ L d τ . (34)

When β h = 500 , N = 10 , under the premise of using Eq. (33) for parameter values, the relative errors of the Weeks inverse Laplace transform for τ = 0.01 , 1 , 100 , 1000 are shown in Figure 5.

Figure 5 shows that, when τ is smaller than 1, the relative error exhibits larger changes from 10⁻⁴ to 1 (Figures 5A,B). When τ is larger than 100, the relative error is smaller than 10⁻¹ and has a downward tendency as σ / τ increases (Figures 5C,D).

Based on the numerical results in Figure 5, an empirical value of σ for the CEB MC90 model is suggested as follows: σ = 0.038 β h τ 2 , for τ < 0.02 β h , σ = 0.038 τ 2 , for 0.02 β h ≤ τ < 100 , σ = τ , for τ ≥ 100 . (35)

For different values of τ , the relative error calculated from Eq. (35) is also obtained as the red cross shown in Figure 5. It is seen from Figure 5 that, when using Eq. (35), the relative error is smaller than 10⁻¹. Particularly when τ is equal to 100 or 1,000, the relative error is smaller than 10⁻³. Therefore, computational accuracy is guaranteed.

For N = 10 and different values of β h , the continuous retardation spectra obtained from different methods are compared in Table 1, where L a τ is calculated from Eqs (33) and (35) while L b τ is calculated from Eq. (27), E r a and E r b are the relative errors between the calculated results obtained by the two calculation methods and the reference solution L d τ , respectively. It is seen from Table 1 that, compared with the Weeks parameter value method—Eq. (27)—the relative errors of the continuous retardation spectra given by Eqs (33) and (35) are smaller in most cases, and the computational accuracy is obviously improved. For the Weeks method, however, the computational accuracy can also be improved by increasing the value of N .

In this section, the Weeks method is applied to the CEB MC90, ACI 209R-92, JSCE, and GL2000 creep models. φ ξ is calculated and analyzed based on the aging Kelvin chain and different methods for the continuous retardation spectrum. N = 10 is used in all cases. The results obtained by a combination of the Durbin inverse Laplace transformation (Durbin, 1974) and an adaptive numerical integration are used as a reference solution.

To further verify the validity of the improved approach, the finite element method is used to compute the long-term creep of concrete by combining the retardation spectrum obtained by the Weeks method with the second-order exponential algorithm (Bažant and Jirásek, 2018). For this purpose, two sets of experimental data, OPC and SF10, were selected from Mazloom et al. (2004). They had different mix proportions, and a pressure of 10 MPa was applied on the OPC and SF10 cylinders on the 28th and 7th days, respectively. Two UMAT user subroutines for material behavior—ACI UMAT and CEB UMAT—of the commercial finite element software ABAQUS were coded. The concrete strain was calculated using the CEB MC90 and ACI 209R-92 creep models. The results are shown in Figure 21, which shows that the finite element results are in good agreement with the experimental results. For the OPC group, the ACI 209R-92 model has higher accuracy, while, for the SF10 group, the CEB MC90 has higher accuracy.

It is noted that one of the purposes of this paper is to ensure that the numerical results agree well with the corresponding expressions of J t , t 0 . As a result, the accuracy of the final numerical results is mainly dependent on whether the analytical expression of J t , t 0 given by design codes is close to the practical situation. Therefore, before the concrete creep is calculated, it is necessary to select a reasonable model according to the practical situation, including the mix proportion, the component shape, loading, and environmental conditions, to ensure the computational accuracy of the numerical results. It should also be noted that, since this paper is mainly concerned with improving on the continuous retardation method, finding a viable, stable strategy to identify the optimal dimension of the approximation in the presence of an error (uncertainty) on data is not addressed. This limitation is expected to be removed in our future work.

Based on the Weeks method, an efficient approximation approach has been developed for the continuous retardation spectra of aging viscoelastic materials. Compared with the existing methods, the approach has several advantages.

(1) It can calculate the continuous retardation spectrum more accurately by only using the first-order derivative of the creep compliance function. The difficulty of calculating the high-order derivatives in the Post–Widder method is avoided.

It should be noted that the proposed approach is only applicable to concrete creep models when the first-order derivative of the compliance function fulfills the requirement of the Weeks method. For concrete creep models with logarithmic compliance functions, such as the fib model (CEB-FIP, 2010), the inversion formula of the Laplace transform has an analytical solution and does not require the Weeks method for the continuous retardation spectrum. In this research, to achieve high computational efficiency and accuracy, the values of σ and b for the CEB MC90, ACI 209R-92, and GL2000 models are determined empirically for N = 10 . As for the JSCE model, a good precision requirement can be achieved when the parameters are taken directly through Eq. (27) for N = 10 . If the proposed approach is applied to other creep models, the parameters can be determined by referring to the analytical process of this paper or directly by increasing the value of N to meet the precision requirements.