Gradient Nonlocal Enhanced Microprestress-Solidification Theory and Its Finite Element Implementation

Time-dependent responses of cracked concrete structures are complex, due to the intertwined effects between creep, shrinkage, and cracking. There still lacks an effective numerical model to accurately predict their nonlinear long-term deflections. To this end, a computational framework is constructed, of which the aforementioned intertwined effects are properly treated. The model inherits merits of gradient-enhanced damage (GED) model and microprestress-solidification (MPS) theory. By incorporating higher order deformation gradient, the proposed GED-MPS model circumvents damage localization and mesh-sensitive problems encountered in classical continuum damage theory. Moreover, the model reflects creep and shrinkage of concrete with respect to underlying moisture transport and heat transfer. Residing on the Kelvin chain model, rate-type creep formulation works fully compatible with the gradient nonlocal damage model. 1-D illustration of the model reveals that the model could regularize mesh-sensitivity of nonlinear concrete creep affected by cracking. Furthermore, the model depicts long-term deflections and cracking evolutions of simply-supported reinforced concrete beams in an agreed manner. It is noteworthy that the gradient nonlocal enhanced microprestress-solidification theory is implemented in the general finite element software Abaqus/Standard with the implicit solver, which renders the model suitable for large-scale creep-sensitive structures.


