Probabilistic chloride diffusion modelling in cracked concrete structures by transient BEM formulation Modelagem

: Reinforcement corrosion is a concern in the structural engineering domain, since it triggers several pathological manifestations, reducing the structural service life. Chloride diffusion has been considered one of main causes of reinforcements' corrosion in reinforced concrete. Corrosion starts when the chloride concentration at the reinforcements interface reaches the threshold content, leading to depassivation, whose assessment of its time of starts is a major challenge. This study applied the transient Boundary Element Method (BEM) approach for modelling chloride diffusion in concrete pores. The subregion BEM technique effectively represented the cracks inherent to the material domain

computation cost is a concern. Nevertheless, the computational effectiveness of BEM enables a probabilistic description via MCS. Very few studies have adopted this coupling scheme [29], [30], [34], thus motivating this research, which deals with the probabilistic modelling of chloride ingress in concrete pores.
The transient BEM formulation for diffusion handles the phenomenon description, whereas its coupling with MCS accounts for the uncertainties modelling. High-order boundary elements perform spatial approximations and the constant interpolation scheme approximates the temporal variables. BEM model accounts the presence of cracks and their influence upon the chloride flux velocity, which is the novelty in this study. Because concrete is cracked by loading and cure actions, the cracks representation enables an accurate and realistic problem modelling. Environmental effects such as chloride binding, temperature, and loadings were included in the modelling, which is a contribution of the study. Three applications demonstrated the accuracy and robustness of the proposed framework. The BEM results were compared against numerical, analytical, and experimental responses from the literature in applications one and two. The third application described the probabilistic modelling for tidal zones and highly aggressive environments.

Integral equations and BEM formulation
The following differential equation governs the transient diffusion problem (Equation 1): (1) where u represents the potential, t is time and k denotes the diffusion coefficient.
The solution of the differential equation requires boundary conditions regarding potential (u) and flux (q), as in Equations 2-3: / q u η = ∂ ∂ . In the present study, the potential represents the chloride concentration at a given point x for time t and q chloride flux.
Wrobel [35] demonstrated Equation 1 provides a boundary integral representation applying one of the following approaches: Laplace Transformation, Dual Reciprocity Method, or Time-Dependent Fundamental Solution. This study has adopted the latter; therefore, Green's second identity is applied to Equation 1 towards the boundary integral representation. The application of weighted residual technique, the integration of the resulting kernels by parts, and the classical limit process lead to the following integral equation: where x indicates the field points and ξ denotes the source points, 0 t is the initial time and F t is the observation time, * u and * q are the time-dependent fundamental solutions to potential and flux, respectively, and c is the BEM free term. The free term value depends upon the source point position; it equals 0.5 for source point on smooth boundaries, such as those used here. Equation 5 contains a domain integral, which can be transformed into a boundary integral by the Dual Reciprocity Method. Nevertheless, for the sake of simplicity, this term has been assumed nil, and Equation 5 can be rewritten as Equation 6.
For two-dimensional problems, the time-dependent fundamental solutions in Equation 6 are as follows (Equations 7-8) [35]: ( ) where F t t τ = − ( F t is analysis time), r indicates the distance between the source ξ and field x points, and , k k r r η η ∂ = ∂ .

Algebraic representation
The numerical solution of Equation 6 requires spatial and temporal discretisation. Dividing boundary Γ into N e boundary elements and time span 0 F t t − into NT time steps, and inverting the order of integration in Equation 6, the following discretised equation can be obtained: where NT i u is the potential value in source point i at time step NT, and 0 n t and n F t are initial and final times at time step n, respectively. Equation 9 can be rewritten assuming a constant approximation along time and the positioning of collocation points at smooth boundaries. Therefore, where: The singular kernels in Equations 11 and 12 were regularized by the semi-analytical method known as Singularity Subtraction Technique [36], whereas the non-singular ones were integrated by the Gauss-Legendre quadrature. Equation 9 also assesses the chloride concentration at the domain. In such case, the free term equals unity since the collocation point is positioned at the domain. Thus, the chloride concentration can be evaluated by the following algebraic representation: Equation 15 can be solved by the time-marching process previously addressed, providing the chloride concentration at internal points ( int u ):

