An image-based meso-scale model for the hygro-mechanical time-dependent analysis of concrete

A new computational framework is developed in this paper for investigating the time-dependent behaviour of concrete including creep, shrinkage and cracking. The developed model aims to explain certain aspects of the time-dependent cracking and creep of concrete that cannot be captured using homogeneous models. The model is based on the scaled boundary finite element method, and it is coupled with a quadtree decomposition algorithm which converts digital images of concrete meso-structures into meshes. Concrete is treated as a two-phase composite which consists of elastic aggregates and mortar that is subjected to time-dependent deformation. The basic creep behaviour is treated as viscoelastic, which is modelled based on a rate-type rheological model corresponding to a Kelvin chain. Drying creep is modelled using a viscous unit which depends on the stress level, and drying shrinkage is stress independent. Both drying creep and drying shrinkage are related to the internal humidity. The humidity distribution within concrete is determined using a diffusion analysis. The moisture movement within mortar is governed by a nonlinear diffusion equation, whereas the aggregates are assumed impermeable. The cracking of concrete is explicitly modelled on the meso-scale through coupling of the continuum damage model for cracking within the mortar phase, and the cohesive zone model for debonding between aggregates and mortar. The proposed model is verified by simulating well-documented experimental studies in the literature. The capability of the proposed model in simulating the time-dependent behaviour of concrete and capturing the crack patterns has also been demonstrated.