INTRODUCTION
Concrete creep is the deformation phenomenon developing over time under the load action (Bažant and Li 2008). The current mechanical models are suitable for describing serviceability limit states (Schlappal et al., 2020). Concrete creep deformation in the early stage mainly comes from the viscoelasticity deformation (Mei and Wang 2020). Meanwhile, the long-term creep behavior can affect the safety and serviceability limit state of concrete structures. Concrete creep, intertwined with damage or cracking, gradually leads to damage evolution and stress redistribution of structural components. One notorious engineering example is the continuously vertical deflections of large-span prestressed concrete bridges (Yu et al., 2012;Tong et al. 2016). Massive concrete cracks at box-shape segments′ bottom slabs and webs appear phenomenally together, which indicates complex interactions between concrete creep and cracking (Tong et al. 2016).
As for concrete creep, microprestress-solidification (MPS) theory successfully describes concrete creep considering the processes of moisture transport and heat transfer (Boumakis et al., 2018). With the MPS theory, the solidification theory assumes the concrete aging as the volume growth of gels related to a hydration degree (Rahimi-Aghdam et al., 2019), and the microprestress theory assumes the relaxation of self-equilibrated stresses in the solid nanostructure of cement gels governing the long-term creep, drying creep (Pickett effect), and transitional thermal creep (Bažant et al., 1997a;Bažant et al., 1997b). The MPS theory supersede empirical formulas in specifications (i.e., ACI, B3/B4, and CEB-FIP models) in terms of the multiphysics approach, which effectively accounts for moisture transport and heat transfer. Consequently, the effects of varying temperature or humidity on concrete creep can be captured.
Accurate prediction of time-dependent responses of these concrete structures necessitates a detailed analysis of interactions between concrete creep and cracking. Some contributions are made to take these interactions in numerical analyses and they are comparatively presented in Table 1. To name a few, the smeared crack approach was combined with the viscoelastic behaviors for concrete in De Borst (1987). Concrete creep subject to uniaxial compression was modeled in Mazzotti and Savoia (2003) and Challamel et al. (2005) with an isotropic damage law being taken into account. Moreover, the viscoelastic behaviors of concrete, realized through a Maxwell chain, were coupled to a microplane constitutive model in Di Luzio, (2009). They further elaborated the concrete creep modeling with the microprestress-solidification theory realized through Kelvin chain (Di Luzio and Cusatis 2013). Recently, Luzio′s continuum method was replaced with a discrete element model in the work by Boumakis et al. (2018). Similar work was also found in Abdellatef et al. (2019). Although exhibiting sufficient accuracy, the lattice discrete particle model adopted in Boumakis et al. and Abdellatef et al. (2019) essentially is not a continuum model and is not suitable for structural analysis. Li et al. (2019) realized the nonlinear time-dependent analysis of prestressed concrete bridges considering cracking, creep, and shrinkage within the framework of Abaqus/Standard. The implicit solver was adopted in their study through viscous regularization provided by Abaqus/Standard. More recently, the coupled effects between creep, damage, and plasticity for concrete were taken into account in Ren et al. (2020). Unfortunately, the concrete damage plasticity model in Li et al. (2019) and anisotropic damage model in Ren et al. (2020) is still a local constitutive law without considering nonlocal effect. Tong et al. (2021) tentatively coupled the localizing gradient damage model to the extended microprestress-solidification theory. The semi-implicit algorithm was adopted and four-field evolutions made this method difficult to be implemented within FE framework.
Although exhibiting apparent advantages, the above literatures still have limitation in the numerical analyses, as follows: • Most of the above analyses are implemented with explicit finite element (FE) solver, and the computational cost is prohibitively high for long-term analyses of large-scale structures. • The softening behavior of concrete is not properly regularized, and, therewith, it would lead to unrealistic results due to meshsensitivity problem (Peerlings et al. 1996). • Most of creep models in the above analyses adopt empirical formulas in specifications and cannot reflect the real environmental conditions, of which varying temperature or humidity cannot be ignored.
It is acknowledged that the continuum-based constitutive law is easy to implement and is appropriate for large-scale concrete structures, especially with implicit finite element (FE) solver (Tong et al., 2016). Nevertheless, coupling concrete cracking with time-dependent behaviors is a nontrivial task. Strain softening in quasibrittle materials is the dilemma that has to been properly treated (Peerlings et al., 1996), which leads to loss of ellipticity of the differential equations and also the so-called mesh-dependent solutions upon element size. A series of regularization methods are continuously proposed targeting for mesh-objective simulations, i.e., the introduction of nonlocality in the constitutive model in either integral forms (Bažant et al., 1984), or gradient-enhanced damage (GED) forms (Schreyer and Chen, 1986;Peerlings et al., 1998;Poh and Sun, 2017), and the so-called phase-field approach (Kou et al., 2019; (Bažant et al., 1984). Among them, the phase-field approach to the progressive events of brittle and quasibrittle fracture is proven to capture all stages of fracture, such as crack nucleation, initiation, propagation, and coalescence (Tanne et al., 2017). However, the computational cost of the phase field approach is prohibitively high and is not suitable for structural analyses hitherto. In this study, the GED model incorporating the high-order deformation terms is utilized. It is admitted that many efforts are made continuously to improve the original one proposed by Peerlings et al. (1996), i.e., the localizing GED model (Poh and Sun 2017) and stress-based GED model (Vandoren and Simone 2018). These improvements effectively overcome the spurious damage growth of the original GED model at large tensile strain. However, the selection of GED models is out of the scope of this study. We deliberatively select the original GED model in Peerlings et al. (1996) as its FE implementation is more easy compared to other GED models.
The proposed framework (termed as "GED-MPS" model hereafter) integrates the GED model for concrete′s mechanical behaviors with the MPS theory describing concrete's creep behaviors. The model is implemented within the general FE software Abaqus/Standard with the implicit solver. In detail, Sect-2 introduces the theoretical background and computational framework of the GED-MPS model. Sect-3 illustrates the FE implementation of the model within Abaqus/ Standard with implicit solver. Sect-4 reveals that the model can obtain mesh-objective solutions, not only for mechanical responses but creep deformations intertwined with concrete damage and cracking. Sect-5 validates the proposed model with the experiments by Gilbert and Nejadi (2004), who recorded the long-term deflections and cracking of a series of simply-supported reinforced concrete beams. Finally, conclusions are derived in Sect-6. To facilitate potential researchers and users, parts of the code and the input file are released at https://github.com/TengTongSEU/Coupled-creepdamage.

CONSTITUTIVE MODELING FRAMEWORKS Thermo-Hydro-Mechanical Coupling for Concrete
The GED-MPS model takes an engineering approach to the coupled time-dependent and mechanical behaviors of concrete, aiming to describe the most important aspects with sufficient accuracy at an affordable computational cost. Complete constitutive law of the GED-MPS model follows a rheological representation, as shown in Figure 1, which consists of • a nonaging spring unit representing instantaneous elastic strain tensor ε el ; • a damage unit with the strain tensor ε dam ; • a solidifying kelvin chain unit with typically more than ten elements in series representing the viscoelastic response with aging effect ε v ; • a dashpot unit with viscosity dependent on microprestress S for the long-term creep strain tensor ε f ; • a shrinkage unit for strain tensor ε sh ; and • a thermal unit for strain tensor ε T .
These units are connected in series, with the identical stress tensor σ being transmitted (Yu et al., 2012;Tong et al., 2016). In the absence of significant plastic and viscoelastic strains that may arise at very high confining pressures, the total strain rate _ ε tot is the sum of individual strain rates (Figure 1), as Concrete is interpreted as a composite material in which the coarse aggregates are embedded as inclusions in cement paste. Aging creep occurs exclusively in cement paste, whereas aggregates behave elastically. To capture the mesoscale behaviors of concrete creep in a more realistic manner, multiscale approaches and more sophisticated micromechanical models are necessary, which require more parameters and detailed information on the individual phases (Pichler et al., 2007;Scheiner and Hellmich 2009).
Theoretically, time-dependent evolutions of all relevant material properties should be taken into account to achieve an accurate damage process, i.e., creep and mechanical properties. Nevertheless, the loading of concrete structures at a very earlystage is not the focus of this study. Therewith, the mechanical properties of mature concrete are assumed to be constant, i.e., tensile strength or fracture energy. Their changes with concrete aging are ignored in this work. The multidecade prediction will provide conservative estimations as the slight increase in the mature concrete's mechanical properties is not considered. The local age-dependent fracture properties should be enriched for the framework for concrete structures subject to early-age loadings.

Creep Evaluation by Microprestress-Solidification Theory
To predict creep deformations, semiempirical models are generally adopted by the engineering community, i.e., ACI, B3/B4, and CEB-FIP models. Numerous material properties act as input variables for these models, i.e., compressive strength, water to cement ratio, cement content, temperature, humidity, etc. Afterward, these engineering creep models predict the time-dependent deformation on a cross-sectional level in accumulated form (Boumakis et al., 2018). Concrete creep can be mainly categorized into the aging of concrete, drying creep (Pickett effect), and transitional thermal creep (Bažant and Jirásek 2018). Bažant et al. (1997a,b) and assume that the continuous formation of C-S-H gels is responsible for the aging of concrete, including short-term chemical aging and long-term nonchemical aging. The first two creep sources, namely concrete aging and drying creep, can be formulated in the B4 model (Hubler et al., 2015), as the compliance function J(t, t′): where t is the current time in days, t′ is the age in days when the sustained stress is exerted, J b (t, t′) is the basic compliance function, J d (t, t′) is the compliance function for drying creep, and q 1 is the instantaneous compliance function (q 1 0.4/E 28 1/E 0 , E 0 is the instantaneous elastic modulus without any creep effect, and E 28 is the 28-day Young's modulus). Note that all the empirical models often yield the predictions far away from the measured responses, which can be attributed to the intrinsic heterogeneity of concrete and the variability in mix designs and environmental conditions. On the other hand, these engineering models lack the capacity to accurately predict the local point-wise response at varying stress/strain, local temperature, moisture content, and curing degree. To this end, the MPS theory is employed to describe the concrete's timedependent behaviors.

Solidification Theory for Aging Viscoelasticity in a Rate-Type Formulation
Concrete creep is found to follow the rules of aging linear viscoelasticity (Bažant et al., 2012a;Bažant et al.,2012b), within the service stress range (σ ≤ 0.4f ′ c , and f ′ c is the compressive strength of concrete). As a consequence, Volterra integral equation can be utilized to evaluate creep deformation under a general stress history σ(t). Unfortunately, the compliance function, which is not of a convolution type due to the phenomenon of concrete aging, requires the history variables in all previous time steps to be stored for the current time step analysis. The computation is prohibitively expensive to predict the long-term performance of large-scale structures. Alongside, it is difficult to couple these memory-dependent variables with other memory-independent phenomena, e.g., concrete cracking, corrosion, and steel relaxation (Teng et al., 2017).
A rate-type formulation can overcome these obstacles. In the rate-type approach, the previous history can be fully represented by internal variables only in the last time step (Yu et al., 2012). Rheological model is popularly adopted for rate-type concrete creep (Jirásek and Havlásek 2014;Yu et al., 2012). To this end, A K chain is employed herein, consisting of more than 10 K units, of distinctive modulus D u (μ 1, 2, /, N) and retardation time τ u (μ 1, 2, /, N); see Figure 1 for details. If the τ μ is arranged to own an infinitely close spacing, the continuous retardation spectrum forms (Bažant and Xi 1995). Widder proposed an approximate inversion formula to uniquely identify the continuous retardation spectrum (Widder 1971). The analytical solution for a given compliance function yields where C (k) is the k-th order derivative on time t of the creep part of the compliance function. In this work C(t, t′) is chosen as with q 2 and q 3 as the adjustable parameters for the viscoelasticity based on the solidification theory, Q(t, t′) is a function that can be referred to the B4 mode, n 0.1 is the fixed exponent, and t − t′ is the duration since the load application t′ and is replaced with kτ u in Eq. 3. In the current application, the limit in Eq. 3 needs not to be realized and it suffices to use k 3 (Yu et al., 2012). A discrete approximation of continuous spectrum gives the discrete spectrum: which corresponds to a discrete Kelvin chain and is required for numerical computations. The exponential algorithm is generally adopted to quantify the time increment Δt, which would effectively avoid the numerical instability problem frequently occurring in the central or backward difference method or Runge-Kutta method.
Utilizing the Kelvin units, one can express the effective incremental modulus E″ of concrete at the current time step t n , as where E 0 is the instantaneous modulus, N is the total number of Kelvin units (see Figure 1), and λ μ is determined by the time increment Δt and the retardation time τ u , as Finally, the solidifying Kelvin chain gives the short-term creep strain tensor rate Δε v at the current time step t n , as where γ (n−1) μ is a state variable at the last time step and it should be updated for each Kelvin chain unit at the end of the current time step for each integration point, as with the rate-type formulation at hand, the values of c u , σ, and ε v in all the previous time steps (from 1 to n−1) need not be stored, which are required in the conventional integral-type creep analysis. These history variables can now be represented only by the state variable c (n) μ in the current time step. The continuous spectrum method and exponential algorithm enable the implementation of the rate-type formulations for creep analysis; for more details, refer to Yu et al. (2012).

Microprestress Theory for Viscous Flow
The rate of viscous strain tensor _ ε f originates in the slippage between adsorbed water layers at the microscale level and can be described by the microprestress theory. Effects of temperature and humidity on the microstructure can be described by introducing three transformed time variables: the equivalent age t e describing the degree of hydration, the reduced time t r characterizing the changes in the rate of bond breakages and restoration on the microstructural level, and the reduced microprestress time t s .
Under standard conditions (i.e., room temperature and sealed specimens), all the three time scales, namely t e , t r , and t s , are equal to the actual age of concrete t. For the general condition, the rates of the transformed time are defined as the product of two functions, which respectively characterizes the effects of temperature and humidity histories. The derivatives of the three transformation times with respect to the actual age of concrete t are (Jirásek and Havlásek 2014) where T is the absolute temperature, h is the humidity, R is the universal gas constant, and Q e , Q r , and Q s are active energies for the hydration, viscous process, and microprestress relaxation, respectively. The bonds across the slip plane representing the nanopore filled with hindered adsorbed water are subject to two stresses: the macroscopic applied stress σ causing shear slip and the tensile microprestress. The rate of flow strain tensor _ where the effective viscosity η is a decreasing function of S. The source of _ ε f lies in the relaxation of disjoining pressure and the rupture of atomic bonds and is formulated by microprestress S. The concept of microprestress is useful for the theoretical justification of evolving viscosity η. With the temperature T and humidity h being taken into account, the drying creep effect and transitional thermal creep can be reflected in the viscous flow strain rate _ ε f . However, the microprestress cannot be directly measured and the calibration of microprestress relaxation is difficult. To this end, by some manipulations and rearrangements, Jirásek and Havlásek (2014) obtain the evolution of the viscous dashpot, as with in which c 0 and k 1 are the constant parameters and P is an exponent usually set equal to 2. The internal variable κ T was introduced to keep track of the previously obtained maximum temperature and to adjust the temperature variation on creep and is defined as where κ Tm and κ Tc are material parameters and κ Tm 0.017 and κ Tc 0.001 are recommended in Jirásek and Havlásek (2014). For a standard choice p 2, the above equations are simplified to By assuming that the evolution of viscosity should be the same as in model B4, the initial condition for differential Eq. 17 can be simply postulated as where q 4 is the variable governing the long-term creep in the B4 model.

Local Damage Model
With the assistance of the effective stress concept (Rokhgireh and Nayebi 2019), the stress-strain relation of quasibrittle materials can be expressed as where σ is the second-order Cauchy stress tensor, d is the scalar damage variable (0 ≤ d ≤ 1), D e is the fourth-order elastic stiffness tensor, and ε′ is the second-order elastodamage strain tensor. The ε′ is the sum of the elastic strain tensor ε el and damage strain tensor ε dam . It can be written in the rate form as To describe the damage evolution, the local scalar-valued equivalent strain ε is adopted to evaluate the elastodamage strain tensor ε′. The details expression of ε(ε′) would be given in Eq. 29. For strain-based damage models, a history variable κ is defined as the maximum ε ever obtained in the previous time steps to guarantee the damage irreversibility as The loading function f ( ε, κ) is introduced to complement the stress-strain relation in Eq. 20, as The loading function f ( ε, κ) and the rate of the history variable, Δκ, have to satisfy the discrete Kuhn-Tucker loadingunloading conditions: The history variable κ describes the damage variable d(κ) and damage irreversibility is automatically satisfied. The κ starts at a damage threshold level κ 0 and is updated by the requirement during the damage growth f 0 (De Borst and Verhoosel, 2016).

Nonlocal Gradient-Enhanced Formulations
In a conventional local damage model, κ is related to the local scalar measure of the strain tensor ε. Nevertheless, in a nonlocal generalization, the (local) equivalent scalar strain ε is replaced by a spatially averaged (nonlocal) equivalent strain ε, which is computed in each material point approximately by the differential form (Peerlings et al., 1996): where ∇ 2 denotes the Laplacian operator and c is a positive gradient parameter of the dimension length squared. The parameter c follows a specific form of the weight function and averaging volume and provides the internal length scale which is needed to regularize the damage localization. Hereafter this parameter is treated as an intrinsic material property and should be carefully calibrated to characterize the material's meso-or microstructure. By integrating the averaging partial differential equation Eq. 25 over the entire domain and adopting Gauss' divergence theorem on the term c Ω ∇ 2 ε dΩ, we can obtain where n is the unit directional vector perpendicular to the boundary.
To guarantee that the average of ε over the entire domain equals that of ε, the natural boundary condition is adopted, resulting in the vanishment of the second term in the left-hand side of Eq. 26.

Damage Evolution
Two relations should be well defined to characterize the macroscopic stress-strain behavior: the evolution law of damage d(κ) and the local equivalent strain ε. The evolution of the damage variable d(κ) is critical to govern the postpeak softening behavior and the degree of brittleness. The damage law proposed in Peerlings et al. (1996) is adopted for simplicity: where κ 0 is the threshold to initiate the damage and α and β are material constants governing the softening behavior. The local equivalent strain ε adopts the modified von Mises equivalent strain, which originates from plasticity models for polymers and proposed (Vogler et al., 2013), resulting in where ε′ represent the summation of elastic and damage strain tensors in Eq. 21, I 1 is the first invariant of strain tensor ε′, J 2 is the second invariant of deviatoric strain tensor ε′, and v is the Poisson's ratio. An important feature of this definition of local equivalent scalar strain is the parameter k which is defined as the ratio of the compressive strength f c to the tensile strength f t of the material, k f c /f t . Compressive failure would vanish as k goes to infinity.

Coupling Damage and Creep
To describe the nonlinear concrete time-dependent behavior, a stress-strain law based on the rate independent isotropic damage model is formulated in Eq. 20, which can be interpreted as Hooke's law linking the strain to the effective stress tensor σ(Bažant and Jirásek 2018): The above equation indicates the stress transmitted by the undamaged continuous solid material, which would be σ, rather than σ with defects being excluded. Considering that the creep stress is much lower than the cracking threshold, this nonlinear creep model is only suitable for the undamaged concrete, not for a cracked body. Thus, the creep strain calculation, similar to the plastic strain, is based on the effective stress tensor σ, rather than the nominal stress tensor σ. As a consequence, the short-term creep strain tensor ε v and the long-term creep strain tensor ε f , obtained in sections Solidification Theory for Aging Viscoelasticity in a Rate-Type Formulation and Microprestress Theory for Viscous Flow based on nominal stress tensor σ, should be amplified by the factor 1/(1 − d).

FINITE ELEMENT FORMULATION Governing Equations and FE Discretization
The model presented in this work is discretized using a Galerkin FE approach. Strong forms of governing equations for the GED-MPS model for quasistatic problems are where b represents an external load vector. The two unknown fields to be calculated are the displacement field u(x, t) and the socalled nonlocal equivalent strain field ε(x, t). These equations are supplemented by the following Neumann and Dirichlet boundary conditions: where t * and u p are prescribed force and displacement at the boundary. Integrating Eqs 31, 32 over the domain Ω and weighting with w u and w ε respectively yields the weak form after applying the divergence theorem and the boundary conditions: In the FE implementation, the displacement field u and nonlocal equivalent strain field ε at an arbitrary point is interpolated by their nodal values (Sarkar et al., 2019): with ṵ and ε being the nodal degrees of freedom, N u and N ε represent the shape functions, and B u and B ε are their partial derivatives with respect to the spatial coordinates. According to the Bubnov-Galerkin method, the corresponding weight functions w u and w ε can be discretized analogously. Consequently, the discretized equilibrium equations of Eqs 34, 35 reduce to To construct a consistent incremental-iterative Newton-Raphson solution procedure, Eqs 38, 39 have to be linearized. The linearization at iteration i with respect to the previous iteration i − 1 is outlined as Of special attention is that the components in the stiffness matrix and the right-hand side vector depend on the ε′ in Eq. 21. The submatrices are defined as and the subvectors in the right-hand sides are defined as Note that the element stiffness matrix becomes nonsymmetric since K uε i−1 is not equal to the transpose of K εu i−1 . The implicit solver with Newton-Raphson algorithm could guarantee the quadratic rate of convergence, even if the nonsymmetric tangential stiffness matrix in Eq. 40 is utilized (Peerlings et al., 1996).

Creep Strain Tensors
The key parameter in implementing the GED-MPS model is to obtain the local equivalent strain field ε(ε′), which is a function of the strain tensor ε′. To this end, now we turn our attention to the numerical treatments of increments of creep strain tensors, namely the aging viscoelasticity Δε v and the viscous flow Δε f .
The aforementioned rate-type formulation can be used for the basic creep compliance J b (t, t′), as presented in Eq. 4. Taking derivative of J b (t, t′) with respect to time gives where reduced time t e (t) and t r (t) are introduced to account for general temperature and humidity conditions, v(t) describes the volume growth function with the solidification process, and Φ(t, t′) is the strain in nonaging solidified matter whose material properties are assumed to be invariable with time.
The v(t) and Φ(t, t′) are defined as As a consequence, the exponential algorithm following the procedures can be implemented to obtain the increment of aging viscoelasticity Δε v : Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 701458 1) At t t 0 , where t 0 indicates the time with the external load being applied, initialize the internal variables γ (0) μ . Select τ u 10 −7+u , with u 1, 2, 3,. . ., 14.
2) Let ξ t − t′, and the continuous spectrum is calculated as with 3) Discretize the spectrum using 4) Calculate λ u according to Eq. 7, and calculate E″ according to Eq. 6. 5) Obtain increment of aging viscoelasticity Δε v , according to Eq. 8. 6) After retrieving the stress increment Δσ by the GED model, update the internal variable c (n) μ ; refer to Eq. 9. 7) Begin the next time step. Now we turn to the increment of viscous flow Δε f , which needs the viscosity Eq. 17 with initial condition Eq. 19. Suppose that the value of viscosity η n−1 at time t n−1 is known from the previous time step or from the initial condition, we wish to obtain the value of viscosity η n η n−1 + Δη at time t n t n−1 + Δt. The exponential algorithm (Bažant 1971) following the procedures is implemented to obtain the increment of viscous flow Δε f .
1) The temperature and humidity of the previous time step T n−1 and h n−1 are known, as well as their increments ΔT and Δh, from the external heat transfer and moisture transport analyses. The viscosity η n at the current time step is explicitly evaluated as for arbitrary large Δt, as 2) Suppose Δ η ≪ η n , the increment of viscous flow Δε f is evaluated as For the creep deformation coupled with concrete damage or cracking, the nominal stress tensor σ or its increment Δσ should be replaced by their effective counterparts σ and Δσ, respectively, to evaluate Δε v and Δε f . Afterwards, the increment of strain tensor Δε′ can be obtained by subtracting Δε v and Δε f from the increment of total strain tensor Δε tot , as well as the increments of hygral and thermal strain tensors (Δε sh and Δε T ). The Δε′ is  subsequently adopted to obtain the stress tensor σ n and damage indicator d n at the current step time t n .

