Abstract

To accurately predict the time-dependent deformation of concrete, a multiscale model with its focus pinned on mesoscale is proposed here to break down the constitutive law of concrete to the mechanics of its different constituent phases. A three-phase unit cell, consisting of one coarse aggregate, mortar matrix, and the interfacial transition zone (ITZ), is employed to represent the basic structural element of concrete on mesoscale. Following Eshelby’s inclusion theory, the Mori-Tanaka homogenization, continuous retardation spectrum method, and isotropic continuum damage model are applied to capture the time-dependent behavior of the unit cell. To take into account the shape effect of aggregate, the explicit Eshelby’s tensor of polygonal inclusion is obtained based on an enhanced approach. The proposed multiscale material model is incorporated into ABAQUS, and its effectiveness and robustness are documented by the simulations of unit cells containing aggregates of different shapes.

1. Introduction

A realistic prediction and statistic evaluation of the time-dependent behavior of concrete attract enormous attention of researchers and engineers in civil engineering because it represents an essential aspect of integrity and serviceability of critical concrete structures expected to provide high-quality service for a long lifespan. Many critical concrete structures are very sensitive to their long-term deformation. If this long-term deformation is primarily caused by concrete creep, we call this type of concrete structure as a creep-sensitive structure [1, 2]. Underestimating a creep-sensitive concrete structure’s long-term deformation generates a profound impact on its structural integrity and serviceability. For example, the excessive deflection will dramatically shorten the service lifetime of a prestressed box girder, and uneven deformation will lead to unexpected stress redistribution triggering severe cracking in a concrete member [3, 4]. Therefore, establishing an accurate and efficient numerical framework for creep structural analysis, which of course must be pinned on realistic constitutive relations, is imperative for the current progress towards durable infrastructures symbolized by a significantly longer lifespan.

Concrete creep, described as the strain growth under sustained stress, originates from its unique microstructure exhibiting both capillary porosity and subcapillary nanoporosity [5]. Therefore, it is not surprising that the phenomenological creep models based on the macro-scale observations, for example, ACI, CEB-FIP (2010), fib (2010), B3, and B4 models, cannot capture the essential time-dependent mechanics of concrete creep from the mesoscale or nanoscale. Recently, the bottom-up approach built upon multiscale and multiphase modeling has become increasingly popular. Its applicability to concrete creep enables us to capture the local response within concrete’s highly heterogeneous microstructure. It is widely accepted now that realistic modeling of the time-dependent mechanics of the different constituent phases and the interactions among them is the key to accurately predicting the multidecade creep concrete [68] (Maekawa 2012).

For multiscale modeling of cementitious materials like concrete, there are currently three popular approaches: (a) multiple-phase representations [6]; (b) micromechanical models based on matrix-inclusion morphology [9, 10]; and (c) fine discretization of concrete microstructure by 2D or 3D elements [11, 12]. The first two approaches usually idealize and simplify the aggregate shape into circle or sphere ones in a representative volume element (RVE) of concrete. As a consequence, analytical or semianalytical predictions for the concrete’s behaviors can be derived. The third approach models the aggregate shape in an explicit way within the framework of finite element (FE) solver. Although with improved accuracy, the computational cost is expensive and renders its wide application in large-scale concrete structures. Despite their differences in the modeling strategies, all three approaches are capable of capturing the concrete’s local responses based on the different constitutive relations of the different mesoscale material phases.

In this study, the characteristic size of the proposed multiscale numerical framework to represent the concrete microstructure resides on the mesostructural modeling. The rationale to select mesoscale (mm) over grain scale (μm) and nanoscale (nm) is to achieve the balance between accuracy and computational cost, especially in the context of full-scale multidecade structural analysis.

In the mesostructural modeling of this study, concrete is represented by a 3-phase unit cell consisting of a single coarse aggregate, mortar matrix, and the interfacial transition zone (ITZ) between the first two. As demonstrated in Figure 1(a), the 3-phase unit cell is analogous to an Eshelby inclusion problem, which is widely studied in micromechanics [8, 13, 14]. The micromechanical solutions based on Eshelby’s tensor can be employed to project the macroscopic strains on the mesostructure (unit cell) after the necessary treatment of ITZ. The constitutive laws corresponding to the time-dependent mechanics of different phases will then be used to obtain the local mechanical response within the unit cell. Then, the average responses of the whole unit cell are determined based on a Mori-Tanaka type homogenization scheme [15].