Introduction
Concrete undergoes time-dependent deformation mainly due to creep and shrinkage, which affects the performance and the design life of concrete structures.Both creep and shrinkage are complex phenomena that are not fully understood yet.Creep represents the continuous increase of deformation under sustained load and for concrete it is generally divided into basic creep and drying creep.Basic creep represents the time-dependent strain developed in a sealed specimen or one that is in hygral equilibrium with the ambient environment.Drying creep is the additional strain developed in a specimen that is subjected to simultaneous drying environment, this is also known as the Pickett effect.Several well accepted mechanisms were proposed for creep including the seepage theory [1], viscous flow [2], and the microprestress-solidification [3] theory.Shrinkage represents the change in volume measured on a load free specimen as a result of drying shrinkage and autogenous shrinkage.For normal strength concrete, the latter one is negligible and total shrinkage is assumed to be equal to the drying shrinkage only which is taken as the strain measured in a specimen that is exposed to a drying environment [4,5].Drying shrinkage is commonly associated with mechanisms include capillary pressure, surface tension, disjoining pressure and interlayer water movement [6,7].Both drying creep and drying shrinkage develop with the presence of a dry atmosphere which is typical under service conditions.Cracking is also commonly seen in concrete structures over time.These cracks can be caused by external loading, or by the time-dependent deformation, especially shrinkage, in inter-nally or externally restrained structural members.Cracking of concrete increases its permeability, making it more susceptible to sulfate [8] and chloride attacks [9,10].
Macro-scale models are typically used for concrete and they consider it as a homogeneous material.Extensive experimental research on creep and shrinkage of concrete has been conducted.A collection of the experimental data in the literature can be found in Hubler et al. [11].By calibrating experimental data, empirical creep and shrinkage models have been developed in design codes, such as the ACI model [4], fib model [12] and B4 model [13] and they allow the prediction of creep and shrinkage from a simple formula.For drying creep and drying shrinkage, these design code models adopted a sectional approach, in which only the mean effects of drying over the cross section of the structural member are considered.Both drying creep and drying shrinkage are assumed to be uniformly distributed within the section.Other more sophisticated models have also been proposed to model drying creep and drying shrinkage [3,[14][15][16][17].Extensive applications of creep and shrinkage models for structural analysis have been made in the past, including on plane specimens [18,19], beams and girders [20][21][22] and other members.
Concrete is a composite consisting of cement paste, fine aggregates and coarse aggregates.The time-dependent deformation observed in concrete is well accepted to originate from the cement paste and it is restrained by the presence of aggregates.The restraint depends on properties of the aggregates, such as their elastic modulus, content, distribution, shape etc.This restraint mechanism which is essential for predicting creep and shrinkage of concrete cannot be explicitly captured by homogeneous models.Attempts have been made to simulate creep and shrinkage using composite homogenisation models [23][24][25][26].These models were developed by making some simplified assumptions on the distribution of aggregates, and they allow the prediction of the deformation of concrete based on the material properties and volume fraction of the constituents.Nevertheless, the random distribution of the aggregates, their size gradation and shape cannot be explicitly modelled.
Meso-scale modelling of concrete has become increasingly popular in recent years due to its capability to explicitly model the different components of concrete.On the mesoscale, concrete is typically treated as a composite of mortar and coarse aggregates.The interfacial transition zone (ITZ) may also be considered with meso-scale models which is a thin layer (around 50μm) of porous region formed around the aggregates due to wall effect [27,28].One key step of meso-scale modelling is the acquisition of the concrete meso-structure.Due to difficulty of obtaining the internal structure of concrete, numerically generated meso-structures have been adopted by most of the studies which was obtained by making certain simplified assumptions regarding to the aggregate shape, size, distribution and etc [29,30].Alter-natively, the internal structure of concrete may be obtained from X-ray microCT scanning, which can be used to run direct image analysis [31,32].
Most of the existing meso-scale models of concrete focused on simulating its crack pattern under load [29,33] without the consideration of the time-dependent deformation.Continuum models based on the finite element method (FEM) are adopted by most of the studies in which the concrete meso-structure is discretised into finite elements with different material properties and constitutive relations assigned to them.The ITZ has been typically modelled as either a thin lay of elements [32,34] or through cohesive zone models [35,36].Other than FEM, discrete modelling approaches such as the beam lattice model [37,38] and the lattice discrete particle model [39] have also been proposed in which the concrete meso-structure is characterised by a grid of beam elements.A comprehensive review of these studies can be found in Thilakarathna et al. [40].Investigations have been made in recent years to model individual aspects of the time-dependent behaviour of concrete, including creep [41][42][43], shrinkage and shrinkage induced tensile cracking [38,44,45].However, only limited attempts [35,46,47] have been made to integrate creep, shrinkage and cracking into a single model.Idiart et al. [35] extended the mesoscale model developed by López et al. [30] to simulate the time-dependent behaviour of concrete.The model accounts for creep using a viscoelastic Maxwell chain and shrinkage as an additional stress-independent strain.The cracking of the mortar phase and ITZ is captured using a cohesive zone model.Despite that cohesive zone models are commonly adopted for ITZ modelling, they may not be suitable for the modelling of cracks within the mortar phase due to their mesh dependency.Havlásek and Jirásek [46] investigated the relation between ultimate shrinkage and the ambient relative humidity using a meso-scale hygro-mechanical model.The creep and shrinkage strain components are modelled using the Microprestress-Solidification theory [3], and cracking of mortar is captured using a scalar damage model.However, the debonding of ITZ has not been accounted for in this study.Another meso-scale hygro-mechanical model is proposed by Ožbolt et al [47].The model accounts for basic creep, shrinkage and cracking.Drying creep was not considered as it is assumed to be caused by the interaction between damage and the time-dependent strain components only, and the debonding of ITZ was not considered either.
Most meso-scale studies have been carried out based on FEM as mentioned above.Meshes generated based on Delaunay triangulation have been adopted by many studies [35,47].However, it can be difficult and time consuming to generate a quality mesh especially when the meso-structure contains aggregates with a large quantity and irregular shapes.To simplify the mesh generation procedure, a uniform finite element mesh was adopted by many studies [32,42].Similar proce-dures have also been adopted in typical beam lattice models for the generation of lattice structure [31,38].Another powerful tool capable of running meso-scale analysis that has not been fully utilised for this purpose is the scaled boundary finite element method (SBFEM).The SBFEM is a continuum based semi-analytical numerical method developed by Song and Wolf [48].The model can be coupled with a fast and automatic quadtree decomposition algorithm that converts digital images of concrete meso-structures into meshes [49].The generated quadtree meshes cannot be directly analysed using FEM due to the presence of hanging nodes.However, the SBFEM is highly complementary with quadtree meshes due to its ability to handle polygonal elements with arbitrary number of sides.By coupling SBFEM with the quadtree algorithm, this method serves as a more efficient alternative for the meso-scale analysis of concrete.
In previous works by the authors, a novel simplified twoscale model based on SBFEM was proposed to predict the creep [50,51] and shrinkage [52] of concrete from the cement paste scale.The two-scale refers to the mortar scale which consists of cement paste and fine aggregates, and the concrete scale which consists of mortar and coarse aggregates.The model was simplified in a sense that drying of the specimens were not explicitly modelled.Instead drying creep and drying shrinkage were assumed to be uniform within the structure.In the present study, a more advanced and integrated hygro-mechanical model based on 2D SBFEM formulation is proposed to account for the coupled effects of basic creep, drying creep, drying shrinkage, aging and cracking simultaneously.The non-uniform internal humidity field within the concrete specimen or member is also considered.In this study, we will focus on the concrete scale.The coarse aggregates are assumed to be elastic, whereas the mortar is subjected to time-dependent deformation including creep and shrinkage.The basic creep response is assumed to be viscoelastic with aging which is simulated using a rate-type Kelvin chain model with time-dependent model parameters.Drying shrinkage and drying creep are addressed using the model proposed by Bažant and Chern [14], in which drying shrinkage is directly linked to the humidity and drying creep is linked to both humidity and stress.The moisture movement in cementitious material is assumed to be driven by a non-linear diffusion equation developed by Bažant and Najjar [53], which leads to the humidity distribution within the structure.In addition to the time-dependent strain components, cracking and softening of mortar under tensile load is accounted for through an isotropic damage model.The modelling of the time-dependent deformation and damage are combined through the use of an effective stress approach.The debonding of ITZ is considered through the cohesive zone model.
In the following sections, the constitutive relations of the materials on the meso-scale are presented in Sect. 2. After-wards, the image-based mesh generation and the SBFEM formulation are introduced in Sect.3.For clarity, section 4 provides a summary of the overall solution procedure for the present hygro-mechanical model.Then in Sect.5, numerical examples are carried out to demonstrate the capability of the proposed model.Finally, this paper is concluded with a summary of the presented findings.

Time-dependent strain components
The proposed meso-scale model accounts for deformation caused by mechanical loading, as well as the change of moisture conditions.Aggregates are assumed to be linear elastic, whereas the mortar matrix is assumed to exhibit timedependent deformation.The total strain in the mortar matrix (t) at a given time t is assumed to be the sum of instantaneous strain ins (t), basic creep strain bc (t), drying creep dc (t) and shrinkage strain sh (t) as expressed in Eq. ( 1).
The rheological model used to simulate the time-dependent deformation of mortar is shown in Fig. 1.The deformation of the rheological model is assumed to be controlled by an effective stress σ .The mortar matrix is assumed to be linear under compression.In this case, the effective stress σ is equivalent to the nominal stress σ .Under tension, mortar undergoes cracking which is captured using a scalar damage model.In this case, the effective stress becomes larger than the actual stress, leading to larger deformation of the entire rheological model as compared to the linear case.Alternatively, the other approach of incorporating damage into a creep analysis of concrete rely on the use of nonlinear viscoelastic models like the modified principle of superposition [54].This approach yields a rheological model with stressdependent spring and dashpot constants and it was used in past studies by the authors [21,55] for the creep analysis of reinforced concrete beams.Nevertheless, it assumes recoverable strain at unloading, which is not a characteristic of the drying creep component of concrete.More details of the damage model is presented in Sect.2.3.
There are four strain components involved in the rheological model as shown in Fig. 1.For conciseness, these strain components are categorised into viscoelastic strain ve (t) and moisture induced strain mis (t).The viscoelastic strain ve (t) consists of the elastic strain e (t) and basic creep bc (t).The viscoelastic behaviour is modelled using a Kelvin chain model consisting of a spring and a number of Kelvin units.The constants of the springs and dashpots vary with time due to aging.The moisture induced strain mis (t) consists of drying creep dc (t) and shrinkage sh (t).Drying creep is treated as viscous and modelled using a dashpot unit.Shrinkage is described with a stress-independent unit.The detailed formulations of these strain components are presented in the subsequent sections.

Viscoelasticity
Following the principal of superposition, the stress and strain relation of a viscoelastic material is expressed as where J (t, t ) is the compliance function that describes the strain at time t caused by a sustained unit stress σ applied at time t .The following form of compliance function is adopted in this study [50] where c 1 and c 2 are constants to be fitted against experimental data and E(t ) is the elastic modulus at loading age.A rate type formulation is adopted to facilitate time-stepping analysis.Following the concepts of continuous retardation spectrum [56], the fitted compliance function is approximated with a Dirichlet series corresponding to a Kelvin chain model (see Fig. 1).The Kelvin chain model is then converted into an incremental form following the exponential algorithm proposed by Yu et al. [57].For a given time increment t k = t k+1 − t k , the incremental stress-strain relation is expressed as where E ve,k is the incremental pseudo modulus for viscoelasticity and ve,k is the strain increment due to basic creep, expressed as The stiffness moduli of the Kelvin chain E 0,k+1/2 and E μ,k+1/2 are evaluated at the middle of each time step t k+1/2 to account for aging, β μ,k and λ μ,k are factors that need to be evaluated for every Kelvin unit, γ μ is the internal variable that needs to be updated at the end of every time step as

Moisture induced strain
Following the work of Bažant and Chern [14], the rate of moisture induced strain ˙ mis is directly linked to the humidity rate ḣ, expressed as where ∞ s is the ultimate free shrinkage at zero humidity, g sh (t e ) characterises the reduction of shrinkage due to hydration and aging.It is defined as the ratio of the elastic modulus at the start of drying t 0 over the elastic modulus at the hydration time t e g sh (t e ) = E(t 0 ) The hydration time t e depends on the moisture condition and its formulation will be discussed later in Sect.2.2, k h is a function that defines the relation between the normalised shrinkage and humidity.Under typical environmental conditions where the humidity is generally above 50%, function k h exhibits an approximately inverse linear relationship to the humidity with the expression of k h = 1 − h.By substituting this into Eq.( 8), it can be simplified to ˙ mis = k sh g sh (t e ) 1 + r σ sgn( ḣ) ḣ (10) in which k sh is the shrinkage ratio defined as the shrinkage strain developed per percentage of humidity drop.Equation (10) may be rearranged to represent the drying shrinkage sh (the stress-independent part) and drying creep dc (the stress-dependent part) component respectively, expressed as ˙ dc = k sh g sh (t e )r σ sgn( ḣ) ḣ By assuming the factor g sh (t e ) to be a constant within each time step, the incremental shrinkage strain sh,k can be simply derived from Eq.(11) as sh,k = k sh g sh,k+1/2 h k (13) where g sh,k+1/2 is evaluated at the middle of each time step t k+1/2 .For drying creep (Eq.12), its rate form is directly related to the stress which is equivalent to a dashpot unit (see Fig. 1) with viscosity of η = (k sh g sh (t e )r sign( ḣ) ḣ) −1 .The incremental formulation of drying creep strain may be k sh g sh (t e )r σ sgn ḣ ḣ dt e (14) where E dc,k is the incremental pseudo modulus for drying creep and dc,k is the incremental drying creep strain.They are expressed as

Integrated incremental formulation
The overall incremental stress-strain relation may be obtained by summing and rearranging Eqs.(4), ( 13) and (15), which leads to where E k and k are the overall incremental pseudo modulus and incremental strain respectively, expressed as

Aging and hygral effects
The mortar phase undergoes a continuous aging process that leads to an increase of its elastic modulus and tensile strength over time.Their variation is assumed to take the following form [12]: where s is a constant that governs the magnitude of aging.
In order to account for the effect of humidity on the aging and viscoelastic process, Bažant [58] proposed the introduction of a hydration time t e and a reduced time t r respectively.The hydration time t e replaces the actual time in Eqs. ( 21) to (23), and the reduced time t r is used in the incremental form which replaces the actual incremental time used in the computation of the factors β μ,k and λ μ,k .They are adopted in this work as follows [58]  4  (24) where β eh characterises the effect humidity on aging, and β ch characterises the the effects of humidity on basic creep; α e and α r are material constants that can be calibrated from experimental data.Following the trapezoidal rule, the incremental form of the hydration time t e,k and reduced time t r ,k may be derived as

Damage criterion
As mentioned in Sect.2.1, the instantaneous stress-strain relation of mortar in compression is taken as linear elastic with aging material properties.In tension it is nonlinear due to cracking.In the present study, a continuum isotropic scalar damage model is adopted to model cracking, with the following stress-strain relation [59] Fig. 2 Stress-strain diagram with exponential softening where ω is a variable that characterises the degree of damage or stiffness reduction.By rearranging this equation, we can see that Hooke's law still applies for a damaged material, except that it links between the effective stress and the strain.
The effective stress σ is expressed as The damage variable ω takes values between 0 when undamaged to 1 when completely damaged.The damage evolution is assumed to follow a typical exponential softening curve [60] (see Fig. 2) where 0 = f t /E is the limit elastic strain under uniaxial tension which changes with time due to aging, f is a parameter that affects the softening response and κ = max{˜ (t) : t = 0, . . ., t k } is the maximum equivalent strain reached in the loading history assuming irrecoverable damage in tension.
The equivalent strain ˜ adopted here is based on the Rankine's criterion of maximum principal stress where E is the elastic modulus and σI denotes the positive parts of the effective principal stress.In the case of unloading or reloading, the stress-strain relation follows a linear function between the origin and the point on the stress-strain diagram corresponding to the maximum equivalent strain κ.

Non-local damage regularisation
It is well accepted that local continuum damage models are sensitive to the domain discretisation.To circumvent this issue, the integral-type non-local damage model [60] or the crack band model [61] have been typically adopted.The former approach is adopted in this study which is based on replacing a state variable with its non-local counterpart through weighted averaging over a given space.In this approach, a local field f (x) in a domain V can be replaced by a non-local field f (x) as where α(x, ξ) is a non-local weight function, defined as α 0 is a monotonically decreasing non-negative weight function of the distance r = ||x − ξ || from each point in consideration.A truncated quartic polynomial function [60] is adopted here: in which R is the interaction radius corresponds to the largest distance from point ξ that affects the non-local average at point x.

Damage with aging material properties
The continuum damage model has been typically applied to cementitious material for short-term analysis.Since the material properties do not vary much in a short period, the damage evolution is governed by a static stress-strain curve as presented in Sect.2.3.1.However, in a time-dependent long-term analysis, such an assumption no longer holds due to aging.As described in Sect.2.2, both the elastic modulus and tensile strength of the mortar increase with time.In the context of incremental analysis, it is assumed that all the parameters are constant within each time step with their values evaluated at the middle of the time step.For two consecutive time steps k and k +1, the stress-strain relation for a mortar specimen that has not been previously loaded is assumed to be governed by its elastic modulus and tensile strength at the middle of these two time steps respectively as shown in Fig. 3a.
In the case when the mortar has been previously loaded, we assume that both the stress and damage state of the material remain unchanged due to aging when moving from time step k to k +1. Figure 3b illustrates the load path of the mortar that carries load from a previous time step.When a change of load happens at time step k + 1, it follows the same aging stressstrain curve presented in Fig. 3a except that the curve has been shifted horizontally to coincide with its previous stress and damage state.In this case, the maximum equivalent strain at the current time step κ k cannot be simply adapted from the previous step.Based on the assumption that the damage state does not change due to aging when moving to a new time step, the starting position of the maximum equivalent strain at time step k + 1 is determined from the damage index at the previous time step ω k by rearranging Eq. (30) as where W ( ) is the Lambert W function.It can also be seen from Fig. 3b that when the material is unloaded back to a stress-free state, it still carries some residual strain.The equivalent strain ˜ k for an aging material is still determined using Eq. ( 31) except that the constant elastic modulus is replaced with the one at the middle of the current time step E k+1/2 .

Cohesive zone model
The cohesive zone model is used to simulate the debonding of the ITZ.The model was first developed by Dugdale [62] and later applied to concrete by Hillerborg [63].It is a phenomenological model in which the crack opening is characterised with a softening traction-separation law that starts once the maximum tensile stress of the material is reached.The model follows an exponential softening as shown in Fig. 4a and described in Eq. ( 36) [50].
where G f is the critical fracture energy, K pen is a large penalty stiffness that minimises the crack opening before reaching f t and ζ = max{λ(t) : t = 0, . . ., t k } is the maximum value of the equivalent crack opening λ in the loading history.The equivalent crack opening λ is defined as a function of the normal displacements u n and tangential displacements u t [64,65] The normal traction T n and tangential traction T t are expressed as a function of the normal crack opening u n and tangential crack opening u t respectively as Similar to the damage model, the traction-separation law also needs to be adjusted to account for aging, which can also influence the critical fracture energy G f .Due to the lack of research on the variation of G f for cementitious materials over time, in the present study, it is assumed to follow the change of the tensile strength over time following Eq.(22). Figure 4b illustrates the load path of the traction separation law that carries load from a previous time step.When loaded at different ages, the traction separation law is assumed to behave following the same concept as in the damage model in which we assume that the stress state and damage level do not change during aging.For ITZ, the damage level is determined as the reduction of secant stiffness with respect to the initial penalty stiffness.

Moisture diffusion
Problems involving moisture or temperature field are typically modelled using diffusion equations.In this study, the non-linear moisture diffusion equation proposed by Bažant and Najjar [53,66] is adopted to determine the moisture profile within a concrete meso-structure.The governing differential equation for moisture transport of cementitious materials is expressed as where C(h) is the moisture diffusivity that depends on the pore relative humidity and ∇ is the differential operator.The dependency of C(h) on humidity is expressed as where C 1 is the diffusivity when fully saturated, α 0 is the ratio of C(0)/C 1 (approximately), h c is the humidity at which ]/2, and r is a parameter that affects the shape of the diffusivity curve.This expression has been adopted in many past studies [18,46,67], and it has also been recommended by the fib model code [12].