Implementation Aspects in Abaqus
The system of equations above is highly nonlinear. To implement the model with implicit FE solver, the general FE software Abaqus/Standard is selected for its built-in implicit incremental-iterative Newton-Raphson algorithm and automatic time-step-ping schemes. The user-defined element (UEL) subroutine enables us to define the element, of which computation of element tangent stiffness and nodal force vectors is realized.
Plane stress is taken into account in the following simulations. To guarantee the consistency requirements, the shape functions and their derivatives for the displacement field u(x, t) are defined over 2D quadrilateral element of eight nodes, whereas these terms are defined over 2D quadrilateral element of four nodes for nonlocal equivalent strain field ε(x, t) (Sarkar et al., 2019). These two element types share the same nodes at the four corners; see Figure 2 for details. Consequently, the nodal displacement vector u ∼ and nodal nonlocal equivalent strain vector ε ∼ are expressed as The corresponding shape functions and their derivatives are and It is noteworthy that the analysis with UEL subroutine is inconvenient in the postprocessing and visualization of the results. Because shape functions are defined by users in the UEL subroutine, the Abaqus cannot extrapolate variables from Gauss points to the element nodes, automatically. To this end, an auxiliary dummy mesh is adopted, consisting of standard Abaqus elements that resemble the UEL elements in terms of number of nodes and integration points. The material response at each integration point in the auxiliary mesh is defined using a user material subroutine (UMAT), which enables the user to define the constitutive matrix and stresses from the strain values.
In this auxiliary mesh, the stress components are deliberatively set as zero so as not to influence the global solution. The data from the UEL for each time increment we want to observe in Abaqus/ Viewer is stored as built-in array SVARS, which allows transferring to the UMAT subroutine by the built-in array STATEV for each corresponding element and integration point. Transferring of values from SVARS array to STATEV array is accomplished by making use of the common statement (Msekh et al., 2015). Figure 3 shows the technique to implement the visualization of the analyses with the UEL subroutines.

