An improved approach for the continuous retardation spectra of concrete creep and applications

Creep is an important physical property of concrete and can lead to additional displacement, stress redistribution, and even cracking in concrete structures, inducing prestress loss of large-scale prestressed concrete structures. When an exponential algorithm is used to calculate the long-term creep of concrete, it is usually necessary to apply the continuous retardation spectra of the material. In the improved approach proposed here, the continuous retardation spectra can be obtained by the Weeks inverse Laplace transform. The CEB MC90 creep model is taken as an example to analyze the computational process, efficiency, and error of the approach. The improved approach is further applied to the ACI 209R-92, JSCE, and GL2000 concrete creep models. Through comparison with other methods, the advantages of the improved approach are illustrated, and some useful conclusions are drawn.

1 Introduction

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)

$\dot{\epsilon }\left(t\right)=\frac{\dot{\mathit{\sigma }}\left(\mathit{t}\right)}{\mathit{E}\left(\mathit{t}\right)}+{\int }_{\mathbf{0}}^{\mathit{t}}\dot{\mathit{J}}\left(t,t_0\right)\mathrm{d}\sigma \left(t_0\right),(1)$

where $t$ is the age of concrete, $t_0$ is the time when the load is applied, $\dot{J}\left(t,t_0\right)$ is the first derivative of the compliance function $J\left(t,t_0\right)$ with respect to $t$, $E\left(t\right)$ is the instantaneous elastic modulus, and $\dot{\epsilon }\left(t\right)$ and $\dot{\sigma }\left(t\right)$ 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\left(t,t_0\right)$ 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.

2 General solution for the continuous retardation spectrum

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 $\epsilon$ is equal to the sum of the deformations of each Kelvin unit and the spring unit.

$\mathit{\epsilon }={\mathit{\epsilon }}_{\mathbf{0}}+{\mathit{\epsilon }}_{\mathbf{1}}+{\mathit{\epsilon }}_{\mathbf{2}}+\dots +{\mathit{\epsilon }}_{\mathit{M}}.(2)$

FIGURE 1

FIGURE 1. Nonaging Kelvin chain model.

For the nonaging Kelvin unit $\mu$ shown in Figure 2, the total stress is equal to the sum of the elastic stress and the viscous stress.

$\mathit{\sigma }\left(\mathit{t}\right)={\mathit{E}}_{\mathit{\mu }}{\mathit{\epsilon }}_{\mathit{\mu }}\left(\mathit{t}\right)+{\mathit{\eta }}_{\mathit{\mu }}{\dot{\mathit{\epsilon }}}_{\mathit{\mu }}\left(\mathit{t}\right),(3)$

where both the spring elastic modulus $E_{\mu }$ and the dashpot viscosity ${\eta }_{\mu }$ do not change with age.

FIGURE 2

FIGURE 2. Nonaging Kelvin unit.

For aging viscoelastic materials such as concrete, the aging Kelvin chain model is needed. For the aging Kelvin unit $\mu$, the stress–strain relationship can be expressed in rate form as

$\dot{\mathit{\sigma }}\left(\mathit{t}\right)={\mathit{D}}_{\mathit{\mu }}\left(\mathit{t}\right){\dot{\mathit{\epsilon }}}_{\mathit{\mu }}\left(\mathit{t}\right)+{\mathit{\eta }}_{\mathit{\mu }}{\ddot{\mathit{\epsilon }}}_{\mathit{\mu }}\left(\mathit{t}\right),(4)$

where the modified age-dependent modulus $D_{\mu }\left(t\right)$ is equal to

${\mathit{D}}_{\mathit{\mu }}\left(\mathit{t}\right)={\mathit{E}}_{\mathit{\mu }}\left(\mathit{t}\right)+{\dot{\mathit{\eta }}}_{\mathit{\mu }}\left(\mathit{t}\right).(5)$

Since $\dot{\sigma }\left(t\right)$ is equal to zero for creep tests, the boundary condition of ${\epsilon }_{\mu }\left(t_0\right)=0$ is satisfied and the strain of the aging Kelvin unit $\mu$ can be obtained as

${\mathit{\epsilon }}_{\mathit{\mu }}\left(\mathit{t}\right)=\frac{\mathit{\sigma }}{{\mathit{D}}_{\mathit{\mu }}\left({\mathit{t}}_{\mathbf{0}}\right)}\left\{\mathbf{1}-\mathbf{exp}\left[-\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)/{\mathit{\tau }}_{\mathit{\mu }}\right]\right\},(6)$

where ${\tau }_{\mu }={\eta }_{\mu }\left(t\right)/D_{\mu }\left(t\right)$ is the retardation time of the Kelvin unit $\mu$ 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 ($\mu =0,1,2,\dots ,M$) and the spring unit are superimposed; the compliance function is given by (Bažant, 1988)

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}\left({\mathit{t}}_{\mathbf{0}}\right)}+{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathit{M}}}\frac{\mathbf{1}}{{\mathit{D}}_{\mathit{\mu }}\left({\mathit{t}}_{\mathbf{0}}\right)}\left\{\mathbf{1}-\mathbf{exp}\left[-\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)/{\mathit{\tau }}_{\mathit{\mu }}\right]\right\},(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):

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}\left({\mathit{t}}_{\mathbf{0}}\right)}+{\int }_{\mathbf{0}}^{\mathbf{\infty }}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left\{\mathbf{1}-\mathbf{exp}\left[-\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)/\mathit{\tau }\right]\right\}\mathbf{d}\mathit{\tau },(8)$

where $L\left(\tau \right)$ 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\left(t,t_0\right)$ is rewritten as

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}\left({\mathit{t}}_{\mathbf{0}}\right)}+\mathit{\phi }\left(\mathit{\xi }\right),(9)$

where $\xi =t-t_0$ is the loading duration, and $\phi \left(\xi \right)$ is defined as

$\mathit{\phi }\left(\mathit{\xi }\right)={\int }_{\mathbf{0}}^{\mathbf{\infty }}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\right]\mathbf{d}\mathit{\tau }.(10)$

By using the inverse Laplace transform method, Tschoegl (Nicholas, 1989) approximated $L\left(\tau \right)$ from the Post–Widder method as

${\mathit{L}}_{\mathit{k}}\left(\mathit{\tau }\right)=-\frac{{\left(-\mathit{k}\mathit{\tau }\right)}^{\mathit{k}}}{\left(\mathit{k}-\mathbf{1}\right)!}{\mathit{\phi }}^{\left(\mathit{k}\right)}\left(\mathit{k}\mathit{\tau }\right),\hspace{1em}\mathit{k}=\mathbf{1},\mathbf{2},\dots ,(11)$

where ${\phi }^{\left(k\right)}\left(k\tau \right)$ is the $k$th order derivative of $\phi \left(k\tau \right)$. When $k\to \infty$, $L_k\left(\tau \right)$ converges to the exact solution:

$\mathit{L}\left(\mathit{\tau }\right)=\underset{\mathit{k}\to \mathbf{\infty }}{\mathbf{lim}}{\mathit{L}}_{\mathit{k}}\left(\mathit{\tau }\right).(12)$

In this method, a sufficiently smooth function $\phi \left(\xi \right)$ 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 ${\alpha }_{\tau }$ to improve the accuracy of calculating continuous retardation spectrum. Thus, the continuous retardation spectrum is defined as

