Numerical modeling of ASR: an updated semi-empirical model

: Alkali-silica reaction (ASR) is a deleterious expansion phenomenon affecting concrete structures worldwide. It occurs when a susceptible chemical composition of concrete components is present, and it is particularly rampant in inherently humid environments. This phenomenon is exacerbated when found in mass concrete structures such as dams and foundations. The amount of volumetric ASR strain used to be deemed as nearly constant, however, recent advances have shown that it is actually affected by the distribution of volumetric stresses. Therefore, this behavior demands an update in the numerical models that have been devised to simulate the anisotropic ASR-driven expansion. This paper deals with Saouma & Perotti’s thermo-chemo-mechanical coupled model, which has been applied in a solely mechanical manner, and updated to account for a varying volumetric strain. A simulation of the experiment that shed new light on the variation of ASR volumetric strain was then carried out with the finite element method package COMSOL. As a result, significantly smaller errors in predicted strains were attained by the updated model in comparison to the original one, and consequently, the new model poses a promising tool for a more accurate simulation of ASR expansion.


The Mechanical Modeling Approach
The Colorado model entails an ASR constitutive model that relies on the assumption that ASR volumetric strain is virtually constant up to the failure region of concrete stresses while also tackling the highly anisotropic stress-strain issue by means of a system of weights to account for the distinct behavior of directional expansion.
In 2017, Liaudat et al. [23] , the latter a creator of the Colorado model himself, published an article entitled 'ASR expansions in concrete under triaxial confinement' [23] on experiments for the most part conducted at the 'alkali-aggregate reaction triaxial machine' built by Saouma ( Figure 1). This apparatus was employed to exert compression through confinement onto a 150x150x150mm approx. concrete cube, by means of 3 active (and 3 reactive) steel plates and rods connected to hydraulic actuators. The machine could deliver axial confining stresses up to -9 MPa while maintaining the specimen humid and at a steady temperature in the 30-70 °C range. Three axial load cases were studied: [-1 -1 -1] with = = =-1 MPa, [-9 -9 -1] with = =-9 MPa, =-1 MPa, and [-9 -9 -9] with = = =-9 MPa. The results of the experiment were analyzed and graphed in that paper [23]. In consequence, a new understanding of the relation between volumetric stress and volumetric strain in the context of ASR was attained, in essence: • Volumetric ASR strain can be altogether restrained by the application of stresses equal to or more intense than the critical ones along the three axes • A biaxial stress state produces ASR strain oriented towards the free direction that seems to be 50% of the total free expansion volumetric strain, when those stresses are both at the critical stress levels or more intense. In such cases, there is no expansion along the (compressed) loaded axes beyond the critical level The above postulations represent a departure from Saouma and Perotti's numerical model implementation of a nearly constant ASR volumetric strain, which allows for the possibility of ASR axial strains occurring even in a direction where the applied compressive stress is above the critical threshold level [1]. Therefore, the conclusion of the article [23] pointed out the need for corrections to the model.
The following relations have been proposed by [23] based on the previous work of Saouma and Perotti [15], adjusting the observed behavior of ASR volumetric expansion rate έ (1/h) under compressive stresses: where έ in Equation 1 is the observed volumetric expansion rate, έ , is the rate of volumetric ASR expansion without applied stresses, in Equation 2 represents the volumetric stress to which the concrete element is being subjected, i.e., = ( )/3 where are the principal directions stresses, � represents the stress under which expansion is totally suppressed, also referred to as critical stress in literature (In all the above, a convention is adopted that negative stress is a compression, positive stress is tension) The above function is a macro mechanical reduction of complex phenomena taking place that include cracking, creep, and chemo-transport [23]. Despite the simplification, the authors deem it accurate in fitting some empirical data.
Regarding critical stress, its value is affected by the concrete mix, as well as environmental factors, such as temperature and humidity [1], [20], [24], [25], usual values are in the range from 7 to 11 MPa. Next in Figure 2, considering a typical 10 MPa critical stress, a curve has been plotted to illustrate Liaudat et al [23] proposed function:   According to figure 3, the total ASR volumetric strain is approximately halved in the biaxial case [ � � 0] with stresses neighboring the critical stress. In this case, this volumetric strain is totally diverted to the free direction with 50% (½) of intensity (weight) instead of a 100% (1).
Moreover, since no expansion can occur in a triaxial stress state when compression in all directions is stronger than the critical value (i.e., < � ), the originally designated weights for such a case must now be changed to null. Table 1, which contains the weights for the interpolation function that governs ASR expansion in the numerical model (NM) by [15] must be updated, affecting nodes #3, #6, # 7, and #8 (column under ≥ 0), as well as the last two columns that must be updated to all zeros (for more details on this Table, see section 2.3 ahead and [15]).  [15].