Concrete meso-structure and quadtree mesh generation
In the present study, the concrete meso-structures used for the analysis are numerically generated.They are generated based on a take-and-place algorithm presented in Guo et al. [68].The aggregates are assumed to be in circular shape.The As highlighted in the introduction section, most of the existing meso-scale studies were carried out on meshes consisting of a uniform element size or elements with similar sizes (from Delaunay triangulartion).Such an approach limits the scale of the analysis as large meso-structures would be discretised into meshes containing a substantial number of elements.In this study, a quadtree decomposition algorithm is adopted for image-based mesh generation which is a powerful alternative for the meshing of composite domains such as concrete meso-structures in 2D [49].It is based on the recursive subdivision of cells into four smaller ones of equal size until only one material exists in one cell.The generated quadtree meshes contain hanging nodes between adjacent cells of different sizes which cannot be modelled directly by traditional FEM.Modifications are needed to enforce their compatibility in FEM which is typically achieved by further decomposing the cells next to the hanging nodes into smaller triangular and quadrilateral elements [69,70].However, the SBFEM serves as an efficient alternative to solve this problem.The quadtree cells can be directly modelled using SBFEM due to its ability to handle arbitrarily sided polygons.The adopted hierarchical quadtree meshes consist of elements of different sizes which reduces the elements needed for computation comparing to the algorithms adopted by existing meso-scale studies.In addition, the quadtree mesh contains only six unique patterns and all the quadtree cells in a mesh can be generated by scaling and rotating these six cells.This feature makes the analysis of quadtree meshes using SBFEM computationally more efficient, as the information of the six unique cells can be pre-computed and quickly retrieved during the analysis.The concept of quadtree mesh generation algorithm outlined in this section may also be applied for 3D problems.For more details of the quadtree algorithm and its coupling with SBFEM, see Saputra et al. [49].