$\mathit{L}\left({\mathit{\alpha }}_{\mathit{\tau }}\mathit{\tau }\right)={\mathit{L}}_{\mathit{k}}\left(\mathit{\tau }\right).(13)$

When the correction factor ${\alpha }_{\tau }$ 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 $\xi$ yields

${\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)={\int }_{\mathbf{0}}^{\mathbf{\infty }}\frac{\mathbf{1}}{{\mathit{\tau }}^{\mathbf{2}}}\mathit{L}\left(\mathit{\tau }\right)\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\mathbf{d}\mathit{\tau }.(14)$

By setting $\delta =1/\tau$, we have

${\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)={\int }_{\mathbf{0}}^{\mathbf{\infty }}\mathit{L}\left(\mathbf{1}/\mathit{\delta }\right)\mathbf{exp}\left(-\mathit{\xi }\mathit{\delta }\right)\mathbf{d}\mathit{\delta }.(15)$

It can be seen from Eq. (15) that ${\phi }^{\prime }\left(\xi \right)$ is the Laplace transform of $L\left(1/\delta \right)$. Therefore, once the inverse Laplace transform of ${\phi }^{\prime }\left(\xi \right)$ is determined, $L\left(1/\delta \right)$ can be obtained.

3 Weeks method for the continuous retardation spectrum

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\left(x\right)$, the analytical expression of the time-domain function $f\left(t\right)$ can be obtained by the Weeks method as

$\mathit{f}\left(\mathit{t}\right)=\mathbf{exp}\left(\mathit{\sigma }\mathit{t}\right){\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{\infty }}}{\mathit{a}}_{\mathit{n}}\mathbf{exp}\left(-\mathit{b}\mathit{t}\right){\mathbf{Lg}}_{\mathit{n}}\left(\mathbf{2}\mathit{b}\mathit{t}\right).(16)$

In Eq. (16), $\sigma$ and $b$ are two parameters that fulfill the conditions of $b>0$ and $\sigma >{\sigma }_0$, where ${\sigma }_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

$\frac{\mathbf{2}\mathit{b}}{\mathbf{1}-\mathit{\omega }}\mathit{F}\left(\mathit{\sigma }+\frac{\mathbf{2}\mathit{b}}{\mathbf{1}-\mathit{\omega }}-\mathit{b}\right)={\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{\infty }}}{\mathit{a}}_{\mathit{n}}{\mathit{\omega }}^{\mathit{n}},\hspace{1em}\left|\mathit{\omega }\right|<\mathit{R},(17)$

where $R$ is the radius of convergence of the Maclaurin series.

$\mathrm{L}{\mathrm{g}}_n\left(x\right)$ in Eq. (16) is the Laguerre polynomial of degree $n$ defined as

$\mathbf{L}{\mathbf{g}}_{\mathit{n}}\left(\mathit{x}\right)=\frac{\mathbf{exp}\left(\mathit{x}\right)}{\mathit{n}!}\frac{{\mathit{d}}^{\mathit{n}}}{\mathit{d}{\mathit{x}}^{\mathit{n}}}\left[\mathbf{exp}\left(-\mathit{x}\right){\mathit{x}}^{\mathit{n}}\right].(18)$

For practical numerical calculation, Eq. (16) can approximately be expressed as

${\mathit{f}}_{\mathit{N}}\left(\mathit{t}\right)=\mathbf{exp}\left(\mathit{\sigma }\mathit{t}\right){\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{N}-\mathbf{1}}}{\mathit{a}}_{\mathit{n}}\mathbf{exp}\left(-\mathit{b}\mathit{t}\right){\mathbf{Lg}}_{\mathit{n}}\left(\mathbf{2}\mathit{b}\mathit{t}\right),(19)$

where $N$ is an integer. The exact solution of $f\left(t\right)$ can be obtained as $N$ approaches infinity

$\mathit{f}\left(\mathit{t}\right)=\underset{\mathit{N}\to \mathbf{\infty }}{\mathbf{lim}}{\mathit{f}}_{\mathit{N}}\left(\mathit{t}\right).(20)$

By using the Bromwich integral and the fast Fourier transform, an approximate expression of $a_n$ in Eq. (17) can be obtained as

${\mathit{a}}_{\mathit{n}}\approx \frac{\mathbf{exp}\left(-\frac{\mathit{i}\mathit{n}\mathit{\pi }}{\mathbf{2}\mathit{N}}\right)}{\mathbf{2}\mathit{N}}{\displaystyle {\sum }_{\mathit{m}=-\mathit{N}}^{\mathit{N}-\mathbf{1}}}\frac{\mathbf{2}\mathit{b}\mathbf{exp}\left(-\mathit{i}\mathit{n}{\mathit{\theta }}_{\mathit{m}}\right)}{\mathbf{1}-\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)}\mathit{F}\left(\mathit{\sigma }-\mathit{b}\frac{\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)+\mathbf{1}}{\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)-\mathbf{1}}\right),(21)$

where ${\theta }_m=m\pi /N$.

In the Laplace transform formula of Eq. (15), the Laplace space function is ${\phi }^{\prime }\left(\xi \right)$, the complex variable is $\xi$, the time-domain function is $L\left(1/\delta \right)$, and the time-domain variable is $\delta$. From Eq. (19), $L\left(1/\delta \right)$ can be obtained through the inverse Laplace transform as

${\mathit{L}}_{\mathit{N}}\left(\mathbf{1}/\mathit{\delta }\right)=\mathbf{exp}\left(\mathit{\sigma }\mathit{\delta }\right){\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{N}-\mathbf{1}}}{\mathit{a}}_{\mathit{n}}\mathbf{exp}\left(-\mathit{b}\mathit{\delta }\right){\mathbf{Lg}}_{\mathit{n}}\left(\mathbf{2}\mathit{b}\mathit{\delta }\right).(22)$

Substitution of $\tau =1/\delta$ into Eq. (22) leads to

${\mathit{L}}_{\mathit{N}}\left(\mathit{\tau }\right)=\mathbf{exp}\left(\mathit{\sigma }/\mathit{\tau }\right){\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{N}-\mathbf{1}}}{\mathit{a}}_{\mathit{n}}\mathbf{exp}\left(-\mathit{b}/\mathit{\tau }\right)\mathbf{L}{\mathbf{g}}_{\mathit{n}}\left(\mathbf{2}\mathit{b}/\mathit{\tau }\right),(23)$

${\mathit{a}}_{\mathit{n}}\approx \frac{\mathbf{exp}\left(-\frac{\mathit{i}\mathit{n}\mathit{\pi }}{\mathbf{2}\mathit{N}}\right)}{\mathbf{2}\mathit{N}}{\displaystyle {\sum }_{\mathit{m}=-\mathit{N}}^{\mathit{N}-\mathbf{1}}}\frac{\mathbf{2}\mathit{b}\mathbf{exp}\left(-\mathit{i}\mathit{n}{\mathit{\theta }}_{\mathit{m}}\right)}{\mathbf{1}-\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)}{\mathit{\phi }}^{\prime }\left(\mathit{\sigma }-\mathit{b}\frac{\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)+\mathbf{1}}{\mathbf{exp}\left(\mathit{i}{\mathit{\theta }}_{\mathit{m}+\mathbf{1}/\mathbf{2}}\right)-\mathbf{1}}\right).(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 $\sigma$ and $b$ can lead to a higher convergence speed. Therefore, it is very important to choose the values of $\sigma$ and $b$ reasonably. For Eq. (19), according to the method provided by Weeks (1966), when solving $f_N\left(t\right)$ in the range of $0<t<t_{\mathit{max}}$, the value can be simply determined as

$\left\{\begin{array}{l}\mathit{\sigma }=\left({\mathit{\sigma }}_{\mathbf{0}}+\frac{\mathbf{1}}{{\mathit{t}}_{\mathit{max}}}\right)\mathit{u}\left({\mathit{\sigma }}_{\mathbf{0}}+\frac{\mathbf{1}}{{\mathit{t}}_{\mathit{max}}}\right),\\ \mathit{b}=\frac{\mathit{N}}{\mathbf{2}{\mathit{t}}_{\mathit{max}}},\end{array}\right.(25)$

$\mathit{u}\left(\mathit{x}\right)=\left\{\begin{array}{c}\mathbf{0}\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}}\mathit{x}<\mathbf{0},\\ \mathbf{1}\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}}\mathit{x}\ge \mathbf{0}.\end{array}\right.(26)$