Node
Weights (wn) # σl σm σk ≥ 0 σk =σu σk = fc ' In the previous Table 1, f c ' is the compressive strength, f t ' is the tension strength and σ u is the critical stress, as designated by the author [15]. From this point on σ u will be the preferred symbol for the critical stress.

General
Attending to the need for an update of the numerical model, improvements to the constitutive model of ASR expansion were proposed and underwent validation consisting in the numerical simulation of some pertinent experiments described in the field literature.
For the sake of simplicity, a mechanical-only approach was assumed, and full focus was allocated to the adjustment of directional expansion weights for the interpolation functions.

Numerical implementation
The simulations reported in this paper were conducted in COMSOL finite element (FE) software and MATLAB processing in conjunction. MATLAB routines were developed in order to evaluate the weights of ASR expansion along the principal directions. These weights are then automatically transferred to COMSOL for evaluation of the strains with FE analysis. Each weight is a directional multiplier of ASR total volumetric strain.
Because of the implicit relation between the ASR strains and the stress state, an iterative procedure is required to solve the problem. The ASR strain is input as the initial strain. The strains are updated at all Gauss points in every iteration performed during the FE analysis. In order to speed up processing, the initial solution (i.e., initial ASR strain state) is input in a way that is consistent with the applied stress state. In the case of the hydrostatic stress state, for example, the initial strain state is taken as in (Equation 3) next: Concrete was modeled as a linear elastic material, a reasonable option since the maximum applied stress on the specimens was approx. only 17% of ′ , in light of a virtually linear stress-strain behavior of concrete up to 30% of ′ . The 10-node quadratic tetrahedral element type was used for the concrete material: an automatic meshing algorithm built the elements mesh using the default normal elements size. Due to its quadratic nature, the displacement field polynomial was integrated exactly by means of a four-point Gauss quadrature.
The numerical model also entails a problem of circular (implicit) relation, which consists of the ASR strains being dependent on the stress state, which in turn, may be dependent on the ASR strains, according to the boundary conditions. This problem is solved using an iterative calculation of weak contributions of stresses and corresponding ASR strains that lead to convergence of the model employing the 'Double Dogleg' trust-region quasi-Newton method ( Figure 4). A relative tolerance of either solution or residual was set up as a termination criterion: 10 -3 for solution-based convergence criterion and 10 -6 for residual error-based convergence criterion.
Next in Figure 4, the analysis steps are outlined. Concrete parameters, such as compressive and tensile strengths, Young's modulus, volumetric strain and creep coefficient must be estimated and assigned to the model within the FE application prior to computation. For more information and all the theory underlying the ASR constitutive model, the reading of [1] is recommended, since the reproduction of it all in this paper was avoided, for length, scope, and focus reasons.