SUBREGION BEM TECHNIQUE
The transient BEM formulation enables the diffusion modelling of homogeneous domains, i.e., with a single diffusion coefficient. Nevertheless, the formulation represents the diffusion in non-homogeneous materials when coupled to the classical subregion BEM technique, which discretises the non-homogeneous domain into homogeneous piecewise regions where the compatibility of potentials and equilibrium of fluxes are enforced along the materials interfaces, as depicted in Figure 1. The algebraic operations required in the subregion BEM technique can be performed straightforwardly in a twodomain problem ( Figure 1). In this case, Equations 11 and 12 lead to the following algebraic representation: where the superscript index indicates the domain and i denotes the material interface.
Since potentials and fluxes are indeterminate along the material interfaces, Equation 17 requires the following compatibility and equilibrium conditions for a proper solution (Equations 18-19): The coupling of Equations 17-19 leads to the following system: Equation 20 can be solved by simply enforcing the boundary conditions, which triggers the classical columns change procedure of BEM. The subregion BEM technique provides generality to the numerical technique because nonhomogeneous materials can be properly handled, and can also be applied for representing cracks. For such a purpose, the material interfaces contain the crack surfaces. Therefore, the entire domain requires a division into subregions, where the compatibility of potentials and equilibrium of fluxes must not be enforced along the crack surfaces, which become ordinary boundaries. The subregion approach does not require hypersingular integral equations, thus simplifying the integration process, and its application is straightforward when the kernels formulations are known.

Transient diffusion modelling in a square domain
A concrete column of square cross-section is exposed to an environment whose chloride ions ( Cl − ) concentration is 1% of the concrete weight, as illustrated in Figure 2a. For the sake of simplicity, the problem symmetry properties were used. Therefore, a quarter of the domain was modelled and the symmetry boundary conditions assured the domain continuity ( Figure 2b). The chloride diffusion coefficient is 12 10 − m 2 /s ( 0.3154 cm 2 /year).
The application assesses the corrosion time initiation at the point depicted in Figure 2b, whose cover thickness is 5 cm for 20 years' structural lifetime. The analytical solution presented by Bitaraf and Mohammadi [13] and the numerical responses based on FEM, Finite Difference Method (FDM), Element-Free Galerkin (EFG), and Finite Point Method (FPM) are references herein and enable the verification of the BEM formulation accuracy.
The BEM mesh is composed of 40 isoparametric lagrangian elements with linear approximation, which leads to 44 collocation points (Figure 2c). The FPM modelling required a grid with 185 nodes, and the FEM employed 328 triangular elements with 185 nodes -the same number of nodes was used by the EFG method. The FDM results were provided by a 11x11 grid with 121 nodes [13], [37]. The threshold chloride content for reinforcement depassivation is 0.10% of the concrete weight [13]. Figure 3 shows the results achieved by BEM, whose time was discretised within 20 time steps and the space kernels were integrated by 10 integration points, and the reference responses. The results in Figure 3 show the superior performance of BEM over other numerical approaches; its responses are in excellent agreement with the analytical predictions. BEM achieved such results with the coarsest mesh among the considered numerical methods, which demonstrates its accuracy.
BEM also accurately predicted the corrosion time initiation, i.e., the time span in which the chloride concentration at 5 cm cover is equal to the threshold content. The structural lifetime assessed by the analytical approach was 10.50 years. BEM predicted 10.65 years, showing a 1.4% difference in relation to the reference. The predictions of EFG, FPM, FEM, and FDM were, respectively, 7.58 years, 6.70 years, 4.67 years, and 4.33 years [13], highlighting the accuracy of BEM in diffusion modelling.