For Eq. (23), let $t_{\mathit{max}}=1/\tau$, and, considering that for all creep models ${\sigma }_0=0$, then $\sigma$ and $b$ can be taken as

$\left\{\begin{array}{l}\mathit{\sigma }=\mathit{\tau },\\ \mathit{b}=\frac{\mathit{N}\mathit{\tau }}{\mathbf{2}}.\end{array}\right.(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 $\sigma$ and $b$ for different models, the analysis process is shown as follows.

To determine the parameters $\sigma$ and $b$, Weideman (1999) proposed a method by minimizing the theoretical error $E$

$\mathit{E}=\mathbf{exp}\left(\mathit{\sigma }\mathit{t}\right)\left({\displaystyle {\sum }_{\mathit{n}=\mathit{N}}^{\mathbf{2}\mathit{N}-\mathbf{1}}}\left|{\mathit{a}}_{\mathit{n}}\right|+\mathit{ϵ}{\displaystyle {\sum }_{\mathit{n}=\mathbf{0}}^{\mathit{N}-\mathbf{1}}}\left|{\mathit{a}}_{\mathit{n}}\right|\right),(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

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathit{\rho }\left({\mathit{t}}_{\mathbf{0}}\right)+{\mathit{\varphi }}_{\mathbf{0}}{\mathit{\beta }}_{\mathit{c}}\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}{{\mathit{E}}_{\mathit{c}\mathit{m}\mathbf{28}}},(29)$

where $E_{cm28}$ is the average elastic modulus of concrete at 28 days, $\rho \left(t_0\right)$ is the correction factor related to the loading time $t_0$, ${\varphi }_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 ${\beta }_c\left(t-t_0\right)$ is equal to

${\mathit{\beta }}_{\mathit{c}}\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)={\left[\frac{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}{{\mathit{\beta }}_{\mathit{h}}+\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}\right]}^{\mathbf{0.3}},(30)$

with ${\beta }_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 $\phi \left(\xi \right)={\beta }_c\left(t-t_0\right)$, the first-order derivative of $\phi \left(\xi \right)$ is given by

${\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)=\frac{\mathbf{0.3}{\mathit{\beta }}_{\mathit{h}}}{{\mathit{\xi }}^{\mathbf{0.7}}{\left({\mathit{\beta }}_{\mathit{h}}+\mathit{\xi }\right)}^{\mathbf{1.3}}}.(31)$

When ${\beta }_h$, $N$, and $\tau$ are given, the error $E$ as a function of $\sigma$ 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 $\sigma$ and $b$ is more flexible.

By setting ${\beta }_h=500$ and $N=10$, $\log \left(E\right)$ is obtained as a function of $\sigma$ and $b$ for different retardation times $\tau$ (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 $\sigma$ and $b$ follow a certain relationship, the error $E$ is kept within relatively low levels as the dark part shown in Figure 3.

FIGURE 3

FIGURE 3. Distributions of $\log \left(E\right)$ (CEB MC90 model) for (A) $\tau =10$; (B) $\tau =100$; (C) $\tau =1000$.

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.

FIGURE 4

FIGURE 4. Relationship between $\sigma$ and $b$ (CEB MC90 model).

If the relationship between $\sigma$ 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={\alpha }_1+{\beta }_{1}i$ and $s_2={\alpha }_2+{\beta }_{2}i$ are two farthest points to the origin of the Laplace function on the complex plane, we have

${\mathit{b}}^{\mathbf{2}}-{\mathit{\sigma }}^{\mathbf{2}}+\frac{{\left|{\mathit{s}}_{\mathbf{2}}\right|}^{\mathbf{2}}-{\left|{\mathit{s}}_{\mathbf{1}}\right|}^{\mathbf{2}}}{{\mathit{\alpha }}_{\mathbf{2}}-{\mathit{\alpha }}_{\mathbf{1}}}\mathit{\sigma }+\frac{{{\mathit{\alpha }}_{\mathbf{2}}\left|{\mathit{s}}_{\mathbf{1}}\right|}^{\mathbf{2}}-{{\mathit{\alpha }}_{\mathbf{1}}\left|{\mathit{s}}_{\mathbf{2}}\right|}^{\mathbf{2}}}{{\mathit{\alpha }}_{\mathbf{2}}-{\mathit{\alpha }}_{\mathbf{1}}}=\mathbf{0}.(32)$

Equation (31) has two branch points: $\xi =-{\beta }_h$ and $\xi =0$. By using Eq. (32), the relationship between $\sigma$ and $b$ for the CEB MC90 model can be obtained as

$\mathit{b}=\sqrt{{\mathit{\beta }}_{\mathit{h}}\mathit{\sigma }+{\mathit{\sigma }}^{\mathbf{2}}.}(33)$

By taking ${\beta }_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 $\sigma$ and $b$.

As seen from Eq. (23), the relationship between $\sigma$ (or $b$) and $\tau$ is also required to perform the inverse Laplace transform and will be determined empirically through error analysis. Thus, a reference solution $L_d\left(\tau \right)$ 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\left(\tau \right)$ 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 $Er$ is defined as

$\mathit{E}\mathit{r}=\frac{\left|{\mathit{L}}_{\mathit{d}}\left(\mathit{\tau }\right)-{\mathit{L}}_{\mathit{N}}\left(\mathit{\tau }\right)\right|}{{\mathit{L}}_{\mathit{d}}\left(\mathit{\tau }\right)}.(34)$

When ${\beta }_h=500$, $N=10$, under the premise of using Eq. (33) for parameter values, the relative errors of the Weeks inverse Laplace transform for $\tau =0.01,1,100,1000$ are shown in Figure 5.

FIGURE 5

FIGURE 5. Relative errors of the Weeks method (CEB MC90 model) for (A) $\tau =0.01$; (B) $\tau =1$; (C) $\tau =100$; (D) $\tau =1000$.

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

Based on the numerical results in Figure 5, an empirical value of $\sigma$ for the CEB MC90 model is suggested as follows:

$\left\{\begin{array}{l}\mathit{\sigma }=\left(\frac{\mathbf{0.038}}{{\mathit{\beta }}_{\mathit{h}}}\right){\mathit{\tau }}^{\mathbf{2}},\text{for}\mathit{\tau }<\mathbf{0.02}{\mathit{\beta }}_{\mathit{h}},\\ \mathit{\sigma }=\mathbf{0.038}{\mathit{\tau }}^{\mathbf{2}},\mathbf{\text{for}}\mathbf{0.02}{\mathit{\beta }}_{\mathit{h}}\le \mathit{\tau }<\mathbf{100},\\ \mathit{\sigma }=\mathit{\tau },\mathbf{\text{for}}\mathit{\tau }\ge \mathbf{100}.\end{array}\right.(35)$

For different values of $\tau$, 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 $\tau$ 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 ${\beta }_h$, the continuous retardation spectra obtained from different methods are compared in Table 1, where $L_a\left(\tau \right)$ is calculated from Eqs (33) and (35) while $L_b\left(\tau \right)$ 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\left(\tau \right)$, 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$.

TABLE 1

TABLE 1. Comparison of $L\left(\tau \right)$ calculated from two different methods (CEB MC90 model).

4 Application of the Weeks method to different creep models

In this section, the Weeks method is applied to the CEB MC90, ACI 209R-92, JSCE, and GL2000 creep models. $\phi \left(\xi \right)$ 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.

4.1 CEB MC90 model

As discussed in the previous section, the continuous retardation spectrum of the CEB MC90 model can be obtained by the Weeks method. When ${\beta }_h=500$, the continuous retardation spectra calculated by the Weeks method, the Post–Widder method with different orders, and the reference solution are shown in Figure 6.

FIGURE 6

FIGURE 6. Continuous retardation spectra calculated by different methods (CEB MC90 model).

Once $L\left(\tau \right)$ is known, $\phi \left(\xi \right)$ can be obtained from Eq. (10). However, the integral involved should be approximated by the Dirichlet series for an aging material. Therefore, the discrete retardation times ${\tau }_{\mu }$ are selected based on accuracy and efficiency. If ${\tau }_m$ is assumed to be a geometric series with an initial value of ${{\tau }_0=10}^{-4}$ and a base of 10, that is,

${\mathit{\tau }}_{\mathit{m}}={\mathbf{10}}^{-\mathbf{4}},{\mathbf{10}}^{-\mathbf{3}},\dots ,{\mathbf{10}}^{\mathbf{6}}\mathit{m}=\mathbf{0},\mathbf{1},\mathbf{2},\dots ,\mathbf{10},(36)$

then Eq. (10) is changed to

$\begin{array}{ll}\mathit{\phi }\left(\mathit{\xi }\right)=& {\int }_{\mathit{\tau }=\mathbf{0}}^{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathbf{0}}}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\right]\mathbf{d}\mathit{\tau }\hfill \\ & +{\displaystyle {\sum }_{\mathit{m}=\mathbf{1}}^{\mathbf{10}}}{\int }_{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathit{m}-\mathbf{1}}}^{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathit{m}}}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\right]\mathbf{d}\mathit{\tau }\\ & +{\int }_{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathbf{10}}}^{\mathit{\tau }=+\mathbf{\infty }}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\right]\mathbf{d}\mathit{\tau }.\end{array}(37)$