The modified triaxial weights in interpolation functions
The anisotropy of ASR expansion is modeled through a system of weights [15] that associates strains in the principal directions to a stress state. A method inspired by a quadrilateral bilinear finite element function is employed for the evaluation of the directional weight, according to a region of orthogonal stresses ( Figure 5). Each region has four nodes with their respective weights for the interpolation performed in search of , according to where the present principal stress is located. See [15] for more details on this procedure.
In this work, we propose a modification of the weights presented by [15]. The modified weights are presented in Table 2, with previous values of modified weights in parentheses. Notice that Table 2, compared to Table 1, dismissed the need for the last column that was present in Table 1.
At nodes #3, #6, # 7, and #8 for ≥ 0, the new ½ value represents the biaxial case where, in the directions other than the one being evaluated ( ), both stresses are at or above the critical threshold level ( ). In such cases, the volumetric strain is considered halved and entirely directed towards the less compressed direction. Likewise, at those same nodes, when the presently evaluated ( ) direction nears the critical value it is ascribed weight zero, which corresponds to a null ASR expansion.
The above weight allocation has been described in full detail in [15] and we briefly summarize it here in the following steps: • Assign as the stress in the present evaluated direction; • Locate the quadrant in Figure 5 that matches and , which represent stresses along the other orthogonal directions, indistinctly; • Pick the 1 , 2 , 3 , 4 weights assigned in Table 2 to the corner nodes of the quadrant, matching the present value; • Interpolate the node weights to determine the resulting present weight by means of a bilinear interpolation function (Equation 5); • Repeat the assignment of and ( ) as in the previous steps to find the corresponding of each i direction.
The above Equation 5 is a bilinear quadrilateral shape function. a and b are the orthogonal dimensions of the rectangular regions presented in Figure 5. and must be subtracted when their values are below this critical value, for consistency in the computation.
Numerical tests are presented in section 3, in order to ascertain if better accuracy is obtained with these changes.

VALIDATION TO EXPERIMENTAL DATA
The experiments that were selected from pertaining literature for the validation of the NM fitted the criterion of having a long-term sustained external or internal (e.g., due to the restraint of rebars) loads, which narrowed somewhat the collection of eligible experiments. The existence of enough data was also a key factor in the selection, for many published experiments lacked enough descriptive quality so volumetric strain, concrete strength, initial tangent modulus, property decay and other data could be inferred.
A coefficient was introduced to account for material properties degradation, including creep and the overtime reduction of the Young's modulus, following [1]. The relation between the damaged Young's modulus and the initial Young's modulus was given by the following Equation 6: where: ∅ is the material property degradation coefficient is the current (damaged) Young's modulus 0 is the initial Young's modulus Because of the orthogonality of applied stresses, the specimen axial directions were taken as the principal directions and shear strains and stresses were neglected.
In general, since the actual kinetics of ASR expansion is not in the scope of this study, the experiments were evaluated at their near final extinguished expansion plateau.
The validation of the updated numerical model (UNM) was first carried out with the source experiment [23]. Regarding the σ u critical stress, instead of the -9.7 MPa value taken in [23], -10 MPa was used for a round number and uniformity with literature [1], [15], [21], [26]. The total volumetric strain was 1.1%, as graphically displayed for 21 days for the free expansion reactive samples average (Figure 3). Creep and shrinkage cannot be straightforwardly measured on reactive concrete and must therefore be evaluated on non-reactive control samples or from analytical expressions. Nevertheless, creep and shrinkage are not the objects of any further consideration here since they have been subtracted from reported strains in the source experiment [23]. The Poisson's ratio was assumed to be 0.2 and constant.
The strains referred to here, as well as in the source experiment [23], are invariably mean axial strains, whose values at day 21 of ASR expansion are displayed in Table 3, according to the specimen identification.  For the [-1 -1 -1] stress state ( Figure 6) only the TM21 sample had all 03 axial strains reported at 21 days in the original [23] paper, so it was the only one simulated, in contrast with the other cases that have a pair of specimens with reported results. According to the data in Table 3, for the present case, the updated numerical model (UNM) and original numerical model (ONM) obtained similar results, despite the latter obtaining results slightly closer to the experimental data. The small difference is due to the UNM having a 'zero ceiling' of ASR volumetric hydrostatic strains at the critical stress level σ u , compared to the ONM in a hydrostatic stress state, which has no ASR volumetric 'zero ceiling' until failure at ultimate stress ( ′ ) (see Figure 7). The magnitudes of all errors, here defined as the difference between values measured in the experiment and values obtained in the numerical simulation, were somewhat low at four decimal digits or 10 −4 .  In the [-9 -9 -9] stress state (Figure 8) there was a considerable advantage of the UNM over the ONM in predicting strains. In fact, the ONM obtained strains are almost an order of magnitude (10x) larger than the experimental ones. Conversely, the UNM achieved more accurate results in this case, with errors as significant as 10 −4 to 10 −5 .  Figure 9 is a simulation plot of total volumetric strains at equal sustained compression stresses at 5 Mpa discretely spaced levels. It reveals a considerable error if the ONM un-updated weight values are employed to predict overall strain, considering that ASR strain is completely curbed at the critical hydrostatic stress levels [23] (see also Figures 3 and 7). For the [-9 -9 -1] ( Figure 10) stress state, strain errors of the UNM were lower than those of the ONM, and closer to the experimental values in 100% of reported strains. The simulation plot of biaxially equal sustained stresses at 5 MPa spaced discrete levels in Figure 11 does also elucidate that the consideration of an unvarying ASR volumetric strain level leads to a very diverse volumetric strain in comparison to the one modeled according to the UNM, which in such cases enforces a maximum ASR volumetric strain that is 50% of the one of stress-free expansion (see also Figures 3 and 12).