Geometry representation
Unlike in FEM where the integral over an element is approximated using Gaussian quadrature, an SBFEM element is integrated analytically along the radial direction and approximated by polynomials in the circumferential direction.The modelling of an arbitrary n-sided polygonal element in SBFEM has been illustrated by Song [71] and it has also been presented in a previous work by the authors [50], so only a brief introduction is provided here.
For each polygonal element, a scaling centre is defined such that it is visible to every point on the boundary.The boundary of a polygon is discretised using line elements with shape function [N (η)] and natural coordinates −1 ≤ η ≤ 1.The surface of the polygonal element which consists of multiple triangular sectors is obtained by scaling the discretised boundary towards the scaling centre along the radial direction.A radial coordinate ξ is introduced, with ξ = 0 at the scaling centre and ξ = 1 at the boundary.The Cartesian coordinates of an arbitrary point (x, y) within a triangular sector are expressed as where T are the nodal coordinate vectors of an M-node line element, and ] is the corresponding one dimensional shape function.

Diffusion equation
The detailed derivation of SBFEM for linear diffusion was originally presented by Song and Wolf [72].In the present study, it is extended to solve for nonlinear diffusion problems.Since most of the formulations remain unchanged, only the essential steps are presented here for brevity.
The governing equation for nonlinear diffusion takes the following form in two-dimension: where {L h } is the differential operator and {q} is the flux vector related to humidity h as with the isotropic diffusion matrix [κ h ] that depends on the humidity

General solution
The humidity field h(ξ, η) at any point along the circumference may be interpolated as By applying the weighted residual method to the governing equation in the circumferential direction η, the SBFEM equation for diffusion in terms of h(ξ ) is expressed as where {} ,ξ and {} ,ξ ξ denote the first and second order derivatives of a variable with respect to ξ respectively.E 0 , E 1 , E 2 and M 0 in Eq. ( 48) are the coefficient matrices, expressed as: in which B 1 and B 2 are the SBFEM strain-displacement matrices.The steady state of Eq. ( 48) can be transformed into a set of first-order ordinary differential equations with twice the number of unknowns as where [Z ] is a Hamiltonian matrix and {Q(ξ )} is the internal nodal flux A block diagonal Schur decomposition of [Z ] is then performed, which results in where [S n ] and [S p ] are upper triangular matrices that contain eigenvalues with negative and positive real parts along the diagonal respectively.For bounded domains, only [S n ] leads to finite humidity at the scaling centre.The variables [ 11 ] and [ 12 ] are transformation matrices that correspond to the humidity, whereas the variables [ 21 ] and [ 22 ] are transformation matrices that correspond to the flux.The general solutions for humidity and flux vector are then derived as where the integration constants c n are determined by the nodal humidity on the boundary {h} = {h(ξ = 1)} as The steady state stiffness matrix [K ] can then be obtained by substituting Eq. ( 59) into Eq.( 58) and set ξ = 1, which is expressed as