When $\tau <{\sqrt{10}\tau }_0$, $\left[1-\exp \left(-\xi /\tau \right)\right]$ is close to unity, and the first term in Eq. (37) can be assumed to be $1/E_0$, which is not affected by the duration $\xi$. When $\tau >\sqrt{10}{\tau }_{10}$, $L\left(\tau \right)$ is close to zero, and the third term in Eq. (37) tends to be zero. Thus, Eq. (37) reduces to

$\mathit{\phi }\left(\mathit{\xi }\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}+{\sum }_{\mathit{m}=\mathbf{1}}^{\mathbf{10}}{\int }_{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathit{m}-\mathbf{1}}}^{\mathit{\tau }=\sqrt{\mathbf{10}}{\mathit{\tau }}_{\mathit{m}}}\frac{\mathbf{1}}{\mathit{\tau }}\mathit{L}\left(\mathit{\tau }\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/\mathit{\tau }\right)\right]\mathbf{d}\mathit{\tau }.(38)$

The integral in Eq. (38) can be approximated by the two-point Gaussian quadrature rule:

$\mathit{\phi }\left(\mathit{\xi }\right)\approx \frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}+\frac{\mathbf{ln}\mathbf{10}}{\mathbf{2}}{\displaystyle \sum _{\mathit{m}=\mathbf{1}}^{\mathbf{10}}}\left\{\begin{array}{c}\mathit{L}\left({\mathbf{10}}^{-\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\mathit{m}}\right)\left[\mathbf{1}-\mathbf{exp}\left(-\frac{\mathit{\xi }}{{\mathbf{10}}^{-\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\mathit{m}}}\right)\right]\\ +\mathit{L}\left({\mathbf{10}}^{\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\mathit{m}}\right)\left[\mathbf{1}-\mathbf{exp}\left(-\frac{\mathit{\xi }}{{\mathbf{10}}^{\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\mathit{m}}}\right)\right]\end{array}\right\}.(39)$

The following definitions are used:

$\left\{\begin{array}{l}{\mathit{\tau }}_{\mathit{\mu }}={\mathbf{10}}^{-\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\left(\mathit{\mu }+\mathbf{1}\right)/\mathbf{2}},\mathbf{for\hspace{0.17em}}\mathit{\mu }\mathbf{\hspace{0.17em}is\hspace{0.17em}}\mathbf{o}\mathbf{dd},\\ {\mathit{\tau }}_{\mathit{\mu }}={\mathbf{10}}^{\sqrt{\mathbf{3}}/\mathbf{6}}{\mathit{\tau }}_{\mathit{\mu }/\mathbf{2}},\mathbf{for\hspace{0.17em}}\mathit{\mu }\mathbf{\hspace{0.17em}is\hspace{0.17em}even}.\end{array}\right.(40)$

Equation (39) then becomes

$\mathit{\phi }\left(\mathit{\xi }\right)\approx \frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}+{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathbf{20}}}\frac{\mathbf{ln}\mathbf{10}}{\mathbf{2}}\mathit{L}\left({\mathit{\tau }}_{\mathit{\mu }}\right)\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/{\mathit{\tau }}_{\mathit{\mu }}\right)\right].(41)$

If $1/D\left({\tau }_{\mu }\right)$ is set to be $\ln 10L\left({\tau }_{\mu }\right)/2$, $\phi \left(\xi \right)$ of the aging Kelvin chain model is changed to

$\mathit{\phi }\left(\mathit{\xi }\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}+{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathbf{20}}}\frac{\mathbf{1}}{\mathit{D}\left({\mathit{\tau }}_{\mathit{\mu }}\right)}\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/{\mathit{\tau }}_{\mathit{\mu }}\right)\right].(42)$

Equation (42) is consistent with the form of Eq. (7), where the number of Kelvin units is set to $M=20$. $1/E_0$ can be obtained by numerical integration algorithms (Jirásek and Havlásek, 2014), but this will cause some errors. In this paper, in order to reduce the effect of numerical integration at low retardation times, $1/E_0$ is solved by subtracting the difference between the numerical solution and the exact solution of $\phi \left(\xi \right)$, which is ${\beta }_c\left(t-t_0\right)$ in Eq. (30) for the CEB MC90 model:

$\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}={\mathit{\beta }}_{\mathit{c}}\left({\mathit{\xi }}^{\prime }\right)-{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathbf{20}}}\frac{\mathbf{1}}{\mathit{D}\left({\mathit{\tau }}_{\mathit{\mu }}\right)}\left[\mathbf{1}-\mathbf{exp}\left(-{\mathit{\xi }}^{\prime }/{\mathit{\tau }}_{\mathit{\mu }}\right)\right],(43)$