ILLUSTRATION OF GED-MPS MODEL
A 100 mm long bar is taken for an illustration, which is subject to uniaxial tensile loading. The Young's modulus is deliberatively reduced by 5% in a 10 mm wide zone in the middle of the bar to trigger localization of deformation.
The first group of analyses embraces the bar discretized with 80 elements, but with three different values of gradient parameter c. Profiles of the damage indicator d corresponding to different c are depicted in Figures 4A,B. For the case with c 1 mm 2 , a clear narrowing of the localization zone is observed, in terms of d. This indicates that the case with c 1 mm 2 is very close to the local damage model, which is effective for brittle materials, but quasibrittle materials. On the other hand, profiles of damage indicator d are different for the other two cases with c 1 and 5 mm 2 , as nonlocality is clearly depicted as the damage overflowing the 10 mm imperfection zone. The value of gradient parameters c clearly affects the width of damage zone  in the middle of the bar. Alongside, the effects of gradient parameters c on the force-displacement relations of the bar are also revealed as the normalized stress-strain curves in Figure 4C. Obviously, with the increase of c, the bar behaves more ductile and the load-carrying capacity also increases slightly. The second group of analyses embraces the bar discretized with different element size, but with identical gradient parameter c 5 mm 2 . Profiles of the damage indicator d corresponding to different element size are depicted in Figures 5A,B. Clearly, although with different element size, all the three cases converge to an identical result, indicating that the GED model satisfactorily eliminates the mesh-sensitive problem frequently encountered in the softening region of quasibrittle materials. Nonlocality is clearly demonstrated with nearly the same three damage zones with the width bigger than the 10 mm imperfection zone. Alongside, effects of element size on the force-displacement relations of the bar are also revealed as the normalized stress-strain curves in Figure 5C. Only a slight difference is observed for the three cases. Figure 5 clearly shows that mesh-dependent solutions of mechanical responses of quasibrittle materials are regularized by the GED model, without creep deformation. However, limited literatures report the mesh-sensitivity problem in the time-dependent analyses of quasibrittle materials, especially when coupled with cracking, according to the authors' knowledge. Intuitively, meshsensitivity solution is a concern, as the creep depends on the stress tensor and the damage state at each material point.
To this end, the proposed GED-MPS model is applied to the same bar. The external loading is applied at the bar's right end at the day of t′ 10. The external displacement is set as 0.08 m and the loading is finished within 0.01 day, and it is kept constant for another 1 day. Loading rate effect, aging effect, shrinkage, and thermal strains are not considered for simplicity. The following parameters are selected: q 1 25 × 10 -6 MPa −1 , q 2 200 × 10 -6 MPa −1 , q 3 20 × 10 -6 MPa −1 , q 4 0 MPa −1 , T 293 K, and h 1.
Similarly, the first group of analyses embraces the bar discretized with identical 80 elements, but with different gradient parameters as c 1, 5 and 15 mm 2 . Profiles of creep strain ε′ in the horizontal direction for the three cases are illustrated in Figure 6A, with effects of damage being involved. For the case with c 1 mm 2 which approaches local damage analysis, the localized creep-induced strain, with peak value of ε′ 3.53 × 10 -2 is observed. For the other two cases, with c 5 and 15 mm 2 , nonlocality is clearly captured for the creep strain ε′, which overflows the 10 mm imperfection zone. Nevertheless, the peak value of ε′ decreases to 1.08 × 10 -2 for the case with c 5 mm 2 and 0.28 × 10 -2 for the case with c 15 mm 2 , respectively.
Alongside, effects of gradient parameters c on the relationship between time and normalized stress are revealed in Figure 6B. The time from day 10.0 to day 10.1is the loading time and various strain-softening responses are captured, with analogy to the phenomenon revealed in Figure 4C. Moreover, relaxations of normalized stresses are illustrated from day 10.1 to day 11.1 in Figure 6C, with effects of c. Interestingly, apparent meshsensitivity is witnessed for the three cases. Although with the peak value of creep strain ε′ 3.53 × 10 -2 , the relaxation of normalized stress is not significant for the case with c 1 mm 2 . It looks like that the localization of creep strain does not significantly affect the relaxation of normalized stress. On the other hand, the loss of normalized stress approaches 0.29 MPa for the case with c 5 mm 2 and 0.75 MPa for the case with c 5 mm 2 ; see Figure 6C. Obviously, the gradient parameter c smears the damage into the neighboring zone, so as for the creep strain intertwined with damage or cracking.