In addition, to approximate better the real morphology of concrete mesostructure, the classic inclusion problem will be extended to take into account inclusions of shapes different from the ellipsoid. In this study, an enhanced approach based on Rodin’s formulation [16] is proposed to efficiently identify the exterior and interior point tensors when the inclusions have different polygonal shapes ranging from triangle to hexagon. The average of the position-dependent interior Eshelby’s tensor needed in approximation can then be identified based on its special property [17] or numerical integral [18].

In this study, the time-dependent deformation is focused on creep in mortar and cracking in ITZ. As for other nonnegligible inelastic phenomena, such as shrinkage, cracking in mortar, and cyclic creep, they will be studied under a similar framework in future investigations. To approximate the creep in the mortar matrix, a rate-type model enriched by the continuous retardation spectrum method and the exponential algorithm [19] will be employed. As for the weakest phase in the unit cell, the ITZ will follow a simplified plastic-damage model where the microcracking is captured by the isotropic degradation of its stiffness.

2. Mesoscale Concrete Modeling Based on a Three-Phase RVE

In multiscale modeling, the selection of representative elements depends on the characteristic size scale where the physical mechanism of interest can be quantitatively identified. While the fundamental mechanisms of creep are encoded in the porous microstructure of concrete, the basic laws of the time-dependent mechanics of concrete can be captured and quantified on mesoscale or grain scale [17]. Therefore, characteristic sizes ranging from grain scale (μm) to aggregate scale (mm) are all acceptable for mesostructural modeling of concrete creep.

Unfortunately, modeling on grain scale demands an extremely high computational cost [21], which is almost unaffordable given the size and time scales of interest in structural analysis. Therefore, mesostructural modeling based on aggregate scale is a rational selection to ensure the desired accuracy at an affordable computational cost. On this size scale, concrete demonstrates a highly heterogeneous mesostructure, which is characterized by randomly distributed coarse aggregates and surrounding mortar matrix.

This mesostructure can be represented by a unit cell consisting of one coarse aggregate and bulk mortar matrix, a composite of fine aggregates (sands) and hardened cement paste. Here, the assumption that mortar is a single homogenous phase is acceptable, even though its microstructure is heterogeneous. The reason is that the scale of micro-heterogenity of mortar is small, compared to the size of the aggregarate, also the meso-scale of concrete. Thus, only its averaged properties and responses need to be considered in mesostructural modeling.

Between the aggregate and mortar, there exists an interfacial layer called the interfacial transitional zone (ITZ). This thin layer, induced by the wall effect during cement hydration, exhibits different properties from the bulk cement paste. Experimental investigation shows that ITZ is thin (∼50 µm), porous, and of a higher CH (calcium hydroxide) concentration. Therefore, ITZ will be weaker than the bulk cement paste. Smadi and Slate [22] found that substantial microcracking was initiated in the ITZ when the macroscopic stress reached only about 60% of the concrete compressive strength. Therefore, the possible cracking in ITZ must be taken into account when the stress redistribution triggered by the long-term deformation is under consideration.

In view of this, the ITZ should be included in the unit cell for mesostructural modeling. A unit cell containing one aggregate, mortar matrix, and ITZ is plotted in Figure 1(a), in which the shape of coarse aggregate is idealized to a circle (sphere for 3D). The different properties of concrete will now be approximated directly based on the corresponding constitutive relations in different material phases to bypass the phenomenological formulation.

To accurately determine the local response of different phases within the unit cell, a straightforward strategy is to discretize the different phases by finer elements in FEM. However, this method, equivalent to using a size smaller than the aggregate scale, will dramatically increase the computational cost and thus impose an enormous computational burden on long-term analysis for large structures. Furthermore, the thin layer of ITZ also raises a challenge in mesh generation and geometry modeling.

An alternative to bypass these obstacles is to take advantage of Eshelby’s formalism. As shown in Figure 1(a), the configuration of the proposed unit cell is analogous to the classic inclusion problem, which means the micromechanical solutions may be applicable. As demonstrated in Figure 1(b), multiple sample points (usually 12 are enough for matrix) are selected in the unit cell. Then, the local stress or strain tensor represented by the sample point can be calculated based on Eshelby’s tensor, which is unique for a given unit cell.

3. Eshelby’s Equivalent Eigen-Strain Method