where ${\xi }^{\prime }$ is the loading duration when ${\beta }_c\left({\xi }^{\prime }\right)$ reaches the maximum value (${\xi }^{\prime }={10}^5$ for the CEB MC90 model). The analysis process is shown in Figure 7.

FIGURE 7

FIGURE 7. Analysis process of $1/E_0$.

By taking ${\beta }_h$ as 500, 800, and 1,200, the values of $\phi \left(\xi \right)$ are obtained by the Weeks and Post–Widder methods with different orders (Figure 8).

FIGURE 8

FIGURE 8. $\phi \left(\xi \right)$ calculated by different methods (CEB MC90 model) for (A) ${\beta }_h=500$; (B) ${\beta }_h=800$; (C) ${\beta }_h=1,200$.

The error is defined as $\left[\phi {\left(\xi \right)-\beta }_c\left(\xi \right)\right]$ for the CEB MC90 model. Thus, the errors of different methods are shown in Figure 9. At the same time, the errors obtained from the method of Jirásek and Havlásek (2014) are also compared in Figure 9.

FIGURE 9

FIGURE 9. Errors of $\phi \left(\xi \right)$ calculated by different methods (CEB MC90 model) for (A) ${\beta }_h=500$; (B) ${\beta }_h=800$; (C) ${\beta }_h=1,200$.

Equation (42) shows that $\phi \left(\xi \right)$ is related to ${\tau }_{\mu }$. For the two-point Gaussian quadrature rule, $L\left({\tau }_{\mu }\right)$ corresponding to different ${\tau }_{\mu }$ is able to describe the retardation spectrum curve as accurately as possible. If $L\left({\tau }_{\mu }\right)$ at two adjacent points ${\tau }_{\mu }$ are connected by line segments, the approximation of $L\left({\tau }_{\mu }\right)$ to the retardation spectrum curve can be observed more intuitively. As shown in Figure 10, the retardation spectrum curve can be described well when ${\tau }_{\mu }$ takes values from Eq. (40).

FIGURE 10

FIGURE 10. Approximation of $L\left({\tau }_{\mu }\right)$ to the retardation spectrum (CEB MC90 model and ${\beta }_h=500$).

In Eq. (42), the number of Kelvin units is taken as 20. If higher computational accuracy is required, the number of Kelvin units can be increased: $M=20k,(k=1,2,3,\dots$). Thus, Eq. (42) becomes

$\mathit{\phi }\left(\mathit{\xi }\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}+{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathit{M}}}\frac{\mathbf{1}}{\mathit{D}\left({\mathit{\tau }}_{\mathit{\mu }}\right)}\left[\mathbf{1}-\mathbf{exp}\left(-\mathit{\xi }/{\mathit{\tau }}_{\mathit{\mu }}\right)\right],(44)$

where

${\mathit{\tau }}_{\mathit{m},\mathit{k}}={\mathbf{10}}^{-\mathbf{4}},{\mathbf{10}}^{-\mathbf{4}+\mathbf{1}/\mathit{k}},{\mathbf{10}}^{-\mathbf{4}+\mathbf{2}/\mathit{k}},\dots ,{\mathbf{10}}^{\mathbf{6}}\mathit{m}=\mathbf{0},\mathbf{1},\mathbf{2},\dots ,\mathbf{10}\mathit{k},(45)$

$\left\{\begin{array}{l}{\mathit{\tau }}_{\mathit{\mu }}={\mathbf{10}}^{-\sqrt{\mathbf{3}}/\left(\mathbf{6}\mathit{k}\right)}{\mathit{\tau }}_{\left(\mathit{\mu }+\mathbf{1}\right)/\mathbf{2},\mathit{k}},\mathbf{for}\mathit{\mu }\mathbf{\hspace{0.17em}is\hspace{0.17em}odd},\\ {\mathit{\tau }}_{\mathit{\mu }}={\mathbf{10}}^{\sqrt{\mathbf{3}}/\left(\mathbf{6}\mathit{k}\right)}{\mathit{\tau }}_{\mathit{\mu }/\mathbf{2},\mathit{k}},\mathbf{for}\mathit{\mu }\mathbf{\hspace{0.17em}is\hspace{0.17em}even},\end{array}\right.(46)$

$\frac{\mathbf{1}}{{\mathit{E}}_{\mathbf{0}}}={\mathit{\beta }}_{\mathit{c}}\left({\mathit{\xi }}^{\prime }\right)-{\displaystyle {\sum }_{\mathit{\mu }=\mathbf{1}}^{\mathit{M}}}\frac{\mathbf{1}}{\mathit{D}\left({\mathit{\tau }}_{\mathit{\mu }}\right)}\left[\mathbf{1}-\mathbf{exp}\left(-{\mathit{\xi }}^{\prime }/{\mathit{\tau }}_{\mathit{\mu }}\right)\right],(47)$

$\frac{\mathbf{1}}{\mathit{D}\left(\mathit{\tau }\right)}=\frac{\mathbf{ln}\left({\mathbf{10}}^{\frac{\mathbf{1}}{\mathit{k}}}\right)\mathit{L}\left(\mathit{\tau }\right)}{\mathbf{2}}.(48)$

It should be noted that, for most concrete creep models (CEB MC90, ACI 209R-92, and GL2000), the retardation spectra are relatively smooth, and taking $k=1$ in Eq. (45) usually satisfies the precision requirement. However, if the retardation spectrum changes sharply in terms of retardation time, it is necessary to increase the value of $k$, which will be discussed in Section 4.3.

4.2 ACI 209R-92 model

The creep compliance function of the ACI 209R-92 model (ACI committee 209, 2008) is

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathbf{1}+\mathit{\varphi }\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)}{\mathit{E}\left({\mathit{t}}_{\mathbf{0}}\right)},(49)$

where $\varphi \left(t,t_0\right)$ is defined as

$\mathit{\varphi }\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}^{\mathit{\psi }}}{\mathit{d}+{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}^{\mathit{\psi }}}{\mathit{\varphi }}_{\mathit{\mu }}.(50)$

In Eq. (50), ${\varphi }_{\mu }$ is the ultimate creep coefficient related to the curing conditions, slump, and ambient humidity of concrete structures and $\psi$ and $d$ are two parameters—usually taken as 0.6 and 10 according to the recommendations in ACI 209R-92, respectively.

To simplify the analysis process, $\phi \left(\xi \right)$ is defined as

$\mathit{\phi }\left(\mathit{\xi }\right)=\frac{{\mathit{\xi }}^{\mathit{\psi }}}{\mathit{d}+{\mathit{\xi }}^{\mathit{\psi }}}.(51)$

Derivation of Eq. (51) with respect to $\xi$ gives

${\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)=\frac{{\mathit{d}\mathit{\psi }\mathit{\xi }}^{\mathit{\psi }-\mathbf{1}}}{{\left(\mathit{d}+{\mathit{\xi }}^{\mathit{\psi }}\right)}^{\mathbf{2}}}.(52)$

Since Eq. (52) has a branch point of $\xi =0$, Eq. (32) cannot be directly used. For the ACI 209R-92 model, when $N=10$, the computational accuracy of the creep compliance function obtained from Eq. (27) is not very satisfactory (Figure 13). With reference to the analysis of the CEB MC90 model, Eq. (28) is used to calculate $E$ and analyze the relationship between $\sigma$ and $b$. Thus $b$ is given by