CALIBRATION WITH TESTS OF SIMPLY-SUPPORTED CONCRETE BEAMS
Engineering community concerns more about long-term deflections of reinforced concrete or prestressed concrete structures, coupled to concrete cracking. In this section, the capability of the GED-MPS model is further extended to the reinforced concrete structures. Longterm tests of simply-supported reinforced concrete beams in Gilbert and Nejadi (2004) are simulated.
Although totally 6 beams and 6 slabs were reported, this paper selected two typical specimens for simulations, namely specimens  B2-a and B2-b. These two specimens owned a 250 × 340 rectangular cross-section and 3,500 mm in length. The concrete cover thickness is 25 mm. Two 16 mm diameter rebars were arranged as the flexural reinforcement. During the test, these two specimens were loaded to 50% (specimen B2-a) and 30% (specimen B2-b) of their flexural strengths, respectively. A constant external loading was applied on the concrete specimens. This loading test began at 14 days and lasted for around 400 days. Concrete properties were measured at different ages and they are summarized in Table 2.

Creep Calibrations
The concrete creep and shrinkage tests were provided in Gilbert and Nejadi (2004), which begun at the 14-day concrete age. No data is provided regarding temperature and humidity, and they are set as T 293 K and h 0.6, respectively. The values are   With respect to parameters for concrete creep, the following creep parameters are selected to fit the experimental creep data, with q 1 47.4 × 10 -6 MPa −1 , q 2 200 × 10 -6 MPa −1 , q 3 20 × 10 -6 MPa −1 , andq 4 46 × 10 -6 MPa −1 ; see Figure 7B. In addition, the following fracture parameters are adopted according to Table 2, with E 28 22.8 GPa, v 0.18, κ 0 0.000122 (f t 2.8 MPa), β 300, α 0.99, and k 24.8/2.8 8.5. The E 0 is the Young's modulus at concrete age of 14 days, and it would rise according to data in Table 2 concerning concrete aging effect.

