## 1 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  
\[\epsilon(t)=\int_{0}^{t}J(t,\tau)\mathrm{d}\sigma(\tau) \tag{1}\]where \(J(t,\tau)\) is the aging creep function and the times \(t\) and \(\tau\) 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).