$\mathit{b}=\sqrt{\mathbf{46}\mathit{\sigma }+{\mathit{\sigma }}^{\mathbf{2}}.}(53)$

When $N=10$, the empirical expression of $\sigma$ is expressed as

$\left\{\begin{array}{l}\mathit{\sigma }=\mathbf{0.05}{\mathit{\tau }}^{\mathbf{2}},\mathbf{\text{for}}\mathit{\tau }<\mathbf{20},\\ \mathit{\sigma }=\mathit{\tau },\mathbf{\text{for}}\mathit{\tau }\ge \mathbf{20}.\end{array}\right.(54)$

The continuous retardation spectra calculated by the Weeks method and the Post–Widder method with different orders are shown in Figure 11.

FIGURE 11

FIGURE 11. Continuous retardation spectra calculated by different methods (ACI 209R-92 model).

From Eq. (42), the values of $\phi \left(\xi \right)$ of the ACI 209R-92 model for different durations $\xi$ are obtained as shown in Figure 12; the errors are shown in Figure 13. For comparison, the errors of $\phi \left(\xi \right)$ calculated from Eq. (27) are also shown in Figure 13 for $N=10$ and $N=50$. Figure 13 shows that, when Eq. (53) and (54) are adopted, the computational accuracy is obviously improved under the same $N$ value.

FIGURE 12

FIGURE 12. $\phi \left(\xi \right)$ calculated by different methods (ACI 209R-92 model).

FIGURE 13

FIGURE 13. Errors of $\phi \left(\xi \right)$ calculated by different methods (ACI 209R-92 model).

4.3 JSCE model

The creep compliance function of the JSCE model (Uomoto et al., 2008) is defined as

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\left\{\mathbf{1}-\mathbf{exp}\left[-\mathbf{0.09}{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}^{\mathbf{0.6}}\right]\right\}{\mathit{\epsilon }}_{\mathit{c}\mathit{r}}^{\prime },(55)$

where ${\epsilon }_{cr}^{\prime }$ is the ultimate creep strain under unit stress and related to the ambient relative humidity, temperature, and the volume/surface ratio. $\phi \left(\xi \right)$ is defined as

$\mathit{\phi }\left(\mathit{\xi }\right)=\mathbf{1}-\mathbf{exp}\left(-\mathbf{0.09}{\mathit{\xi }}^{\mathbf{0.6}}\right).(56)$

Derivation of $\phi \left(\xi \right)$ with respect to $\xi$ gives

${\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)=\mathbf{0.054}{\mathit{\xi }}^{-\mathbf{0.4}}\mathbf{exp}\left(-\mathbf{0.09}{\mathit{\xi }}^{\mathbf{0.6}}\right).(57)$

Since Eq. (57) has only a branch point $\xi =0$, Eq. (32) cannot be directly used to evaluate the values of $\sigma$ and $b$. Instead, Eq. (27) is adopted to determine the values of the two parameters. For $N=10$, the continuous spectra calculated by the Weeks and Post–Widder methods with different orders are shown in Figure 14.

FIGURE 14

FIGURE 14. Continuous retardation spectra calculated by different methods (JSCE model).

It is seen from Figure 14 that, compared with the CEB MC90 and ACI 209R-92 models, the continuous retardation spectra calculated by the Weeks method are steeper at the peak. The continuous retardation spectra of the CEB and ACI models show significant changes in the range from $\tau ={10}^{-2}$ to $\tau ={10}^5$, while that of the JSCE model mainly shows changes in the range from $\tau ={10}^{-1}$ to $\tau ={10}^3$. If $k=1$ is set in Eq. (45), $L\left({\tau }_{\mu }\right)$ cannot effectively describe the retardation spectrum curve. As shown in Figure 15A, the curve around $\tau ={10}^2$ is not well reproduced. Therefore, it is necessary to increase the value of $k$. By setting the number of Kelvin units to $M=40$— $k=2$ in Eq. (45)—it is seen from Figure 15B that the computational accuracy is greatly improved. Therefore, for the JSCE model, it is necessary to adopt this retardation time value method. The numerical results and the errors for different methods are shown in Figures 16, 17, respectively.

FIGURE 15

FIGURE 15. Approximation of $L\left({\tau }_{\mu }\right)$ to the retardation spectrum (JSCE model) for (A) $k=1$; (B) $k=2$.

FIGURE 16

FIGURE 16. $\phi \left(\xi \right)$ calculated by different methods (JSCE model).

FIGURE 17

FIGURE 17. Errors of $\phi \left(\xi \right)$ calculated by different methods (JSCE model).

4.4 GL2000 model

The creep compliance function of the GL2000 model (Gardner and Lockman, 2002) is defined as

$\mathit{J}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)=\frac{\mathbf{1}}{{\mathit{E}}_{\mathit{c}\mathit{m}\mathit{t}\mathbf{0}}}+\frac{{\mathit{\varphi }}_{\mathbf{28}}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)}{{\mathit{E}}_{\mathit{c}\mathit{m}\mathbf{28}}},(58)$

where $E_{cmt0}$ is the elastic modulus at loading time and ${\varphi }_{28}\left(t,t_0\right)$ is the creep coefficient at 28 days. The expression of ${\varphi }_{28}\left(t,t_0\right)$ is

${\mathit{\varphi }}_{\mathbf{28}}\left(\mathit{t},{\mathit{t}}_{\mathbf{0}}\right)={\mathit{\varphi }}_{\mathit{c}}\left\{\begin{array}{c}\frac{\mathbf{2}{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}^{\mathbf{0.3}}}{{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}^{\mathbf{0.3}}+\mathbf{14}}+{\left(\frac{\mathbf{7}}{{\mathit{t}}_{\mathbf{0}}}\right)}^{\mathbf{0.5}}{\left[\frac{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)+\mathbf{7}}\right]}^{\mathbf{0.5}}\\ +\mathbf{2.5}\left(\mathbf{1}-\mathbf{1.086}{\mathit{h}}^{\mathbf{2}}\right){\left[\frac{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)}{\left(\mathit{t}-{\mathit{t}}_{\mathbf{0}}\right)+\mathbf{0.12}{\left(\mathit{V}/\mathit{S}\right)}^{\mathbf{2}}}\right]}^{\mathbf{0.5}}\end{array}\right\}.(59)$

In Eq. (59), ${\varphi }_c$ is the correction term for the drying effect before loading, $h$ is the relative humidity, and $V/S$ is the volume/surface ratio of specimens. To simplify the calculation, $\phi \left(\xi \right)$ is expressed as

$\mathit{\phi }\left(\mathit{\xi }\right)=\mathbf{2}\frac{{\mathit{\xi }}^{\mathbf{0.3}}}{{\mathit{\xi }}^{\mathbf{0.3}}+\mathbf{14}}+{\left(\frac{\mathbf{7}}{{\mathit{t}}_{\mathbf{0}}}\right)}^{\mathbf{0.5}}{\left(\frac{\mathit{\xi }}{\mathit{\xi }+\mathbf{7}}\right)}^{\mathbf{0.5}}+\mathbf{2.5}\left(\mathbf{1}-\mathbf{1.086}{\mathit{h}}^{\mathbf{2}}\right){\left[\frac{\mathit{\xi }}{\mathit{\xi }+\mathbf{0.12}{\left(\mathit{V}/\mathit{S}\right)}^{\mathbf{2}}}\right]}^{\mathbf{0.5}}.(60)$