Long-Term Deflections and Cracking Distributions
The GED-MPS model is verified with experimental data of specimens B2-a and B2-b, which were subject to external loading of 18.6 and 11.7 kN, respectively. The identical calibrated parameters listed in Sect-5.1 are utilized without any modification. Some artificial defects were deliberatively introduced at the bottoms of FE models to trigger the cracks.
ffects of different gradient parameter c 4 and 14 mm 2 on the long-term behaviors are investigated with the other parameters being the same. Comparisons of FE and experimental crack profiles for specimens B2-a and B-2b are illustrated in Figures  8A,B, respectively. For the FE results with c 4 mm 2 , crack patterns correspond well with the test in height, region, and spacing of cracks. On the other hand, cracking profile for the case with c 14 mm 2 yields less accurate prediction in the cracking height. Absent data are provided in the report regarding the cracking profiles at the commencement of the test.
Long-term vertical deflections of Specimen B2-a with different c are compared with experimental data in Figure 9. Both simulation results yield satisfactory agreements. Interestingly, little difference is observed for the deflections with different c. Possible reasons are attributed to the following two reasons. Shrinkage deformation plays a dominant role over creep deformation for the specimen, as merely 50% peak strength is applied and the stress level is low. Another possible reason is that although higher c reduces the damage height, it broadens the damage zone with respect to the case with lower c.
Furthermore, conventional linear viscoelastic analysis considering only creep and shrinkage deformations, without intertwined effect with concrete damage and cracking (termed as MPS in Figure 10), is compared to the GED-MPS model in Figure 10A for specimens B2-a and Figure 10B for specimen B2b. Clearly, conventional linear viscoelastic analysis underestimates their long-term vertical deflections, which can be satisfactorily remedied by the GED-MPS model.
Considering the effect of different temperatures on the long-term deflection of these specimens, four different temperatures are selected, and they are T 0.1°C (273.1 K), 20°C (293.1 K), 40°C (313.1 K), and 60°C (333.1 K). The long-term deflections subject to different temperatures of specimens B2-a and B2-b are presented in Figure 11, while all the other parameters are kept constant.
The effects of different humidity on long-term deflection of specimens B2-a and B2-b are further revealed in Figure 12, with different humidity of h 0.6, 0.8, and 1.0, whilst all other parameters are fixed. It is noteworthy that h 1.0 illustrates the effect of creep on the long-term deflection only, without shrinkage effect.