Diffusion modelling in a cracked material
The application deals with the chloride diffusion modelling in a cracked material. The responses provided by BEM were compared against the experimental results of Şahmaran [38], who conducted tests in prismatic specimens (355.6 x 50.8 x 76.2 mm) molded with mortar ( / w c = 0.485) and reinforced with three levels of a 1 mm diameter steel reinforcement mesh in a 6 mm grid. After 43 days' curing in water, the specimens were pre-cracked using a 4-point bending test, which led to a single crack in the mortar samples. The unloaded specimens were immersed in a 3% NaCl solution for 30 days, after which their chloride concentration along the cracked area was quantified. According to Şahmaran [38], the widths (w) of the cracks were narrower at a deeper level of the specimens than at the surface (V-shaped). An optical microscope measured such widths at five different points, providing a 29.4 μm average width and an 18.7 mm crack depth.
BEM used a two dimensional specimen of 356.0 mm x 76.0 mm cross-section for numerically reproducing the experiment. The crack was positioned at the cross-section centre and modelled by the subregion technique as a straight one parallel to the specimen's height (Figure 4a). The chloride content along the specimen's lower surface was 0.51% (per weight of cementitious materials), as suggested by Du et al. [39], who analysed the experimental data presented in [38]. The chloride concentration at the upper surface was nil, whereas the vertical surfaces were sealed with resin, leading to nil flux values. Therefore, such boundary conditions were enforced in the numerical modelling. The chloride diffusion coefficient in the mortar was 12 4.746 10 x − m 2 /s ( 0.41 mm 2 /day), considering ratio / w c = 0.485 [39].
Two subregions discretised the cross-section, requiring 150 isoparametric lagrangian elements of quadratic approximation and 320 nodes. The crack discretisation was composed of 7 nodes and positioned at the subregions interface. Figure 4b illustrates the boundary mesh.
Because of the capillarity effect along the crack surfaces, which appear due to their V-shape, the chloride concentration prescribed along the crack surfaces started with 0.51% at the crack mouth and reduced progressively at a 0.07% rate for each node, finishing with 0.09% at the crack tip. For the sake of simplicity, the reinforcements were not represented in the numerical model. The red dots in Figure 4b represent the measuring points, which are equally spaced by 5 mm along the cross-section. The experimental analysis accounted for an uncracked specimen (w=0), which served as testimonial; the case was considered in the numerical modelling. The same BEM mesh was used in the testimonial modelling, in which both compatibility of chlorides concentration and equilibrium of fluxes were enforced along the crack nodes. Figure 5 displays the experimental results of Şahmaran [35] and the numerical responses achieved by BEM. The numerical modelling required 10 integration points for spatial integrations and 30 time steps in a 30-day time exposure to the NaCl solution. The numerical results provided by BEM are in good agreement with the experimental responses of Şahmaran [35] for the cracked specimen, although they underestimated the chloride concentration values for the 5 < depth 15 < mm cover range in the non-cracked specimen. Nevertheless, the numerical modelling represented the crack influence on the chloride concentration distribution along the material domain, and the non-requirement of domain mesh by BEM provided high accuracy for the internal fields assessment. [40] that uses random numbers and limit state function simulations for uncertainty quantification purposes. In structural engineering applications, it usually assesses the probability of failure ( f P ) of complex structures and structural systems subject to randomness. The probability of failure of a system can be determined by the integral of the joint probability density function ( X f ) evaluated over the failure domain ( f Ω ). However, the X f function is implicit in complex engineering problems, hampering an explicit evaluation of the probability of failure, which can be assessed by simulation techniques that evaluate limit states punctually. In MCS, the simplest simulation technique, an indicator function [ ] I x indicates if a simulated condition belongs to failure or safe conditions. The function is unitary when the simulated point belongs to the failure domain, and nil otherwise. Therefore, the probability of failure can be assessed by the following integral:

PROBABILISTIC MODELLING: MONTE CARLO SIMULATION Monte Carlo Simulation (MCS) is a numerical simulation technique developed by Metropolis and Ulam
Equation 21 can be estimated by a finite number of samples, s n , according to: in which f n is the number of samples in the failure domain and ˆf P is the estimated probability of failure.
Each evaluation of the indicator function involves a simulation of the limit state function, which defines the frontier between failure and safe domains. In the present study, the mechanic-probabilistic modelling is based on the limit state of corrosion initiation ( G ), defined by the difference between the threshold chloride content ( th C ) and the chloride concentration ( C ) on the reinforcement interface, at point x and time t , as follows: (23) The chloride threshold has been often represented by deterministic values accounting for the class of environmental aggressiveness and concrete composition. Since the quantity displays an important randomness behaviour, th C was modelled as a random variable in this study.
The probabilistic modelling performed here accounts for the environmental influence (e.g., temperature and chloride binding) on the concrete, as well as that of loads on the structure. Such effects were considered in the evaluation of the diffusion coefficient of concrete ( 0 k ) by introducing penalization factors corresponding to chloride binding ( 1 F ), temperature ( 2 F ), and acting loads ( F ), as follows [41]: • Temperature factor ( ) (27) where: U = activation energy of the diffusion process (J. • Acting loads factor ( ) ε ε (28) where:

Probabilistic scenarios of chloride diffusion in a typical bridge beam
The application assesses the probability of reinforcement depassivation accounting for the influence of climatic conditions (temperature), chloride binding capacity, and effects of acting loads on an I shape cross-section commonly employed in bridge beams. It considers the beam a part of reinforced concrete structures exposed to tidal splashes.
The cross-sections largest dimensions are 90 cm x 40 cm, reinforced by eight 20 mm diameter rebars. Because of the high environmental aggressiveness, the concrete cover is 50 mm [43]. The beam connects the structural system by its upper boundary, which composes the structure deck; therefore, the chloride flux at this boundary is nil. The remaining cross-section boundaries, i.e., the lower and vertical are exposed to the chloride attack ( Cl − ).
The subregion BEM technique represents the macrocracks in the numerical modelling. Their lengths are h and h 2 , idealized as perfectly straight. For the sake of simplicity, the chloride concentration at the macrocrack surfaces equals that at the beam surface ( 0 C ). Figure 6 illustrates the cross-section dimensions and macrocracks positions. The cracks were positioned at the lower cross-section boundary towards representing the mechanical degradation in the concrete caused by tensile stress. The diffusion coefficient in region 1 accounts for the load factor (Equation 29); therefore, the analysed beam can be considered simply supported. The other effects (temperature and chloride binding) are equally considered in the diffusivity of regions 1, 2 and 3.
The probabilistic modelling accounts for three random variables. The w/c ratio data were adopted from Stewart et al. [44], who suggested 0.40 as the mean value for exposure to tidal or splash zones based on maximum w/c ratio specified in AS 5100.5. The coefficient of variation (C.O.V) assumed was 10% for concrete manufactured under controlled conditions, and the mean and C.O.V values of the surface chloride concentration ( 0 C ) for splash zones were established by Val and Stewart [45], from data on offshore and onshore RC structures along the Australia coast reported by Collins and Grace [46] and McGee [47]. The chloride threshold content ( th C ) was based on data provided by Vassie [48], who studied the percentage of corrosion cases in UK bridges in relation to the total chloride concentration at the reinforcement level [45]. Table 1 shows the input parameters of the model. The numerical modelling adopts the transient BEM approach, whose mesh at the external boundaries is composed of 172 isoparametric lagrangian elements of quadratic approximation, leading to 384 collocation points. Since the reinforcements are assumed non-porous media, the nil flux condition is enforced along their boundaries. Eight isoparametric boundary elements with quadratic approximation discretise each reinforcement boundary and 10 Gauss-Legendre points perform the space integrations, whereas 15 time steps discretise the analysis time.
Eight probabilistic scenarios enabled the assessment of isolated and joint influences of the environment and loadings on the chloride diffusion in concrete pores. Table 2 shows the scenarios considered, where dots indicate the presence of effect in the modelling. The even scenarios account for crack length with h=2 cm.  Figure 7 shows the evolution of the probability of failure along time for a 75-year structural service-life. Scenario 1, which disregards the influence of the effects, overestimates the probability of reinforcement depassivation when compared to scenarios that include the chloride binding capacity isolated and jointly with the influence of temperature or macrocracks on the concrete. At 50 years, the overestimation is 33%, 15%, and 6% for scenarios 3, 4, and 5, respectively, which is expected, since the effects of cracks and temperature accelerate the chloride diffusion, whereas chloride binding reduces the process, although in different proportions. The results also show the importance of representing environmental and mechanical influences on the diffusion process of chlorides in concrete for avoiding overestimation or underestimation of the design service life. Realistic frameworks such as those proposed here are fundamental for an accurate structural service life assessment.
The probability of failure significantly increases when the modelling accounts for the chloride binding capacity (see the curves for scenarios 1 and 3). This effect represents the concrete capacity for binding the free chlorides into its pores, and can be explained by the approximately 30% reduction in the 1 F factor (Equation 26). On the other hand, the temperature growth increases 28% the concrete diffusivity, whereas acting loads penalize this variable by 396%. Such environmental actions increase the probability of failure, as observed in scenario 5 in comparison to scenario 3, and scenario 7 in comparison to scenarios 3 and 5. The loading effects influenced the diffusion phenomenon more strongly than the temperature in the simulations, since the damage was complete. Nevertheless, different behaviours can be displayed if usual structural service loads are considered.
The loading effects assume a major importance in the diffusion modelling (see Figure 7), since they uniformly penalize the concrete diffusivity along the cross-section domain. However, cracks also increase the chloride diffusion velocity. Because their discontinuities introduce preferential ingress paths, the probability of failure increases when cracks are considered. The presence of 2 cm long cracks reduces 20% of the probabilistic structural service, which demonstrates their importance in the modelling.
The cracks influence can be further investigated. For this purpose, scenario 5 (h=0) are compared against scenario 6, in which the probabilistic modelling accounts for crack lengths increasing progressively with 1 ≤ h ≤ 4 cm. Figure 8 displays the numerical responses provided by the transient BEM approach. The numerical results in Figure 8 show cracks with h equals 1 cm and 2 cm exert a small influence on the evolution of the probability of failure. However, the probability of failure significantly increases when h assumes 3 cm and 4 cm, which are higher than half of the cover thickness. Obviously, this behaviour accounts for severe environmental aggressiveness and 5 cm cover thickness. For 50-year time span, the probability of failure for 4 cm long crack is 35% higher than the case in which cracks are disregarded. Fib Bulletin 34 [53] recommends an 10% maximum depassivation probability for exposition to seawater. In such a context, the avoidance of cracks would extend the structural service life by 10 years, on average.