If ${\phi }_1\left(\xi \right)$ and ${\phi }_2\left(\xi ,\gamma \right)$ are defined as

${\mathit{\phi }}_{\mathbf{1}}\left(\mathit{\xi }\right)=\frac{{\mathit{\xi }}^{\mathbf{0.3}}}{{\mathit{\xi }}^{\mathbf{0.3}}+\mathbf{14}},\hspace{1em}{\mathit{\phi }}_{\mathbf{2}}\left(\mathit{\xi },\mathrm{\gamma }\right)={\left(\frac{\mathit{\xi }}{\mathit{\xi }+\mathrm{\gamma }}\right)}^{\mathbf{0.5}},(61)$

the derivation of $\phi \left(\xi \right)$ with respect to $\xi$ yields

$\begin{array}{ll}\hfill {\mathit{\phi }}^{\prime }\left(\mathit{\xi }\right)& =\mathbf{2}{\mathit{\phi }}_{\mathbf{1}}^{\prime }\left(\mathit{\xi }\right)+{\left(\frac{\mathbf{7}}{{\mathit{t}}_{\mathbf{0}}}\right)}^{\mathbf{0.5}}{\mathit{\phi }}_{\mathbf{2}}^{\prime }\left(\mathit{\xi },\mathbf{7}\right)+\mathbf{2.5}\left(\mathbf{1}-\mathbf{1.086}{\mathit{h}}^{\mathbf{2}}\right){\mathit{\phi }}_{\mathbf{2}}^{\prime }\hfill \\ \hfill & \left(\mathit{\xi },\mathbf{0.12}{\left(\mathit{V}/\mathit{S}\right)}^{\mathbf{2}}\right).\hfill \end{array}(62)$

Because the function form of ${\phi }_2$ is similar to Eq. (30), and the function form of ${\phi }_1$ is similar to Eq. (51), the parameter value forms of the CEB MC90 and ACI 209R-92 models can be referred to respectively. For $N=10$, the values of $\sigma$ and $b$ recommended for ${\phi }_1^{\prime }\left(\xi \right)$ and ${\phi }_2^{\prime }\left(\xi ,\gamma \right)$ are

${\mathit{b}}_{\mathbf{1}}={\mathit{\sigma }}_{\mathbf{1}},{\mathit{\sigma }}_{\mathbf{1}}=\mathbf{0.6}\mathit{\tau },(63)$