Shape function
By substituting Eq. ( 59) into Eq.( 57), the humidity field along the radial direction is expressed as By substituting Eq. ( 61) in Eq. ( 47), the humidity field at any point within the element can be interpolated as where [N ] is the scaled boundary shape function, expressed as

Governing equation for diffusion
The diffusion equation of an SBFEM element takes the same form as the standard FEM, expressed as where [K ] is the stiffness matrix and [M] is the mass matrix with unit density.The formulation for stiffness matrix has been derived in Eq. ( 60).The mass matrix is formulated as where [m] satisfies the Lyapunov equation The time derivation in Eq. ( 64) shall be discretised in order to run incremental analysis.For time step t k = t k+1 − t k , the diffusion equation in Eq. ( 64) is approximated using finite differences and the backward Euler method as The backward Euler method is adopted here as it was found to be unconditionally stable.For clarity, Eq. ( 67) can then be rearranged to separate the humidity field at different time steps as

General solution
Similar to the humidity field, the displacement field u(ξ, η) at any point within an SBFEM element is interpolated as where [N u (η)] is the shape function for displacement.The SBFEM equation for displacement shares the same form as the steady state diffusion equation, expressed as For displacement problems, the diffusion matrix [κ h ] in Eqs.(49) to ( 51) is replaced with the elasticity matrix [D].The SBFEM equation for displacement in Eq. ( 70) can be solved following the same procedures presented in Eqs.(53) to (63).
The general solutions appeared in Eqs. ( 57) and ( 58) in this case are for the displacement {u(ξ )} and internal force vector {q(ξ )} respectively.The stiffness matrix formulation remains the same as presented in Eq. ( 60)

Shape function and B-matrix of SBFEM elements
For displacement problems, the strain-displacement relationship is governed by where [L u ] is the differential operator for displacement and its expression in terms of the scaled boundary coordinate system takes the following form.
where [b 1 ] and [b 2 ] are matrices that depend on the geometry of the elements only.By substituting Eqs. ( 69) and ( 72) into Eq.( 71), the strain field is expressed as where [B] is equivalent to the B-matrix in FEM, [B 1 ] and [B 2 ] are the SBFEM strain-displacement matrices.

Governing equation for time-dependent displacement
The incremental SBFEM formulation is derived from Eq. ( 18) by applying the principal of virtual work.It shares the same form as the one presented by the authors for the simulation of creep and shrinkage [50][51][52].For an undamaged SBFEM element, the incremental force-displacement relation at a given time step k with time increment of t k = t k+1 − t k is given as where [K ] k is the stiffness matrix as expressed in Eq. ( 60).In this case, the coefficient matrices in Eqs. ( 49) to (51) that are used to construct the stiffness matrix is formulated using the pseudo modulus E k (see Eq. ( 19)) instead of [κ].{ F ext } k is the external load vector, and { F } k is the equivalent load vector due to time-dependent strain components, expressed as From the above equation, [B] is the strain-displacement matrix, [D] k+1/2 is the material constitutive matrix evaluated at the middle of each time step.For two-dimensional problems, d is expressed as The equations presented above shares the same form as the ones presented by the authors for the simulation of creep and shrinkage [50][51][52], except that the formulation has been extended in the present study to explicitly account for drying creep.In addition, the variation of the prescribed strain { } k within an SBFEM element was approximated by a polynomial function in previous works.However, due to the complexity of the problem in this study, it is assumed that { } k is a constant within each element and it is evaluated at the scaling centre of each element.

SBFEM formulation for coupled time-dependent deformation and damage
Previous works conducted by the third author accounted for damage using SBFEM [73].Here, the same concepts are used and further extended to cover damage in structures that undergo time-dependent deformations.Due to the complexity of the problem, the damage index ω is assumed to be uniform within each element.Such an approach significantly simplifies the formulation and it is justified as long as sufficiently fine elements are used.By combining Eqs. ( 18) and ( 29), the total force-displacement relation for coupled isotropic damage and time-dependent deformation analysis is formulated as where the external load vector {F ext } k and effective equivalent load vector { F } k are obtained by summing their incremental values up to time step k The term 'effective' depicts variables formulated using undamaged material properties, and { F } is determined following Eq.( 76).To facilitate time-stepping analysis, Eq. ( 78) needs to be converted into an incremental form.For a small time step, the incremental force-displacement relation is linearised as where [K dmg ] k is the stiffness matrix used for the iterative procedures and the tangent stiffness is adopted typically.{ Fint } k−1 is the effective internal load vector at step k − 1 with expression of Newton-Raphson iterative procedures are used to solve the nonlinear problem.Due to complexity of the derivation of the tangent stiffness matrix and the computational cost associated with its formulation at every time step and iteration, the secant stiffness matrix is adopted for the iterative procedures in this study.At step k, the secant stiffness matrix for each element can be formulated as The secant stiffness matrix can be quickly computed in each iteration by directly scaling the pre-computed stiffness matrix.

Interface elements
Interface elements are used to accommodate the cohesive zone model which is used to simulate the debonding of the ITZ.A secant modulus approach is also used for the interface elements.Their secant stiffness matrix K z is assembled as where [D z ] k+1/2 is the material's secant stiffness evaluated at the middle of each time step following Sect.2.4, [R] is the rotation matrix that transforms the vectors from global coordinate system to local coordinate system and [B z ] is the shape function matrix.Their formulations can be found in Zhang et al. [50].The internal load vector of interface elements at any time step k is expressed as The Newton-Cotes integration scheme is adopted in the formulation of the interface elements because the standard Gaussian integration scheme leads to spurious oscillations when a high penalty stiffness is used [74].