In the classic Eshelby’s inclusion problem, namely an ellipsoidal inclusion in an infinite matrix domain, the distribution of elastic fields in the unit cell can be effectively evaluated via the equivalent eigenstrain method [13, 14]. As demonstrated in Figure 2, this elegant method treats the inhomogeneity caused by the inclusion in the matrix as an equivalent eigenstrain based on the fact that in an ellipsoidal inclusion, there exists a uniform eigenstrain, which induces uniform disturbance. The mechanism in this inclusion problem can be described aswhere C and represent the fourth-order elasticity tensors of matrix and aggregate, respectively; , , and are the far-field strain, strain disturbance, and eigenstrain, respectively; Sin is the interior Eshelby’s tensor. The equivalent eigenstrain can be obtained by solving the equation system of equations (1) and (2). After that, the stress and strain distribution in the matrix can be identified with the help of an exterior Eshelby’s tensor.where Sex is the exterior Eshelby’s tensor.

In plain strain problem, for an ellipsoidal inclusion, the position-independent interior Eshelby’s tensor can be expressed aswhere is the Poisson’s ratio of the matrix. Unlike the interior Eshelby’s tensor, the exterior one is position-dependent. As shown above, by using the classic micromechanical solutions, the strain and stress tensors can be efficiently calculated based on the given far-field strain of this unit cell.

If there did not exist the thin layer of ITZ, the configuration of the unit cell would be the same as a classic 2-phase inclusion problem for which the classic Eshelby’s tensors can be directly applied. Although some 3-phase models have been developed to extend Eshelby formalism for considering an interfacial layer [23, 24], the constraints imposed in their formulations limit their application for general cases. Therefore, in this study, a Mori-Tanaka homogenization scheme is adopted to downgrade the 3-phase configuration to the classic 2-phase inclusion problem by smearing the response of ITZ on aggregate. It has been shown that this method may be ineffective when the inclusion is very stiff [25]. Fortunately, this is not the case for concrete because the moduli of its 3 phases are generally in the same order of magnitude.

4. Micromechanical Solutions of Polygonal Inclusions

In practice, the shape of aggregate used in mix design is far from circle and sphere. Indeed, inclusion displaying a polygonal shape is a more realistic representation of concrete mesostructure. Thus, investigation on the micromechanical solutions of the polygonal inclusions is of significant importance for modeling concrete mesostructure. As shown in Figure 3, unit cells of circular, hexagonal, squared, and triangular inclusion inside should be investigated and then implemented in the proposed mesostructural modeling.

To solve the polygonal inclusion problem, the key is to seek the evaluation of related Eshelby’s tensor, which maps the eigenstrain to disturbance. Following the universal expression provided by Mura [14] for Eshelby’s tensor at point x for uniform eigenstrain as shown in equation (2), Rodin [16] proposed an algorithmic closed-form solution for Eshelby’s tensor of a polygonal inclusion.

In equation (2), the harmonic potential and biharmonic potential of the inclusion are defined as

In Rodin’s approach for plain strain problem, the integral over polygonal inclusion is evaluated piece-wisely by dividing the p-sided polygon with arbitrary shape into a series of triangles sharing the same vertex x. In this partition, the base of each triangle is an edge of the polygon (see Figure 4).

In each triangle, the local coordinates and their basis pair are defined as shown in Figure 5. Accordingly, the positions of the vertices at the base are represented by the pairs and , which can be evaluated by vector production as shown in equation (4).where stand for the position vectors of the vertices at the base. When evaluating the integral over the whole triangle, it is much easier to further partition the duplex into two simplexes [16]. Then, the harmonic and biharmonic potentials for each triangle can be calculated via integral in local coordinates, which will provide the following results (see equations (8)–(10)).

The potentials obtained above will be differentiated with respect to x to get the contribution of each triangle to Eshelby’s tensor. Eshelby’s tensor for a polygon is obtained based on the assembly of the triangles.

It is worth to mention here that although Rodin’s algorithm is conceptually straightforward, it is not convenient, sometimes tedious in computation, which results from the repetitive coding of different triangles. If examining the partition carefully, we are able to find that the common geometry shared by the triangles can be utilized to simplify the calculation. Based on this observation, an improved algorithm is proposed.

As shown in Figure 5, for the ith triangle, by denoting as the angle rotating from x1 axis to the local axis, Ti (the transformation matrix for this triangle) can be expressed as follows:

It is not difficult to see and now in the transformed Cartesian system , and this will make the calculation of vector production in equation (4) much easier since the pairs are coordinate-independent. Besides, there exist the following relations:

For the ith triangle, the criteria used to judge whether the point of interest x falls into the inclusion can be now easily expressed as an inequality.

With the help of transformation, for each triangle, only the unique rotation angle and the coordinates of the two vertices at the base need to be prepared for the calculation. Therefore, it enables us to use the same program for every partition.