${\mathit{b}}_{\mathbf{2}}=\sqrt{{\mathit{\sigma }}_{\mathbf{2}}^{\mathbf{2}}+\mathrm{\gamma }{\mathit{\sigma }}_{\mathbf{2}}},\left\{\begin{array}{l}{\mathit{\sigma }}_{\mathbf{2}}=\left(\mathbf{15}/\mathrm{\gamma }\right){\mathit{\tau }}^{\mathbf{2}},\mathbf{\text{for}}\mathit{\tau }<\mathbf{0.04}\mathrm{\gamma },\\ {\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.6}\mathit{\tau },\mathbf{\text{for}}\mathbf{0.04}\mathrm{\gamma }\le \mathit{\tau }<\mathbf{100},\\ {\mathit{\sigma }}_{\mathbf{2}}=\mathit{\tau },\mathbf{\text{for}}\mathit{\tau }\ge \mathbf{100}.\end{array}\right.(64)$

Two cases are considered, and the parameters are listed in Table 2. The continuous retardation spectra calculated by the Weeks and Post–Widder methods for different orders are shown in Figure 18.

TABLE 2

TABLE 2. Parameters for the two cases.

FIGURE 18

FIGURE 18. Continuous retardation spectra calculated by different methods (GL2000 model) for (A) case 1; (B) case 2.

It is apparent from Figure 18 that the continuous retardation spectra begin to grow significantly from $\tau ={10}^{-4}$ and decrease gradually at a long retardation time. Therefore, more Kelvin units are needed to cover the range of the continuous retardation spectra. Consequently, ${\tau }_0$ is set to ${10}^{-5}$, the maximum value of $m$ in Eq. (36) is taken as 18, ${\tau }_{18}={10}^{13}$ for this model, and the number of Kelvin units in Eq. (42) changes from 20 to 36. The values of $\phi \left(\xi \right)$ of the GL2000 model for different duration $\xi$ are obtained as per Figure 19, and the errors are shown in Figure 20. Since the GL2000 model was not discussed by Jirásek and Havlásek (2014), the error obtained from their method is not compared in Figure 20.

FIGURE 19

FIGURE 19. $\phi \left(\xi \right)$ calculated by different methods (GL2000 model) for (A) case 1; (B) case 2.

FIGURE 20

FIGURE 20. Errors of $\phi \left(\xi \right)$ calculated by different methods (GL2000 model) for (A) case 1; (B) case 2.

4.5 Discussion

For all the concrete creep models discussed in this section, the proposed approach based on the Weeks method improves the computational accuracy and efficiency of the continuous retardation spectra and the creep compliance. It should be noted that, for the case of long duration, the errors of $\phi \left(\xi \right)$ decrease with the increase of $\xi$ and approach zero. There are two main reasons for this phenomenon. First, the creep compliances used have negligible variation for a long duration and have a clear upper limit. Thus, the approximation through the Dirichlet series will maintain the same characteristics as long as the retardation time range selected is wide enough. Second, by calculating $1/E_0$ based on the most accurate compliance value—Eq. (43)—the effect of numerical integration for a short duration on the final results is reduced.

For the creep models considered, when $\xi \ge 1$, the absolute error of $\phi \left(\xi \right)$ given by the improved approach (using the Weeks method) is smallest and maintains a lower level (usually smaller than 0.01). For the ACI 209R-92 and JSCE models, the results from the proposed approach achieve higher accuracy than the Post–Widder method. When calculated by the improved approach, $\phi \left(\xi \right)$ of the ACI 209R-92 model is closest to the reference solution, while $\phi \left(\xi \right)$ of the JSCE model exhibits the most significant improvement on the Post–Widder method.

The errors of $\phi \left(\xi \right)$ obtained from the improved approach and the Jirásek–Havlásek method are both controlled at a relatively low level. For the CEB MC90 model, the results obtained from Jirásek–Havlásek are more accurate for low retardation times, while those obtained from Weeks are more accurate for high retardation times. As shown in Figure 9, the errors corresponding to the Weeks method are smaller compared to the Jirásek–Havlásek method when $\xi >10$. For the ACI 209R-9 and JSCE models, the results obtained from Weeks are more accurate and the errors of $\phi \left(\xi \right)$ are smaller.

For the other creep models which are not discussed in this paper, the proposed approach is still applicable when they fulfill the requirement of the Weeks inverse transform. In general, a large value of $N$ can be used to achieve the required computational accuracy and a proper choice of $\sigma$ and $b$ can significantly improve computational efficiency. In general, Eq. (27) can be used to determine the parameters, and $10\le N\le 50$ is recommended. The larger the value of $N$, the more accurate the calculation results, but computational efficiency will be reduced. If ${\phi }^{\prime }\left(\xi \right)$ has two or more poles and branch points, Eq. (32) can also be used to determine the parameters, and precision requirements are usually satisfied when $N=10$.

5 Comparison with experimental results

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.

FIGURE 21

FIGURE 21. Comparison between experimental results and numerical results for (A) OPC; (B) SF10.

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\left(t,t_0\right)$. As a result, the accuracy of the final numerical results is mainly dependent on whether the analytical expression of $J\left(t,t_0\right)$ 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.

6 Conclusion

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.

(2) Unlike the method proposed by Jirásek and Havlásek (2014), in which the correction factor is empirically determined for each concrete creep model at a given derivative order, the proposed approach is based on a solid theoretical foundation and can be conveniently applied to various concrete creep models.

(3) Better computational accuracy can be achieved for a long loading duration. As illustrated by different concrete creep models, the error of the creep compliance function obtained by the proposed approach is controlled within 0.02 for a loading duration of ${10}^{-1}<\xi <{10}^6$. Therefore, the approach is applicable to long-term creep analyses of aging viscoelastic materials, such as concrete.

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 $\sigma$ 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.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

XZ: Conceptualization, Methodology, Writing–review and editing. LB: Formal Analysis, Validation, Writing–original draft. HR: Data curation, Software, Writing–original draft. XF: Funding acquisition, Project administration, Supervision, Writing–original draft. JZ: Funding acquisition, Resources, Writing–review and editing. YG: Investigation, Visualization, Writing–original draft.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by the National Science Foundation for Excellent Young Scholars (Grant No. 52222808), the Zhejiang Provincial Natural Science Foundation (Grant No. LY20E080027), and the National Natural Science Foundation of the People’s Republic of China (Grant Nos 52008413, 52078509, 52178255, and 52278279).

Conflict of interest

Authors HR, XF, and YG were employed by Metallurgical Group Corporation of China.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor J-GD declared a past co-authorship with the author JZ.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmats.2024.1340883/full#supplementary-material

References

ACI Committee 209 (2005). Report on factors affecting shrinkage and creep of hardened concrete. Farmington Hills: American Concrete Institute.

Google Scholar

ACI Committee 209 (2008). Guide for modeling and calculating shrinkage and creep in hardened concrete. Farmington Hills: American Concrete Institute.

Google Scholar

Baumgaertel, M., and Winter, H. H. (1989). Determination of discrete relaxation and retardation time spectra from dynamic mechanical data. Rheol. Acta 28, 511–519. doi:10.1007/BF01332922

CrossRef Full Text | Google Scholar

Bažant, Z. P. (1988). Mathematical modeling of creep and shrinkage of concrete. New York: John Wiley and Sons.

Google Scholar

Bažant, Z. P., Hauggaard, A. B., Baweja, S., and Ulm, F. J. (1997). Microprestress-solidification theory for concrete creep. I: aging and drying effects. J. Eng. Mech. 123, 1188–1194. doi:10.1061/(asce)0733-9399(1997)123:11(1188)(1997)123:11(1188)

CrossRef Full Text | Google Scholar

Bažant, Z. P., and Jirásek, M. (2018). Creep and hygrothermal effects in concrete structures. Netherland: Springer Dordrecht.

Google Scholar

Bažant, Z. P., and Xi, Y. (1995). Continuous retardation spectrum for solidification theory of concrete creep. J. Eng. Mech. 121, 281–288. doi:10.1061/(asce)0733-9399(1995)121:2(281)(1995)121:2(281)

CrossRef Full Text | Google Scholar

CEB-FIP (1993). Design of concrete structures, CEB-FIP model-code 1990. London: British Standard Institution.

Google Scholar

CEB-FIP (2010). Fib Model Code for concrete structures 2010. London: British Standard Institution.

Google Scholar

Durbin, F. (1974). Numerical inversion of laplace transforms: an efficient improvement to dubner and abate's method. Comput. J. 17, 371–376. doi:10.1093/comjnl/17.4.371

CrossRef Full Text | Google Scholar

Elster, C., and Honerkamp, J. (1991). Modified maximum entropy method and its application to creep data. Macromolecules 24, 310–314. doi:10.1021/ma00001a047

CrossRef Full Text | Google Scholar

Emri, I., and Tschoegl, N. W. (1995). Determination of mechanical spectra from experimental responses. Int. J. Solids Struct. 32, 817–826. doi:10.1016/0020-7683(94)00162-P

CrossRef Full Text | Google Scholar

Gardner, N. J., and Lockman, M. J. (2002). Design provisions for drying shrinkage and creep of normal-strength concrete. ACI Mater. J. 98, 159–167. doi:10.14359/10199

CrossRef Full Text | Google Scholar

Hubert Rüsch, D. J., and Hilsdorf, H. K. (1983). Creep and shrinkage. New York: Springer.

Google Scholar

Jirásek, M., and Havlásek, P. (2014). Accurate approximations of concrete creep compliance functions based on continuous retardation spectra. Comput. Struct. 135, 155–168. doi:10.1016/j.compstruc.2014.01.024

CrossRef Full Text | Google Scholar

Kaschta, J., and Schwarzl, R. R. (1994). Calculation of discrete retardation spectra from creep data-Ⅰ. Method. Rheol. Acta 33, 517–529. doi:10.1007/BF00366336

CrossRef Full Text | Google Scholar

Mazloom, M., Ramezanianpour, A. A., and Brooks, J. J. (2004). Effect of silica fume on mechanical properties of high-strength concrete. Cem. Concr. Compos. 26, 347–357. doi:10.1016/s0958-9465(03)00017-9

CrossRef Full Text | Google Scholar

Mead, D. W. (1994). Numerical interconversion of linear viscoelastic material functions. J. Rheology 38, 1769–1795. doi:10.1122/1.550526

CrossRef Full Text | Google Scholar

Nicholas, W. T. (1989). The phenomenological theory of linear viscoelastic behavior. Heidelberg: Springer.

Google Scholar

Park, S. W., and Kim, Y. R. (2001). Fitting prony-series viscoelastic models with power-law presmoothing. J. Mater. Civ. Eng. 13, 26–32. doi:10.1061/(asce)0899-1561(2001)13:1(26)(2001)13:1(26)

CrossRef Full Text | Google Scholar

Ramkumar, D. H. S., Caruthers, J. M., Mavridis, H., and Shroff, R. (1997). Computation of the linear viscoelastic relaxation spectrum from experimental data. J. Appl. Polym. Sci. 64, 2177–2189. doi:10.1002/(SICI)1097-4628(19970613)64:11<2177::AID-APP14>3.0.CO;2-1

CrossRef Full Text | Google Scholar

Schapery, R. (1962). "A simple collocation method for fitting viscoelastic models to experimental data", in: Graduate Aeronautical Laboratory. (California: California Institute of Technology).

Google Scholar

Uomoto, T., Ishibashi, T., Nobuta, Y., Satoh, T., Kawano, H., Takewaka, K., et al. (2008). Standard specifications for concrete structures-2007 by Japan society of Civil engineers. Concr. J. 46, 3–14. doi:10.3151/coj1975.46.7_3

CrossRef Full Text | Google Scholar

Weeks, W. T. (1966). Numerical inversion of laplace transforms using Laguerre functions. J. ACM 13, 419–429. doi:10.1145/321341.321351

CrossRef Full Text | Google Scholar

Weideman, J. a.C. (1999). Algorithms for parameter selection in the weeks method for inverting the laplace transform. Siam J. Sci. Comput. 21, 111–128. doi:10.1137/s1064827596312432

CrossRef Full Text | Google Scholar

Zienkiewicz, O. C., Watson, M., and King, I. P. (1968). A numerical method of visco-elastic stress analysis. Int. J. Mech. Sci. 10, 807–827. doi:10.1016/0020-7403(68)90022-2

CrossRef Full Text | Google Scholar