A multi-physics modelling based on coupled diffusion equations to simulate the carbonation process

ISISE, Guimarães, Portugal Abstract: Carbonation is widely recognized as a cause of significant pathologies in reinforced concrete structures and different modelling strategies are presented in literature the simulate the phenomenon evolution. In opposition to the deleterious effect in reinforced concrete, for historical mortar made with aerial lime, the carbonation is essential for the hardening process. For both materials, carbonation process presents similarities. This work presents the background/implementation of an algorithm for a multi-physics simulation of the main fields associated with the carbonation process. This modelling was previously validated in literature. A 1D algorithm is implemented, using the Finite Difference Method. Its feasibility is demonstrated through the simulation of results presented in the literature. A parametric study is also shown considering the main parameters involved, important observation regarding the influence of the parameters on the carbonation finite difference method. Resumo: A carbonatação é amplamente reconhecida como causa de patologias significativas em estruturas de concreto armado e diferentes estratégias de modelagem são apresentadas na literatura para simular a evolução do fenômeno. Ao contrário do efeito deletério do concreto armado, para argamassas históricas feitas com cal aérea, a carbonatação é essencial para o processo de endurecimento. Para ambos os materiais, o processo de carbonatação apresenta semelhanças. Este trabalho apresenta a implementação de um algoritmo para uma simulação multi-física dos principais campos associados ao processo de carbonatação. Essa modelagem foi validada anteriormente na literatura. Um algoritmo 1D é implementado, usando o Método das Diferenças Finitas. A viabilidade é demonstrada através da simulação dos resultados apresentados na literatura. Também é mostrado um estudo paramétrico, considerando os principais parâmetros envolvidos. São detalhadas observações importantes sobre a influência dos parâmetros na profundidade da carbonatação. Palavras-chave: modelagem multi-física, carbonatação, concreto, cal aérea, método


INTRODUCTION
In reinforced concrete structures, over their service life, external deleterious substances may penetrate (e.g. carbon dioxide, chloride, sulfate) into concrete and find the rebars, altering the pore solution composition into aggressive conditions [1], [2]. Specifically considering the penetration of CO2, the main consequence is the occurrence of carbonation process. Carbonation is a major cause of concrete structures deterioration, therefore its evaluation should be carefully considered in the durability design of reinforced concrete structures [2]- [8].

Multi-physics models for simulation of carbonation process
In this section, some general information about different multi-physics models are presented. Different strategies were chosen with a large number of parameters involved. A model based on thermo-hygro physics is presented by Ishida and Li [29], Ishida and Maekawa [30] and Ishida et al. [31]. The model coupled with moisture equilibrium/transport, temperature dependent parameters were also assumed [29].
The modelling recommended by Bary and Sellier [2] is based on macroscopic mass balance equations for water, the carbon dioxide contained in gaseous phase and the calcium present in the pore solution [2]. These equations rule the diffusion and permeation processes of the three variables: saturation degree, carbon dioxide partial pressure and calcium concentration in pore solution [2].
Coupling between carbonation and chloride diffusion is explored in the context of both homogeneous and heterogeneous concrete models [32], this model considers the coupling of heat, relative pore humidity, chloride, and carbonation fields [32].

Models based on simple diffusion equations
The multi-physics models from literature summarized in this section are based on the CO2 modeling using a simple diffusion equation.
The study done by Burkan Isgor and Razaqpur [4] presents a general diffusion model implemented in Finite Element Method) for thermo-carbo-hygro simulations, the modelling strategy is decoupled with mechanical aspects. In another research, Steffens et al. [33] adopted a diffusion model to study the carbonation. The model proposed by the authors [33] combines results of extensive studies by Bunte and Rostasy [34] on diffusion of CO2 in different types of concrete and the modeling of the reaction kinetics of carbonation by Saetta et al. [35] with a coupled temperature and moisture modelling strategy, it considers that the CO2 penetrates into the concrete mainly gaseous by diffusion through air-filled pores [10], [33]. Additionally, Meier et al. [36] presents a similar mathematical, in continuity, Peter et al. [37] include the hydration reactions in strategy presented by Meier et al. [36].
The mathematical modelling of carbonation process for concrete and mortars is complex [10]. The first model adopted to simulate the aerial lime mortar carbonation was shown by Ferretti and Bažant [21], for that reason, and considering the complexity involved and the requirements of information about the material, the strategy presented by Ferretti and Bažant [21] has been chosen for this work (description may be see next section). This modelling strategy has been chosen because it is the only modelling presented in literature to simulate the carbonation process in aerial lime mortar.