CONCLUSION
The GED-MPS model is proposed and implemented with implicit FE software Abaqus/Standard. The model integrates the GED model capable of circumventing the mesh-sensitivity in the conventional mechanical analyses for quasibrittle materials exhibiting strain softening behaviors, and the MPS theory capable of predicting point-wise creep responses.
The GED-MPS model anchors at the rate-type formulation, and, therefore, is capable of incorporating other memorydependent processes. After being enriched with continuous spectrum method and exponential algorithm, the model is efficient and stable for large-scale structural analysis. The model successfully regularizes the mesh-sensitivity problem, which is already concerned with mechanical analyses of materials exhibiting softening, but less with the nonlinear creep intertwined with cracking. The GED-MPS model is applied to simply supported reinforced concrete beams. It is proven that the proposed model is capable of capturing not only the long-term vertical deflection trends, but also the time-dependent cracking propagations affected by creep and shrinkage. Further improvements to polish the model could be conducted, i.e., implementing more accurate anisotropic constitutive model for concrete, instead of the isotropic one, replacing the current GED model with stress-based GED model to address the issue of spurious damage growth. With the GED-MPS model being implemented in the implicit FE algorithm, authors are expected to extend this model to accurately predict long-term behaviors of real large-scale creep-sensitive structures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation. To facilitate potential researchers and users, parts of the code and the input file are released at https://github.com/TengTongSEU/ Coupled-creepdamage.