Solution procedure
As mentioned in Section.3.1, a two-dimensional quadtree mesh only has six unique cell patterns.For each material that could appear in the meso-scale analysis, a set of six master cells is defined by considering the cells to have a unit edge length and a unit Young's modulus.The stiffness matrix [K ], and mass matrix [M] can be pre-computed for each master cell and each material type before the analysis.Other important matrices including B 1 , B 2 , [S n ] and [ 11 ] can also be pre-computed to speed up the formulation of the equivalent load vector { F } k .During the analysis, these matrices can be quickly retrieved and scaled to account for cells with different side length and material properties.
The present hygro-mechanical model carries out analysis on two layers including the hygral layer which determines the distribution of humidity within a structure and a mechanical layer which determines the deformation due to time-dependent behaviour.For clarity, the computational procedures involved in the two layers are presented separately.

Hygral layer
Box 1 shows the procedures for the analysis of moisture movement.Due to the dependence of the stiffness matrix on diffusivity C(h), iterations are carried out (see steps b to f).In each iteration, the stiffness matrix [K ] of each element is updated based on the current humidity at the scaling centre and the respective diffusivity C(h).A residual vector R h is defined by rearranging Eq. ( 68) as and the iterations stop once the norm of the residual vector |R| is smaller than a tolerance value R h,tol

Mechanical layer
Box 2 shows the procedures for analysis of the timedependent deformation.This is again a nonlinear problem that requires iterations due to material damage and ITZ debonding.Steps a to f are computed only once at the beginning of each time step.The iterations happen between steps h to k.The damage index ω k is updated every iteration from which the secant stiffness and residual load vector are determined.The residual load vector for SBFEM elements R u,e and interface elements R u,z are determined as Box 1 Computational procedure for hygral layer Similar to the hygral layer, the iterations stop once the norm of the residual vector |R u | is smaller than the tolerance value R u,tol .After convergence is attained, the internal variables γ μ,k are updated, which will be used in the subsequent time step to compute the incremental strain due to basic creep ve,k+1 .

Numerical examples
This section starts with a numerical example that simulates the experiment by Bažant and Xi [75].This example aims to demonstrate the capability of the proposed model in modelling drying creep with account of changes in the humidity distribution over time.This is followed by the simulation of the experiment by Idiart et al. [76] which demonstrates the capability of the proposed model in capturing shrinkage-induced cracks in concrete that is subjected to uniform drying.This section concludes with a numerical simulation of the differential drying of a concrete slab.The effect of boundary conditions and aggregate size on shrinkage-induced cracks are investigated through a parametric study.

Simulation of drying creep
Bažant and Xi [75] carried out a set of experiments to demonstrate that both drying creep and microcracking contribute to the apparent strain measured on drying creep specimens.The experimental setup of Bažant and Xi [75] is illustrated in Fig. 5 which includes testing of specimens under eccentric loads.The equivalent stress distributions of the eccentric loads have also been shown in the figure.The specimens tested were square prisms with length of 406.4 mm and side of 101.6 mm.Two load eccentricities were considered, a small one of 5.3 mm and a large one of 24.6 mm.For each set of tests, sealed and unsealed specimens were tested.The sealed specimens undergo basic creep only, whereas the unsealed specimens undergo basic creep, drying creep and drying shrinkage.For these specimens, drying was only allowed in one dimension that is on two of the opposite side surfaces.The specimens loaded with small eccentricity were under compression, whereas the specimens that were loaded with large eccentricity were subjected to tension on the left side that is less than the tensile capacity of concrete f t .Nevertheless, this tensile stress when combined with internal restrained shrinkage and differential drying can cause microcracking that can be simulated with the proposed meso-scale model, but cannot be captured using macro-scale homogeneous models.The environmental humidity was kept at 50% throughout the experiment.
The concrete used for the test was made from weight ratio of water: cement: sand: gravel = 0.5:1.0:2.5:3.0.The coarse aggregates used for the mix have a maximum size of 19 mm.Based on the mix proportion and the typical weight density of the materials, the proportion of mortar and coarse aggregates for the meso-scale analysis is estimated to be 61.1% and 38.9% respectively.The concrete meso-structure generated for the analysis is shown in Fig. 6 which is represented by an image with a resolution of 4096 by 1024 pixels.The quadtree decomposition criteria is set to give element sizes between 1 pixel and 4 2 pixels.The produced quadtree mesh contains 488,209 elements with a total number of 553,712 nodes.A small section of the quadtree mesh is presented in Fig. 6.The ratio of the number of elements generated to the total number of pixels is around 11.6%.An overall fine mesh is adopted in this study for the convergence and accuracy of the simulation because the crack pattern is not known a priori and it propagates or change with time due to the time-dependent deformations.It should be noted that the SBFEM approach allows for adaptive local mesh refinement without the need to regenerate the entire mesh [77], which can further reduce the computational cost.However, this feature of the SBFEM approach is not adopted in this study.
The material parameters adopted for the meso-scale analysis are shown in Table 1.For an ideal validation of the proposed model, all of the input parameters should be obtained or calibrated from experimental data on the same or similar mortar mix.However, due to the lack of experimental data on mortar, some of the parameters adopted in this example are back calibrated in a way that they give a good fit against the experimental results of the concrete specimens.In existing meso-scale studies, the elastic modulus of mortar is well accepted to be in the range of 20-30 GPa, whereas the aggregates have an elastic modulus of 60-70 GPa [29,35,46].The Poisson's ratio of both mortar and aggregates are 123 Box 2 Computational procedure for mechanical layer Eq. ( 7) Fig. 5 Experimental set up of Bažant and Xi [75] commonly assumed to be 0.2 [35,43].The aging coefficient s is adopted from Zhang et al. [51] which was calibrated for a similar mortar mix.The coefficients that relate moisture to hydration and creep are adopted from Bažant and Chern [14].
The basic creep properties of mortar are calibrated such that the simulated deflection agrees well with the experimental results.The validity of meso-scale models for the prediction of basic creep of concrete from mortar has already been demonstrated in a previous work by the authors [51] through comparison with accompanied experimental program.Both drying creep and drying shrinkage are related to the humidity which is obtained from the nonlinear diffusion analysis.For diffusion, the parameters α 0 , h c and r take the values from the fib model [12].C 1 is back calculated from that of the concrete which is estimated based on the fib model [12].A simple formula based on the proportion of mortar is adopted to determine C 1 for mortar assuming that aggregates are impermeable.For parameters relating drying shrinkage and drying creep, k sh is estimated based on the shrinkage results presented in [78].Due to the lack of relevant data, the parameter r governing drying creep is back calibrated to ensure a good agreement between the simulation and experimental results.For damage parameters, the tensile strength of mortar is assumed to be 10% the compressive strength of typical mortar mixes and the remaining damage parameters are back calibrated.The tensile strength and fracture energy of the ITZ are assumed to be f t = 1.5 MPa and G c = 0.01 N/mm following [46].It should be noted that due to the fact that some of the parameters are not adopted from relevant experiment data, the comparison between test and numerical results as shown subsequently is not a quantitative validation In the work of Bažant and Xi [75], eccentric loads were applied through pins at the two ends of the specimen with steel plates attached between the specimen and the pins for load transfer.Since the exact dimensions of the plates were not reported, we assume that the eccentric point load is transferred to the entire section of the specimen at the edges, which is simulated through a linearly distributed load on the top and bottom surfaces.The simulated results under sealed condition (basic creep only) along with the experimental results are shown in Fig. 7. Figure 7a shows the results under load with small eccentricity and Fig. 7b shows the results under load with large eccentricity.The calibrated model parameters for basic creep is able to replicate the test results under both loading scenarios as observed from the figure .For the specimens that are subjected to drying, all the timedependent strain components need to be considered in the simulation, including basic creep, drying creep and shrinkage.The last two are directly related to the internal humidity.Figure 8 shows the simulated change of humidity over time in the mortar phase.Mortar is assumed to be in fully humid condition at the start of drying and aggregates are assumed to be impermeable.After 1 day, it can be observed from the figure that drying only occurred close to the exposed surfaces and the majority of the mortar remains fully humid.Drying propagates to the interior of the specimen over time and the centre of the specimen reaches a humidity of around 75% after 200 days of simulation.Please notice that the humidity at the surface is equal to the external relative humidity of 50%.
The effect of damage caused by tension at the meso-scale was found to be insignificant for the specimen that was loaded with small eccentricity.Figure 9a shows the damage level of the mortar phase at 200 days after loading for the specimen that was loaded with large eccentricity.It can be observed that cracking is more significant on the left edge.This is due to the combined effects of the eccentric load and differential drying.Concrete is subjected to differential shrinkage due to the non-uniform distribution of humidity.Shrinkage is higher on the surface of the specimen than the interior, which introduces large tensile stress on the surface and cracking.When these tensile stresses are combined with eccentric loading, the cracking pattern becomes different at the two edges of the specimen.In the examined case, cracking is more severe at the left edge.Most of the observed cracks are surface cracks that are developed between the exposed surface and the outer layer of aggregates, and they are generally perpendicular to the drying surface.On the left edge, some of the cracks are also observed to penetrate further into the specimen.It can be seen that cracks generally propagate next to the aggregates which restrain the shrinkage deformation of mortar and introduce stress concentrations.Figure 9b shows the decohesion of the ITZ in the deformed meso-structure.Only minor ITZ decohesion is observed on the left side of the specimen.The above observations were not reported in the tests, and in general, they can hardly be measured.
Figure 10 shows the comparison between simulated and measured mid-span deflections when the specimens are subjected to drying.A good agreement is obtained for the specimen loaded under small eccentricity.For the specimen loaded with large eccentricity, good correlations are observed up to about 40 days.After that, the rate of the increase of deflection over time observed from the experimental results decreases which leads to some minor difference between the simulated and test results.Nevertheless, an overall good agreement is obtained.
It should be noted the presented study is in 2D which is not accurately representative of real concrete specimens.Concrete specimens and aggregates are generally 3D in nature which can affect the stress field and damage pattern.It has been demonstrated in previous works by the authors [51,52] that the use of 2D meso-scale models is sufficient for the prediction of creep and shrinkage strain in typical laboratory size concrete specimens that are loaded in one direction only as the out-of-plain restraint is generally small.This was shown through good comparison between simulated and measured responses.However, the accuracy of 2D simulations in predicting crack patterns and humidity distribution requires further validation.Nevertheless, the present model serves as a meaningful tool for clarifying the effect of various factors on the time-dependent response, and as a basis for the development of 3D models in the future.