• Model of Ferretti and Bažant [21]
This section presents the mail information about the modelling strategy presented by Ferretti and Bažant [21] This model is a multi-physics coupled model involving four fields, humidity, heat, pollutant flow (CO2) and reaction [10]. The numerical model for deterioration was developed, considering the characterization of the concrete and the environmental conditions [10]. The present mathematical model is based on studies conducted by Saetta [38], Saetta et al. [35] and Saetta et al. [39], Saetta and Vitaliani [9]. The modelling presented herein was also experimentally validate in previous studies conducted by Saetta and Vitaliani [9]. The equations used to model the phenomena are presented below.
The humidity field [21], [28], [39], [40] is governed by: where: h reads as the humidity (%), α2 means the parameter related to the water generation during the carbonation process (more details are presented in sequence), Cw is the diffusion of water (mm 2 /day) and R is the degree of chemical reaction (%) [9]. The coefficient α2 is related to the maximum content of calcium carbonate [CaCO3]max = Pmax (the term P represents the formation of CaCO3), which depends mostly on the material composition. The angular coefficient of the sorptiondesorption isotherm k, dependent fundamentally on temperature [9]: where: PM represents the molecular weight of the molecule.
This equation derives from the kinetic of the carbonation reaction [9], [10]. For each CaCO3 molecule produced by the reaction, a molecule of water is also created, consequently, considering the molecular mass [9], [10]: where: M denotes the molecular mass of the molecule specified.
Accordingly, in terms of mass per unit volume, Equation 3 may be redefined as: where: Vcls is the considered volume of element, consequently [9]: where: dw represents the water content variation per unit volume and unit time, while dP* is the variation of the calcium carbonate concentration (such variables are both expressed in kg/m 3 ) [9].
Assuming, for illustration, that Pmax = 0.0096 kg/m 3 and kh =1 m 3 /kg, then α2 = 0.0017. The analytical determination of the coefficient α2 proves somewhat uncertain because it is difficult to unequivocally assign the coefficients kh and Pmax. The carbon dioxide diffusion field is governed by: where: α3 is the parameter related to gas consumption during the carbonation process, Dc represents the diffusion of the gas (CO2 for instance) and c is the gas concentration (%). Equation 7 is based on the second Fick's law. Parameter α3 is influenced by the carbonation, and based on a considerations on the chemistry of this reaction. The reaction field is governed by: where: T means the temperature (K).
where: Cw,rif and Dcrif are the diffusivities in standard conditions for water and CO2 respectively (mm 2 /day). The function f1*(h) is defined as: where: α, hc and n were already defined in Equation (α = 0.05, hc = 0.75 and n = 6).
For gas diffusion phenomena (for instance CO2), the following expressions are given: where: T0 means the reference temperature (296 K), Ea is defined as the activation energy (kJ/mol), R reads the universal gas constant (J/mol×K), and T is temperature. Function f3 (te) is related to the concept of equivalent age, and it is defined as: where: te is the equivalent age (days). The diffusion process is reduced with the decreasing of the porosity, to simulate this phenomenon a function f4(R) may be defined as [43]: where: parameter ζ varies between 0 and 1, and measures the slowing of diffusion phenomenon due to reduction of the porosity [38]- [40]. Ferretti and Bažant [21] adopted ζ = 0.3, meaning that a reduction of 30% for the diffusivity value occurs with the total reaction, further information and other values for this value can be seen in literature [10], [28], [40]. Function R 4 F describes the influence of temperature on the progress of the chemical reaction, it may be expressed as [21], [33]: where: A is the impact number [33] and Ea means the activation energy (kJ/mol).
After the description of the modelling approach, some considerations should be addressed. Considering these facts, in literature the number of works that use the multi-physics modelling strategies are raising, however more development are necessary [10]. In general terms, the multi-physics modelling strategy may be still considered a complex simulation, since it required advanced computational tools and a deep understand of the material behavior and phenomena coupling and the process of obtaining parameters is still a challenge.
Specifically considering the modelling strategy implemented in the present work [21], it may be considered a powerful tool, however some limitations are implicates: e.g. it does not consider others couplings with aggressive agents (e.g. chlorides, sulfates) and the movement of fluids are based are based on simple processes. A general review about the modelling strategies and the complexities involved may be seen in literature [10].

Implementation of multi-physics modelling: A brief discussion
This section describes the strategy of implementation used to simulate the carbonation process. The coupled model is based on the work of Ferretti and Bažant [21], and it is implemented using the Finite Difference Method (FDM) [44].
The decoupled humidity field for 1D using FDM is previously presented in Oliveira et al. [45]. For carbon dioxide field, the mathematical implementation is similar, thus for the sack of brevity, the implementation of both fields are herein omitted, further details may be seen in literature [45].
All simulations of present work assume constant temperature, for the sake of brevity, no further details will be herein. The model has been implemented in 1D, axisymmetric and 2D conditions, further information may be found in and Oliveira [10]. For brevity and due to the similarity with the 1D implementation, the 2D and axisymmetric implementations are not presented herein, more details can be seen in Oliveira [10]. The chosen notation for reaction field R is the same adopted in Oliveira et al. [45], where "i" represents the node, and "n" reads as the time step [10].
The evolution of reaction field over time in Finite Difference Method may be expressed as:

Substituting Equation 22 in Equation 21
, it may be obtained: This equation may be reorganized as: Rearranging the equation, the expression for inner nodes, it reads: The definition of ( ) R 4 F T will not be the used herein, because the temperature is assumed constant.
For coupled humidity field, considering the development done in Oliveira [10] (decoupled humidity field) and the humidity coupled equation with reaction field (α2(∂R/∂t)), the final equation for humidity field for inner nodes (except in the boundary and symmetry) is given by: As already cited, the effective humidity and CO2 diffusivity depends on the other fields (R, h and c) and Cw, from Equation 26 is needed to solve numerically the problem. Herein, details about this development is omitted, for the sake of brevity. More details may be seen in literature [10].
Considering the equation for carbon dioxide, with the coupling term related to reaction field (term α3(∂R/∂t)), the following expression is obtained in case of the inner nodes (except in the boundary and symmetry): For the boundary and symmetric nodes, Equations 26 and 27 may be adapted as also done for the decoupled humidity field in appling Neumman boundary condition, or an imposed/fixed value adopting the Dirichlet formulation [44].
The terms of equation expressed in time "n+1" are in the left side, while terms in time "n" are shown in right side. In sequence, general considerations regarding the numerical simulation of humidity field for 2D condition in FDM are following discussed. The current implementation considers the field equation, which may be revised adopting the application of the chain rule [46]: After this discussion, Figure 2 schematic represents the solution of the coupled system of equations. The symbol (^) denotes a matrix or a vector, while the superscript/subscript denotes the iteration/step. It should be highlighted that, for tolerance, three different values (tolerhumidity, tolerreaction and tolercarbon_dioxide) were chosen, since the order of magnitude of each filed was dissimilar, further information may be seen in literature [10]. With the purpose of validate the implemented modelling, the 1D results presented in the work of Ferretti and Bažant [21] were numerically simulated. The authors [21] studied the failure of Pavia Tower in Italy [10], [47]- [49]. In order to validate the independence of the analysis from time and spacing discretization, diverse values of time steps and nodal distance were tested. The results presented small discrepancies further details may be seen in Oliveira [10].

RESULTS: SENSITIVITY ANALYSES
Several sensitivity analyses were done during the study, changing the input model parameters, more information can be seen in Oliveira [10]. For these, a 1D wall with two meters (2 m) length, in contact with the environment in the two boundaries was simulated. The obtained results are presents for one meter, because of the symmetry condition (at 1 meter). For the numerical simulations, the distance "zero" (the first node) represented the node in contact with the environment, and the node located at distance one meter (1 m) was the symmetric one, more information can be seen in Oliveira [10]. In order to illustrate the obtained results, the analyses for 50 years are shown. The parametric analyses regarding the CO2 initial diffusivity (Dc,rif) was studied and wholly parameters were based on the Ferretti and Bažant [21].
The parametric results for CO2 initial diffusivity are shown, since the large influence of this parameter on the final carbonation results. A range from one hundred times higher and smaller the value adopted by Ferretti and Bažant [21] was considered (Dc,rif = 2400 mm2/day, this is cited as the next figures as "typical value"). Time and nodal discretization were adopted with constant values (Δt = 1 day and Δx = 1.25 cm). For the initial conditions, the considerations adopted above were maintained. Figure 3 shows the results for CO2 concentration and Figure 4 shows the results for the reaction field.
The carbonation profile was expressively affected by the initial diffusivity of CO2. This effect was expected. After 50 years, the wall presented distinct reaction profiles, according to the used parameter. For instance, the carbonation front was located on the first 5 cm of the wall, if the diffusivity Dc,rif = 24 mm 2 /day (the smallest value); while the wall presented around 80 cm with R = 1, for Dc,rif = 240000 mm 2 /day. These cases with extreme values illustrated, the importance of the parameters' selection for the numerical model, further information can be seen in Oliveira [10].  Results for the humidity field are shown in Figure 5. The influence of Dc,rif was on the carbon dioxide and reaction fields. On the contrary, the parameter did not affect significantly the results in terms of humidity field. The tendency for 50 years was also observed for different ages, further information can be seen in Oliveira [10]. The same range of values for the initial carbon dioxide diffusivity (Dc,rif) was adopted, and in order to illustrate the behavior for longer ages, further information can be seen in Oliveira [10]. A node located at 40 cm from the boundary in contact with the environment was selected. The results for reaction field, for this node, in analyses over 500 years (~182500 days) are shown in Figure 6. For horizontal axis, the upper scale in Figure 6 indicated in "years" and the bottom is indicated in "days" to facilitate the understanding, further information can be seen in Oliveira [10].
The results presented in Figure 6, indicated the significant dependence of the reaction profile for this depth, according with the value for the initial the CO2 diffusivity. The same dependence happened for the others depths, further information can be seen in Oliveira [10].
Unambiguously for the results presented in Figure 6, when the highest value for the carbon dioxide diffusivity was adopted, Dc,rif = 240000 mm 2 /day, or equivalently one hundred time higher than the value cited by Ferretti and Bažant [21], after around 30 years (~1095 days), the reaction was completely (R = 1). On the other hand, for the Dc,rif = 2400 mm 2 /day or equivalently the value cited by Ferretti and Bažant [21], the reaction reached values close to one (R ≈ 1) after around 400 years (~146000 days), further discussion about this issue can be seen in Oliveira [10].