Report of simulations of other experiments
Other experiments that fitted the criteria of having a long-term sustained external or internal as well as sufficient descriptive information were used in the research validation process. Next is the list of experiments as referred to in the research and corresponding sources: • Multon and Toutlemonde [27] • Liaudat et al. [23] • Gautam and Panesar [28] • Fan and Hanson [29] • Wald et al. [30]

The Multon experiment simulation [27]
This ASR expansion experiment was performed with cylindrical specimens (13 cm diameter, 24 cm height) subjected to sustained compressive stress states. These were induced in the three orthogonal directions, either actively by means of a hydraulic jack along the axial direction, or passively by means of steel rings that offered a constraint in the radial direction. For this simulation the values of parameters were taken at the 450 days age. For this experiment, the parameters were E = 37.2 GPa, ′ = 36.5 MPa and, ′ = 3.7 MPa, ∅ = 4.0. The active axial compressive load was -10.0 MPa.

The Gautam experiment simulation [28]
In this experiment concrete cubes of 254 mm sides were subjected to multiple sustained stress states along the three longitudinal directions. The experiment went on until the exhaustion of ASR reactions. The results of this study confirmed that the multiaxial stress state orients the ASR swelling of concrete in an anisotropic way. The data of this experiment were collected from the graphed values at the reaction exhaustion age at 485 days. The parameters employed were E = 35 GPa, ′ = 43.4 MPa, ′ = 3.5 MPa, = 2.0 (the Young's modulus and tensile strength were not declared and were inferred according to [31] ). The property degradation coefficient was also unknown and taken as 2, following the usual value range found in [32], where the deformation module is estimated to be reduced to approximately a third in ASR affected concrete as compared to sound concrete. The volumetric strain was assumed to be 0.00434, according to the freely expanding specimens' graphic data.

The Fan experient simulation [29]
In this experiment, six 150 x 250 x 1500 mm beams were cast with reactive coarse aggregates and subjected to an accelerated expansion setup for 360 days. Some companion cylinders were also cast for evaluation of the mechanical properties. The beams were reinforced with either two 10 mm or two 16 mm bars centered 46 mm above their lower face. Parameters for this experiment were taken at the 360 days age. For this experiment, the parameters were E = 31.0 GPa, ′ = 35.6 MPa and, ′ = 3.5 MPa. The ASR volumetric strain was assumed to be 0.002505, according to the freely expanding specimens' graphic data. The two modeled specimens were not subjected to external loads. References to creep cannot be found in the article, hence the ' 'creep coefficient was estimated according to [32] to be 2.0.

The Wald experient simulation [30]
This experiment was performed with 33 (480x480x480mm) reinforced concrete cubic specimens. Different uniaxial, biaxial and triaxial arrangements of rebars were used. The nominal bar diameters also varied: 13, 19, and 22 mm were employed to impose varying degrees and configurations of forces and stresses within the specimens. In this experiment, stresses were applied passively by the rebar as the concrete material expanded. The data of this experiment were collected from the graphed values at the reaction exhaustion age at 450 days. The used parameters were E = 35.6 GPa, ′ = 35.6 MPa, ′ = 2.0 MPa, = 2.0 (the creep coefficient was not declared and was inferred according to [32], the other parameters were not declared in the cited source article but were found in Wald's PhD's thesis where the experiment is also described [21]). According to the freely expanding specimens' graphic data, the volumetric strain was assumed to be 0.01875. Table A.1 in Appendix A outlines the input data, as well as the results of the simulations carried out with the aforesaid experiments: Figure 13 shows the sum of absolute values of errors of both the ONM and the UNM according to the validation experiments. Since the order of magnitude of errors varied somewhat according to experiment source, the values have been normalized around their average for comparison in a same graph.

RESULTS
Regarding the uniaxially confined Multon experiment, as expected, the errors were identical. This is due to the fact that the proposed directional weights did not change for uniaxially confined/loaded models.
As for the Fan experiment, it had very minute strain and errors when compared to the other experiments. In these samples, loading was a 'chemical post-tension' reaction imposed by the reinforcement bars that resisted expansion at the bottom of the specimens.
The figure indicates that the other experiment models, which even had more samples (see Table A.1 in Appendix A), favored the UNM over the ONM regarding absolute errors.   The above figure suggests that, in this case, both the ONM and the UNM may tend to slightly overestimate strains in numerical simulations.
The following graph ( Figure 16) indicates that the UNM has a distribution of results that is more 'centered' around the experimental values than the ONM: Considering a t-student distribution, according to the data presented in Table A.1 in Appendix A, the hypothesis that the expected errors of the original numerical model (ONM) are not smaller than those of the updated numerical model (UNM), but are actually the same, must be rejected at a < 0.0000001 significance level, in a paired t-test with 77 degrees of freedom (for more details, check Table B.1 in Appendix B).
The correlation coefficient (R) of 0.94 was obtained for the UNM, compared to 0.87 for the ONM.

CONCLUSION
The updates to the Colorado numerical model undertaken in the present article incorporated new notions on volumetric strain behavior and directional ASR anisotropy, particularly when concrete is subject to sustained biaxial or triaxial compressive stresses. In accomplishing this, a more accurate numerical model was devised, according to the reported results.
The approach of this article was to apply the numerical model for mechanical purposes only, in order to introduce the mentioned updates in the simplest manner possible and as an example of a practical structural field application, with a small number of inputs that can be estimated. The conjunction of all the updates in a more comprehensive approach, including the kinetics of expansion, which is already part of the Colorado Model, is nevertheless feasible when necessary or convenient.
Finally, more numerical tests are needed and therefore recommended, through the production of new bespoke experiments to validate the numerical model updates here proposed, as well as its scope of application. The advancement of the understanding of the ASR phenomenon must be endeavored along with the consolidation of existing constitutive models that depict it, such as the one that was the object of this paper.    Table B.1 next presents the parameters of the paired t-test performed on the absolute errors of the UNM and ONM as presented in Table A.1 in Appendix A, to evaluate if the average of UNM errors is lesser than the ONM errors at a significant level. where: 0 value zero, used to test the hypothesis that the averages are no different ̅ is the average of the differences between pairs is the standard deviation of the differences between pairs is the number of data pairs is the number of degrees of freedom is the t statistic = ̅ − 0 √ is the p-value, here denoting the probability of wrongly rejecting the hypothesis that the average of the two samples is the same. The p-value was obtained in [33].