Simulation of shrinkage-induced cracks under uniform drying
The experimental program of Idiart et al. [76] is used for this simulation.The experiment involves testing of 2 mm thin plate of cementitious composite subjected to drying.The plate is made of cement paste with cylindrical aggregates made from stainless steel.The cement paste was made with a water-to-cement ratio of 0.5 and the specimens were allowed to dry after 28 days of wet curing.The crack patterns were recorded at the end of the drying period and were reported in the study.The exact positions of the aggregates have also been given which allows a qualitative comparison with the experimental results.The idea of using thin samples in the tests was to simulate uniform drying scenario which in turn produces 2D crackpatterns.Under uniform drying, shrinkage induced cracks The environmental humidity started at 100% and was gradually decreased to 60% humidity over the course of 7 days and then kept constant for an additional 31 days.Due to the specimens being very thin, we can safely assume that the internal humidity of the specimen follows that of the environment.The gradual decrease of humidity was not specified in the study, and so it is assumed to take the shape of Fig. 11 with a drying rate of 1% per hour between the various humidity levels.
The experiment considered two aggregate distributions with 10% and 35% aggregate volume fraction respectively.In addition, four different sizes of the sample are considered The aggregate sizes were scaled along with the samples with diameters of 2 mm, 3 mm, 4 mm and 6 mm respectively, to ensure a uniform geometry on the meso-scale.The distribution of aggregates adopted by the experiment [76] are shown in Fig. 12.Digital images with a resolution of 2048 by 2048 pixels are generated to represent the adopted mesostructures.The quadtree decomposition criteria is again set to give element sizes between 1 and 4 2 pixels.The produced quadtree mesh contains 291094 elements with a total number of 300574 nodes for the meso-structure with 10% aggregate volume fraction, and 363880 elements with 394694 nodes for the meso-structure with 35% aggregate volume fraction.The ratio of the number of elements generated to the total number of pixels is 6.9% and 8.7% for the two meso-structures respectively.Similar to the previous example, a fine mesh is adopted here to ensure the damage zone is accurately captured.
The crack patterns were reported for meso-structures with 2 mm and 10% aggregates, 6 mm and 10% aggregates, and 4 mm and 35% aggregates.Despite that the specimens were thin and the humidity was carefully controlled, a uniform crack pattern across the thickness was not obtained in the tests.Two types of crack patterns were reported, including middle-section cracks and through-going cracks.In the present study, because of the 2D characteristics of the model, the simulated results are compared with the measured cracks at middle-section.
The material properties adopted in this example are summarised in Table 2.The elastic properties of cement paste and steel aggregates are adopted from Idiart et al. [76].The parameters related to the time-dependent deformation are calibrated from experimental data reported in Zhang et al. [51] for a similar type of mix.Due to the lack of relevant data, the damage parameters are calibrated such that the computed cracks patterns agree reasonably well with the measured cracks.It is worth noting that the crack path is mainly governed by the elastic properties of the steel aggregates and the time-dependent properties of cement paste, whereas the damage parameters mainly affect the length of the developed cracks.It was found that the predicted response will only slightly change with varying the input damage parameters within the logical range, which provides a level of validation to the proposed model.Figure 13 shows the comparison between computed and measured crack patterns.In the simulated results, crack patterns are presented using damage levels and they are manifested as a band of damaged region.It can be seen that the simulated crack patterns agree reasonably well with the measured crack patterns.It should be noted that an exact prediction of the crack patterns is very difficult considering that cement paste is a highly heterogeneous material.