No explicit expression of the global Eshelby’s tensor for a polygon is given in Rodin’s paper, which might be attributed to the difficulty of carrying out the extremely complicated differentiation of potentials with respect to the global coordinates. To the best knowledge of the authors, there exists no such explicit expression based on Rodin’s evaluation of potentials. Therefore, they are addressed here for the purpose of coding convenience. Note that equations (14)–(19) should be combined with equation (2) to provide the explicit Eshelby’s tensor.where

From the obtained tensor, it can be seen that if the inclusion is not ellipsoidal; for example, polygonal, the algorithm based on Eshelby’s equivalent eigenstrain, is no longer effective. The interior Eshelby’s tensor is position-dependent; therewith, the uniform eigenstrain is transformed into a nonuniform one. Therefore, it eliminates the possibility to solve equations (1) and (2) simultaneously.

To overcome this obstacle preventing the use of the equivalent eigenstrain method, Nozaki and Taya proposed an approximate approach that is proved effective for polygonal inclusion with sides more than square [26]. This approximate approach assumes equation (20) can replace equation (2) when solving the equivalent eigenstrain.where ()ave stands for the average over inclusion. Moreover, there exists a special property of Eshelby’s tensor for regular polygons except for square [17]. It is found that the average of the interior Eshelby’s tensor is equal to the one at the center of this polygon. Therefore, in regular polygonal inclusion problem, except for the square case, equation (20) can be replaced with the following equation:

As for the square inclusion, the numerical approach has to be employed to find the average of Eshelby’s tensor. Special treatment of numerical integral must be employed in this case because of the local singularity. After obtaining the approximate eigenstrain, the subsequent steps will be the same as the classic equivalent eigenstrain method. By taking advantage of this enhancement, the equivalent eigenstrain method can be applied to the polygonal inclusions.

5. Time-Dependent Deformation in Mortar Matrix

In the unit cell, the mortar matrix accounts for the major portion of the hydration products of concrete. In view of this, it will be a realistic approximation that creep happens exclusively in the mortar phase. The reason is that creep has been found to occur predominantly in the cement gel of hardened cement paste [17, 27, 28].

The creep in the mortar phase is represented by the sample points inside the matrix. Their strains will change with time due to the creep. Generally, they can be expressed as εi(t) = J(t, t')σi, where J(t, t') is the compliance function, which represents the strain at time t induced by a unit sustained stress applied at age t'; εi is the strain, and σi is the sustained stress at the ith sample point. For service load, the concrete creep generally follows an aging linear viscoelastic stress-strain relation, which means the use of the principle of superposition is effective. However, the integral-type formulation resulting from the direct use of the principle of superposition brings two main disadvantages in large structural analysis: (1) computational inefficiency, and (2) incompatibility with other time-independent phenomena that do not obey the principle of superposition. To overcome this obstacle, it must be replaced with an equivalent rate-type law [2].

Recently, Yu et al. [2] proposed an efficient rate-type algorithm for full-scale creep structural analysis, in which the viscoelastic stress-strain relation is approximated by the Kelvin chain model consisting of a series of Kelvin units [29, 30]. The continuous retardation spectrum method [31] equipped by the Widder's approximate inversion formula based on Laplace transform inversion [3234] is used to determine the spectrum of Kelvin chain units uniquely. Unlike the conventional numerical fitting, the spectrum identified by this approach will satisfy the thermodynamic constraints imposed by the concrete aging. After that, an unconditionally stable algorithm, called exponential algorithm [35, 33], is employed to convert the first-order ordinary linear differential equations expressed by the Kelvin chain to quasi-elastic incremental stress-strain relation.

6. Damage Initiation and Accumulation in the ITZ

It is not difficult to imagine that significant stress redistribution will be generated because of the creep in the mortar matrix. As the weakest material phase in the unit cell, the ITZ is especially vulnerable to the stress concentration induced by the material mismatch. In order to capture the possible material degradation in the ITZ, a realistic constitutive law should be introduced for it.

Considering the thickness of ITZ, the contact surface model and cohesive crack model have been frequently used to approximate the behavior of ITZ recently. In this multiscale model, a simplified Mazars type damage model [37] is adopted, in which a damage indicator D is applied to describe the isotropic degradation of the material. Therefore, the constitutive law in ITZ can be written as , where is the initial elastic modulus of ITZ, which can be set as about 70% of mortar modulus based on the nanoindentation investigation. The isotropic damage indicator D is defined aswhere εcr stands for the critical equivalent strain of damage initiation; and stands for the largest value of the equivalent strain that has ever been achieved in the deformation history up to the current step. In this study, the equivalent strain in ITZ is determined by the positive principle strains as follows:

Here, the Macaulay operator, , is introduced to ensure the assumption that compression has no contribution to the damage initiation and accumulation in ITZ.

The tangential matrix Et, which is needed in a rate-type algorithm, can then be easily obtained according to the relation Et = /. In light of the possibility that the loading history in the ITZ may not be monotonic due to the stress redistribution, the loading and unloading cases are considered in this damage model. Here, Et is set equal to E0 when unloading and reloading happen.

7. Homogenization Scheme and Its Implement in FEM

In the proposed multiscale model, the macroscopic response of the unit cell should be obtained to carry on the calculation for the next step, which necessitates the use of the Mori-Tanaka type homogenization scheme for estimating the average property of the composite. In general, the Mori-Tanaka scheme is carried out in the form of an effective field, which results from Eshelby’s solution for the inclusion problem. The essential part of the homogenization scheme is to obtain the effective stress and strain field of the unit cell.

The effective stress and strain field is evaluated as

Within the framework of elasticity, the relationship between stress and strain in the r-th phase and the whole domain is defined as

Finally, the average property will be obtained based on the volume fraction and elastic tensor of each phase. One thing that needs to be noted is that the Mori-Tanaka homogenization will be employed twice due to the need of downgrading the 3-phase unit cell to the classic 2-phase inclusion problem.

8. Numerical Examples of Unit Cell Containing Polygonal Inclusion

In this study, the proposed unit cell model is implemented in the material subroutine UMAT of ABAQUS for FEM analysis. A proper arrangement of the exponential algorithm, micromechanical formalism, and Mori-Tanaka homogenization scheme is involved in this implementation. In the numerical simulations, the following material properties of mortar, aggregate, and ITZ are assigned: EM = 30 GPa, = 0.18, EA = 33 GPa, = 0.2, EITZ0 = 21 GPa, = 0.18, and εcr = 0.0005.

To verify this multiscale model, as shown in Figure 6, the simulation is carried out for the creep of unit cells equipped with different shaped inclusions. Regular polygons inscribed in the same circle as the circular inclusion are considered. The results of FEM analysis in ABAQUS based on the same, but fine-meshed, unit cells are set as a benchmark (see Figure 6(b)).

In Figure 7, it can be seen that for circular, hexagonal, and squared inclusions, the results of the deformation history of the unit cell under sustained loading agree well with the FEM results, only an error greater than 5% happening for triangular inclusion. This observation confirms Nozaki and Taya’s conclusion [26] that the approximate method provides good results for the polygonal inclusions with sides more than a square. Although the result for triangular inclusion is not as good as other polygons, the maximum error of 8.97% is still acceptable, and the multiscale result is conservative.

After the model is verified by comparing it to a fine-meshed model, the shape effect of the inclusion on the creep analysis is studied. A series of simulations are carried out for a concrete specimen under sustained loading with the size shown in Figure 8(a). Inclusions of different shapes, with a constant volume fraction of 7.1%, are considered in the model.

As shown in Figure 8(b), for hexagonal and circular inclusions, their results are very close to each other. This is reasonable because a hexagon is assumed a good approximation of a circle in terms of shape. However, the squared and triangular ones behave very differently from the circle and hexagon. It is worth to mention here that when simulating the time-dependent behavior of concrete, ignoring the shape effect of inclusions might cause errors greater than 10%, as shown in Figure 8(b).

9. Conclusions

Considering the shape effects of aggregates, an improved multiscale model based on the mesostructure of concrete is proposed in this study. Damage accumulated in the ITZ is also included in the calculation. To capture the aging viscoelasticity of mortar, rate-type formulation based on an exponential algorithm, coupled with continuous retardation spectrum, provides a more efficient approach than the primitive integral-type formulation. Furthermore, the high computation cost imposed by the structure size and time duration in long-term deformation analysis can be ameliorated by utilizing the micromechanical solutions based on Eshelby’s equivalent eigenstrain method. Numerical examples show the efficiency and effectiveness of this model in simulating the nonlinear time-dependent behavior of concrete. To accurately estimate the long-term deformation of concrete, the shape effect induced by aggregate must be considered. The present study shows that a 10% error may appear if the aggregate shape is ignored in the mesostructural modeling.

Data Availability

All data, models, or code generated or used during the study are available from the corresponding author by request.

Conflicts of Interest

The authors declare that they have no known conflicts of interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation for Young Scientists of China (51808113) and the National Natural Science Foundation for Young Scientists of Jiangsu Province (BK20180389).