Simulation of shrinkage-induced cracks under differential drying
In this example, a concrete specimen subjected to differential drying is considered in order to simulate real scenarios in construction where only the upper surface is typically exposed to the environment.Under this condition, shrinkage-induced cracks develop due to the combined effect of aggregate restraint and differential drying.A parametric study is carried out to investigate the effect of various parameters on shrinkage-induced cracking.A component of a concrete slab with thickness of 100 mm by 200 mm length is considered.
For brevity, we will consider the same concrete mix and material properties as adopted in Sect.5.1.However, typical The effect of ambient humidity is investigated first.The generated digital image of concrete meso-structure is shown in Fig. 14 with a resolution of 1024 by 2048 pixels.The quadtree decomposition criteria is set to be the same as previous examples.Three constant environmental humidity levels are considered, including 30%, 50% and 70%.The slab is supported on the bottom surface against vertical movement and only one point is restrained in the horizontal direction.The simulated crack patterns under the three different cases are shown in Fig. 15.Only the top 30 mm of the slabs are presented in the figure as no damage is observed in the inner part of the slab.It can be clearly observed from the figure that a lower environmental humidity leads to more pronounced cracking.This can be attributed to the higher tensile stress induced on the drying surface as a result of the steeper humidity gradient between the exposed surface and the interior of the slab.It can be seen that for the case with 30% humidity, a more severe crack pattern is observed in which the cracks penetrated deeper inside the slab and cracks between aggregates start to develop.
In the second part, the effect of aggregate size on shrinkage-induced cracking is investigated.Three concrete  16 shows the simulated crack pattern on the top 30 mm of the slabs with different aggregate sizes.More cracks are observed in the concrete with smaller aggregate.This can be attributed to the increase of the quantity of aggregates which created more locations with stress concentrations for cracks to develop.On the other hand, concrete with larger aggregates is observed to have wider damaged zone or in other words larger cracks.A similar observation has also been made in the numerical study by Grassl et al. [80].
Another simulation is carried out here to investigate the effect of horizontal restraint on shrinkage-induced cracks.In addition to the restraint on the bottom surface, the slab is restrained from horizontal movement on the left and right surface in this case.The slab is subjected to an environmental humidity of 50%.The simulated crack patterns in this case are shown in Fig. 17.When the slab is restrained from horizontal movements, one localised cracking path is observed in addition to the typical surface cracks observed in Fig. 15.This is due to the tensile stress introduced in the horizontal direction as a result of the structural restraint.In this case, the debonding of the ITZ is also clearly visible especially around the localised crack as shown in Fig. 18.As already highlighted in the first numerical example, the observations obtained from this example are qualitative ones due to the 2D nature of the proposed model.

Conclusions
An SBFEM based 2D meso-scale hygro-mechanical model has been presented in this paper for the modelling of the time-dependent behaviour of concrete.The model takes digital images of concrete meso-structures as input and converts them into meshes through a quadtree decomposition algorithm that is compatible with SBFEM analysis.Concrete is treated as a two-phase composite consisting of aggregates and mortar, and the ITZ between mortar and aggregates are modelled using the cohesive zone model.
The model accounts for basic creep, drying creep, drying shrinkage and cracking.Basic creep is treated as viscoelastic, and it is modelled using a rate-type Kelvin chain model.Drying creep is modelled using a viscous unit which depends on the current stress level.Both drying creep and drying shrinkage are related to the internal humidity level.The internal humidity distribution within the concrete meso-structure is obtained from a diffusion analysis in which moisture movement in mortar is governed by a nonlinear diffusion equation and aggregates are assumed to be impermeable.Cracking of concrete on the meso-scale is modelled through the combination of a continuum damage model for mortar, and a cohesive zone model for the ITZ.
The capability of the proposed model in predicting the time-dependent deformations and cracking patterns was verified through comparison with experimental observations from the literature.A parametric study was also carried out which demonstrated the effect of ambient humidity, aggregate size and restraining condition on the development of shrinkage-induced cracks in a concrete specimen subjected to differential drying.It was found that an increase of aggregate size reduces the number of surface cracks observed, but it generates larger cracks.Given the capability of the proposed model, it serves as an effective tool for the simulation and understanding of the time-dependent behaviour and cracking of concrete under various realistic conditions, i.e. variable loading and variable drying condition.The model may also serve as a powerful tool to examine the role of various components of concrete and their properties on the time-dependent behaviour of concrete, which can assist in the optimisation of concrete mix for serviceability.In future works, the proposed meso-scale hygro-mechanical mode will be extended to simulate the time-dependent behaviour on 3D concrete specimens in which a more accurate representation of the humidity and damage fields can be obtained.The model will also be coupled with an adaptive mesh refinement technique to reduce the computational cost.

Funding Open Access funding enabled and organized by CAUL and its Member Institutions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 1
Fig. 1 Illustration of the constitutive model

Fig. 9 Fig. 10
Fig. 9 Fracture of concrete under large eccentric load (t = 200 days) Fig. 10 Mid-span deflection over time under drying condition

Fig. 18
Fig. 18 Simulated ITZ debonding under 50% humidity and restrained horizontal movement 1. Initialise the humidity {h} 1 and diffusivity C(h 1 ) Eq. (41) 2. For time step k = 1 : n with time increment of t k = t k+1 − t k a. Construct the stiffness matrix [K ] and the mass matrix [M] If |R h |> R h,tol , return to step b, else, go to the next time step 1. Initialise the internal variables γ μ,1 and damage index ω 1 2. For time step k = 1 : n with time increment of t k = t k+1 − t k a. Determine the incremental hydration time t e,k and the incremental reduced time t r ,k Eqs.(24, 25) b.Compute the elastic modulus E k+1/2 and the tensile strength f t,k+1/2 Eqs.(21, 22) c.Determine the incremental pseudo modulus for basic creep E ve,k and drying creep E dc,k Eqs.(5, 16) d.Determine the incremental strain due to basic creep ve,k , shrinkage sh,k and drying creep

Table 1
Adopted parameters for mortar and aggregates in Example 1

Table 2
Adopted parameters for cement paste and steel aggregates -Idiart's experiment