Corrosion Prediction with Parallel Finite Element Modeling for Coupled Hygro-Chemo Transport into Concrete under Chloride-Rich Environment.

The prediction of the chloride-induced corrosion is very important because of the durable life of concrete structure. To simulate more realistic durability performance of concrete structures, complex scientific methods and more accurate material models are needed. In order to predict the robust results of corrosion initiation time and to describe the thin layer from concrete surface to reinforcement, a large number of fine meshes are also used. The purpose of this study is to suggest more realistic physical model regarding coupled hygro-chemo transport and to implement the model with parallel finite element algorithm. Furthermore, microclimate model with environmental humidity and seasonal temperature is adopted. As a result, the prediction model of chloride diffusion under unsaturated condition was developed with parallel algorithms and was applied to the existing bridge to validate the model with multi-boundary condition. As the number of processors increased, the computational time decreased until the number of processors became optimized. Then, the computational time increased because the communication time between the processors increased. The framework of present model can be extended to simulate the multi-species de-icing salts ingress into non-saturated concrete structures in future work.

1. Introduction

One of the main long-term durability problems of reinforced concrete structure is the corrosion of reinforcing bars (rebar) in concrete. The corrosion can be resulted from several necessary conditions, and one of them is a high concentration of chloride ions near rebars, where the chloride ions come from deicing chemicals used in the winter on roadways and parking structures or seawater for the off-shore structures. Once the chloride content on the surface of steel reinforcement reaches a threshold value and the moisture and oxygen are sufficiently provided, the corrosion of steel bar is initiated [1,2]. The corrosion initiation period is the time during which substances such as water, chloride ions, oxygen and carbon dioxide penetrate through the concrete cover. The length of period depends on the resistance of concrete to the transport processes and the severity of the environmental conditions they are exposed to [3,4,5,6]. In practical engineering, the chloride contaminated RC structures have many internal cracks due to rebar corrosion and the prediction of crack width and propagation is very important for safety and serviceability [7,8]. Recently, the severity of corrosion damage was demonstrated through the corrosion-induced modeling with a statistical approach on bridge structures [9,10,11]. In concrete bridges, there are various concentrations of moisture and chemicals on the top surface of bridge decks due to the scatter of de-icing chemicals when they are applied. Especially, moisture and chloride concentrations have spatial distributions and seasonal variations. The service life of concrete is affected by a combination of material properties and microclimate. The microclimate is a term for the climatic conditions at the concrete surface or very close to the surface. This condition on the surface has a more decisive effect on the conditions inside the concrete than most other parameters [12]. The model of microclimate such as environmental humidity and temperature was proposed by Bazant et al. [13].

To simulate more realistically the durability performance of reinforced concrete structures, sophisticated numerical methods for fully coupled moisture and chloride transport processes and reliable material models for the transport parameters are needed. In addition, in order to predict accurately the moisture and chloride distributions in a large reinforced concrete structure, a large number of fine finite element meshes are needed. For the use of the large meshes in numerical analysis, in this study, a parallel algorithm will be employed to a parallel processing system with up to 2000 processors. To consider more realistic boundary conditions on concrete structures, humidity model applied with random stochastic process will be demonstrated and approximately 1.5 million nodes and three million elements will be used. Furthermore, multi-boundary conditions to describe the actual chloride concentrations and humidity distribution will be applied instead of only constant boundary conditions.

2. Basic Diffusion Formulation of Unsaturated Concrete

2.1. Governing Equation

The two fully-coupled partial differential equations governing the coupled chloride and moisture diffusion through non-saturated concrete are simply derived by employing the mass balance equations and Fick’s law [14,15].

First, the flux of chloride ions (J) through a unit area of porous media depends on the gradient of chloride ions as well as the gradient of moisture as. Thus, Fick’s law is modified in this study to include the coupling effects between moisture and chloride transfer. The governing equations are shown in Equations (1) and (2). where D = chloride diffusion coefficient (cm2/day); C = free chloride concentration (in gram of free chloride per gram of concrete, g/g); D = humidity diffusion coefficient; H = pore relative humidity; = humidity gradient coefficient, which represents the coupling effect of moisture diffusion on chloride penetration; and = chloride gradient coefficient, which represents the coupling effect of chloride ions on moisture diffusion.

When chloride ions ingress into the concrete, some of them are bound to the internal surface of the cement paste and aggregates, which are called bound chlorides and freely go through the concrete. Steel corrosion is related only to the free chloride content but not to the total chloride content [14].

The mass balance of chloride ions and moisture can be expressed using Fick’s second law as Equations (3) and (4), where

: total chloride concentration (in gram of free chloride per gram of concrete, g/g);

: water content;

: chloride binding capacity; and

: moisture binding capacity.

Equations (3) and (4) can be rewritten as, where is coupling parameter, , and is coupling parameter, .

The general boundary conditions are as below, where = convective chloride coefficient, = convective relative humidity coefficient, = ambient chloride ions, and = ambient relative humidity.

and are the part of boundary with constant chloride ions and relative humidity and and are the part of boundary subjected to specified chloride ions and relative humidity flux, respectively. and form the complete boundary surface for the chloride diffusion problem, and and form the moisture diffusion problem [16].

2.2. Material Parameters

To numerically solve the chloride diffusion problem, many material parameters must be determined: the moisture capacity (), chloride binding capacity (), humidity diffusion coefficient (), and chloride diffusion coefficient () (in Equations (5) and (6)) [15].

2.2.1. Moisture Capacity

The moisture capacity of the concrete was developed based on the multiphase and multiscale model [17]. Assuming that the effect of the shrinkage of concrete can be evaluated simply by the average of the moisture capacities of the aggregate and the cement paste as shown in Equation (11), where and = weight percentages of the aggregate and cement paste, and and = moisture capacities of aggregate and cement paste, which can be calculated based on the model developed (Xi et al., 1994a, b) and Xi (1995a, b) [18,19,20,21].

2.2.2. Chloride Binding Capacity

A modified relationship between the bound chloride , and the free chloride , was established by Tang and Nilson (1993) based on Freundlich isotherm, and was proposed by Xi and Bazant (1999); which is differentiated with respect to and then one would obtain the chloride binding capacity of concrete as defined by Xi and Bazant (1999) to be; where A and B = chloride adsorption related constants, 0.3788 and 1.14, β = ratio of pore solution volume to concrete weight, L/g, and = weight ratio of C-S-H gel to concrete (g/g) [14,22].

Based on the definition of β, the following Equation (14) can be easily derived; where and = weight percentages of cement paste and aggregate in concrete mix (g/g), = cement paste adsorption isotherm, = aggregate adsorption isotherm, and = density of pore solution measured in g/L.

Specific gravities of concrete and C-S-H are similar, therefore, Xi and Bazant (1999) assumed that the weight fraction of C-S-H in concrete, , is equal to the volume fraction of C-S-H in concrete, , then, where = volume fraction of anhydrous pores of cement particles, and = volume fraction of capillary pores of cement paste [14].

2.2.3. Humidity Diffusion Coefficient

The humidity diffusion coefficient of concrete depends on the diffusion coefficients of aggregate and cement paste. Using the composite theory (Christensen, 1979), the effective diffusion coefficient of concrete can be evaluated as following Equation (16), where = aggregate volume fraction, = humidity diffusion coefficient of the cement paste, and = humidity diffusion coefficient of the aggregates [23].

The humidity diffusion coefficient of aggregates in concrete is very small due to the fact that the pores in aggregates are discontinuous and enveloped by cement paste and so it can be neglected. The humidity diffusion coefficient of cement paste can be predicted by using the empirical formula described in Xi et al. (1994b) [19].

2.2.4. Chloride Diffusion Coefficient

Chloride diffusion coefficient in saturated concrete was studied by Xi and Bazant (1999) as following Equation (17),

In Equation (17), the functions of the chloride diffusivity consist of the concrete curing time (), gravel volume fraction (), relative humidity (H), temperature (T), and free chloride concentration () [14].

First factor of chloride diffusion coefficient accounts for the effect of the water–cement ratio () and curing time () in Equation (18),

The second influence factor is to consider the effect of composite action of the aggregate and the cement paste in Equation (19), where = volume fraction of aggregate in concrete, and and = chloride diffusion coefficient of aggregate and cement paste.

The third factor, f3(H), in Equation (20) is to consider the effect of relative humidity level on the chloride diffusion coefficient. A model proposed by Bazant et al. (1972) can be used, which was developed initially for moisture diffusion [24]. where = critical humidity level, 0.75.

Arrhenius law was used by Xi and Bazant (1999) to introduce the temperature effect of the forth factor in chloride diffusion coefficient as shown in Equation (21), where

and = reference and current temperatures in Kelvin, T0 = 296 K;

R = gas constant, 8.314 J/mol·K; and

U = diffusion process activation energy, depending on w/c ratio [14].

The detail description of diffusion process activation energy, U can be found in the paper of Ababneh et al. (2003) [15].

The dependence of chloride diffusion coefficient on free chloride concentration for the fifth factor is presented in Equation (22); where and m = 8.333 and 0.5, respectively, according to Xi and Bazant (1999) [14].

3. Finite Element Formulation

In order to solve a time-dependent coupled moisture-chloride diffusion problem, a large linear system equation was derived with finite element method. The finite element formulation will be briefly introduced in this section. The continuous variables in the coupled chloride and moisture diffusion equations, free chloride () and relative humidity () are spatially discretized over the space domain, . The domain discretization can be described as shown in Equation (23). in which nel is the total number of elements in space domain and is an element. It is also defined as the boundary of computational domain and the boundary of subdomain.

The unknown variables in Equations (24) and (25) are defined in terms of nodal values, and , where is the triangle element shape function. The notations and are row and column vectors, respectively. The element shape functions are expressed as following Equation (26), in which N is the shape function for node i and n is the total numbers of nodes in an element. The unknown vectors of free chloride in Equation (27) and relative humidity in Equation (28) can be defined as,

The nodal free chloride concentrations and relative humidity are solved by substituting the approximated values of Equations (24) and (25) into governing equations of Equations (5) and (6), and applying the Galerkin procedure to the weak forms, then the finite element matrix can be obtained as shown in Equation (29): where the element matrices and vector are as following Equations (30)–(32),

In detail, the components in element matrices are as Equations (33)–(38),

Finally, Equation (39) is also discretized in time space with time interval as following,

The value of parameter is related to the solution method adopted in the program. Typical values of are 0, 1/2 and 1 correspond to fully explicit, semi-implicit and fully implicit methods, respectively. The semi-implicit method called Crank–Nicholson method is used in this study.

Equation (39) is simplified as linear system equation. where

4. Implementation of Parallel Finite Element Method

4.1. Various Programs Adapted in Parallel Finite Element Program

In order to implement the parallel finite element program for the coupled moisture and chloride problem, various programs were employed such as Triangle for mesh generation, ParMETIS, PETSc (Portable, Extensible Toolkit for Scientific Computation), and MPI (Message Passing Interface) [25,26,27,28,29,30,31].

In this study, Triangle was used for the mesh generation of triangle element, which was created at Carnegie Mellon University. Triangle generates exact delaunay triangulations, and are suitable for finite element analysis [25].

PETSc (3.0.0 p8) is a large and versatile package integrating distributed vectors, distributed matrices in several sparse storage formats, Krylov subspace methods, preconditioners, and Newton-like nonlinear methods with built-in trust region or linesearch strategies and continuation for robustness. It is designed to provide the numerical infrastructure for application codes involving the implicit numerical solution of PDEs, and it uses on MPI for portability to most parallel machines. The PETSc library is written in C, but may be accessed from user codes written in C, Fortran, and C++ [32].

MPI (Message Passing Interface) is a standardized and portable message-passing system designed to function on a wide variety of parallel computers. The standard defines the syntax and semantics of library routines and allows users to write portable programs in the main scientific programming languages (Fortran, C, or C++) [27,28,29].

ParMETIS extends the functionality of METIS and includes routines based on a parallel graph-partitioning algorithm that are especially suited for parallel computations and large-scale numerical simulations involving unstructured meshes. In typical FEM computations, ParMETIS dramatically reduces the time spent in interprocess communication by computing mesh decompositions such that the number of interface nodes/elements is minimized [30].

4.2. Overlapping Domain Decomposition Method with Additive Schwarz Preconditioner

A parallel program is typically developed by dividing the program into multiple fragments that can execute simultaneously, each on its own processor. In the finite element analysis, this can be accomplished by applying a domain decomposition method. Domain decomposition method is the method usually used for solving large scale system equations and it is also suitable for parallel programming because of data locality [15]. There are two types of domain decomposition methods, overlapping and non-overlapping methods. In this study, the overlapping method was employed to solve the linear sparse matrix. That is because the advantage of overlapping domain decomposition is easier to setup in algebraic approach and faster convergence than non-overlapping domain decomposition. Furthermore, the boundaries of extended subdomains are smoother than non-overlapping subdomains.

Iterative solver must be used in the iterative domain decomposition method. In this study, for the iterative solver, GMRES (Generalized Minimal Residual method) was chosen for both global and local matrix. GMRES is mainly chosen because of its ability to solve non-symmetric linear system as in the case of our problem. To improve the convergence of this problem, the additive Schwarz method preconditioner was applied. Figure 1 shows the flow chart of parallel pre-process and FE solver.

Figure 1

Framework of parallel FE method based on PETSc: (a) flow chart for parallel pre-process; and (b) flow chart for parallel FE solver.

5. Numerical Results

5.1. Applied Bridge Overview

The Castlewood Canyon bridge is located on Highway 83 in the Black Forest of central Colorado. The original two-lane reinforced concrete arch bridge was built in 1946. The arches are 1.93 m wide by 1.78 m deep at the base with the depth tapering down to 1.0 m at the highest point. This bridge was severely dilapidated and was in need of repairing, enlarging, and strengthening. Parts of the concrete had spalled off of the deck, columns, and arches, and the steel rebar under concrete were severely rusted. In 2003, the original arch was repaired, the spandrel columns and decks were replaced and widened from about 0.9 m to 1.0 m including railings, and the overall length of the bridge was increased from 114 m to 123 m as shown in Figure 2.

Figure 2

Castlewood Canyon bridge, Franktown, Colorado.

5.2. Parallel Finite Method of Large-Scale Concrete Structure

5.2.1. Modeling of Castle Wood Canyon Bridge with Large Meshes

In order to analyze the penetration of chloride and humidity in a concrete bridge, a large number of meshes are needed to capture the diffusion phenomenon within thin layer from concrete top surface to steel rebar. For depicting in detail concrete cover depth on rebar, the information of nodes and elements are created by mesh generator. Triangle as mesh generator is specialized for creating two-dimensional finite element meshes.

For castle wood canyon bridge, approximately 1.5 million nodes and three million elements were generated for input files as shown in Figure 3. To visualize and check the mesh size and shape, ParaView was employed as illustrated in Figure 4 [33]. ParaView is an open-source, multi-platform data analysis and visualization application. It is developed to analyze extremely large datasets using distributed memory computing resources and can be run on supercomputers to analyze datasets of terascale as well as on laptops for smaller data. For large scale information, VTK file format is adopted because this format is easy to read and write by hand or programmatically. This file format is automatically created when the program runs. For partitioning origin meshes, Parmetis was used to be embedded in the parallel program to automatically partition. Parmetis is an MPI-based parallel library that implements a variety of algorithms for partitioning unstructured graphs, meshes, and for computing fill-reducing orderings of sparse matrices. ParMETIS extends the functionality provided by METIS and includes routines that are especially suited for parallel AMR computations and large scale numerical simulations. The algorithms implemented in ParMETIS are based on the parallel multilevel k-way graph-partitioning, adaptive repartitioning, and parallel multi-constrained partitioning schemes developed.

Figure 3

Modelling of Castlewood Canyon bridge.

Figure 4

Partitioning of Castlewood Canyon bridge with 16 processors.

5.2.2. Environmental Humidity Model

Environmental humidity model plays an important role of diffusion of chloride in concrete. The humidity model used in this study was proposed by Bazant and Xi (1993) [13]. Even though the moisture and temperature transport in concrete were coupled, humidity fluctuation was considered and temperature was assumed to be constant during analysis time. To establish a practicable and realistic humidity input model, real climate record in Chicago was employed as shown in Figure 5.

Figure 5

Environmental humidity Record, midway station, Chicago (Bazant et al., 1993) [13].

The environmental humidity model consists of three different components as shown in Equation (43) and demonstrated in Figure 6. First, is the mean value, standing for the stable horizontal trend. Second, is a random phase process corresponding to the harmonic variation of humidity. Third, is a random normal distribution with a specific mean and variance [13]. where , mean value of environmental model, , the random phase process of a one day period with = 0.1 and the random phase process of a one year period with = 0.08, and have uniform distributions with constant density , , normally distributed random numbers are generated with specific mean and variance, 0.0046 [13].

Figure 6

Environmental humidity model: (a) humidity model without random noise; and (b) humidity model with three components (including random noise).

For the chloride concentration on the top surface of concrete, various boundary conditions were used to actualize the realistic boundary conditions as illustrated in Figure 7. One can see that the interference between the interfaces of different concentration.

Figure 7

Comparison of single and various boundary conditions: (a) single boundary condition; and (b) two boundary conditions.

5.2.3. Speed-Up for Parallel Algorithm

Speed-up is measured for the performance of the parallel implementation of finite element program. The definition of speed-up is a process for increasing in performance between two systems processing the same problem. With regard to speed-up in this study, 8–2048 processors were used in this study and the meshes with the number of about 1.5 million nodes were analyzed. When using up to eight processors, the program was not operated due to memory capacity when input data can be assigned on each memory which means the problem was too huge to solve with a couple of processors.

As the number of processors used in the analysis increased, speed-up was improved up to 512 processors because speed-up was better than ideal condition. After more than 512 processors were used, speed-up gradually increased because the communication time between the processors increased. As shown in Figure 8, one can see the optimum point of number of processors.

Figure 8

Speed up over number of processors.

5.2.4. Effect of Boundary Condition

(1) Constant boundary condition

The physical and material models of the castle wood canyon bridge used in this paper were mentioned above. The validation of these models was secured by Ababneh et al. and Suwito et al. [14,15]. Ababneh et al. proved the model with alternating-direction implicit (ADI) finite-difference method and the numerical solutions were compared with the experimental results obtained by the 90-day ponding test (AASHTO T 259-80) [14]. Suwito et al. verified the accuracy of the model implemented with finite element method and numerical results agreed well with the experimental data [15].

It was assumed that the concrete bridge contained initially no chloride ions and has 60% relative humidity (RH). The top surface of concrete bridge was exposed to 3% NaCl and 100% RH. To analyze the diffusion of coupled chloride and humidity, 8 to 2048 processors and approximately 1.5 million nodes were used with parallel finite element method.

Obviously, this example is not to simulate the concrete behavior under service condition. The triangle element meshes were employed and the whole domain was partitioned into 16 sub-domains to the same as the number of processors shown in Figure 3 and Figure 4. The sub-domains were divided with ParMetis so that each sub-mesh had almost the same number of elements.

Figure 9 and Figure 10 show the contour of chloride and humidity distribution into the entire domain. The total and free chloride concentration at 50 mm below from the top surface of concrete is Figure 11 and the variation of relative humidity is Figure 12 over time up to 550 days. As a result of analysis, the concentration of chloride and humidity is continuously accumulated inside concrete. Specially, the variation of humidity inside concrete might not be increased in reality due to the seasonally change of the humidity distribution, so that periodic humidity model should be applied as the boundary condition.

Figure 9

Contour of chloride distribution.

Figure 10

Contour of humidity distribution.

Figure 11

Chloride concentration over time at 50 mm depth.

Figure 12

Humidity over time at 50 mm depth.

(2) Periodic boundary condition

In order to simulate the real diffusion phenomenon in concrete, the periodic boundary condition should be required. Figure 13 and Figure 14 describe the result of chloride concentration and relative humidity at different depth over time. As shown in Figure 13, total concentration is gradually increasing until almost one year. That means the chloride concentration can be affected by the periodic humidity boundary condition. After one year it suddenly increases because of the influence of the constant chloride initial condition. However, the total chloride concentration at 50 mm is reduced up to 15% rather than the concentration using constant boundary. In order to predict and demonstrate the realistic and seasonal chloride ingress through concrete cover, seasonal change information of chloride concentration on the top surface of concrete deck.

Figure 13

Chloride concentration over time at 50 mm depth.

Figure 14

Humidity over time according to depth.

In Figure 13, the change of humidity distribution is influenced from the humidity model inputted as boundary condition. The short period random noise, one of random noise terms in environmental humidity model, can only reach a shallow portion near the top surface of concrete deck. As the concrete depth increased, the effect of short period random noise decreases. However, the seasonal variation of humidity concentration has an influence on 50 mm depth of concrete deck. At 50 mm from the top surface of concrete deck, the maximum value of humidity is about 77%, which is decreased up to 20% comparing with constant boundary condition. In addition, the result of about 500 days shows that the humidity trend is reversed between 10 mm and 50 mm, which means that evaporation of humidity just near concrete surface is easier than evaporation at 50 mm. Therefore, this prediction model is reasonable and effective to describe the diffusion phenomenon of chloride and humidity.

The parallel finite element program was developed based on the robust mathematical material model. The program can be used to simulate the coupled moisture and chloride penetration into non-saturated concrete structures. Chloride ion is one of sources causing the steel corrosion in reinforced concrete structures. The material parameters related to chloride and moisture diffusion in concrete are taken into account. These parameters include chloride and moisture diffusion coefficients, moisture capacity, chloride binding capacity, and the coupling parameters reflecting the coupling effects between moisture and chloride transfer in concrete.

For the implementation of parallel FE analysis, Triangle for mesh generation, ParMetis, PETSc, and MPI are employed. This program also used the overlapping domain decomposition method with additive Schwarz preconditioner. As the result of simulation, the computation time decreased until the number of processors became optimized when the number of processors increased. Then, the computational time increased because the communication between the processors increased.

The present model can be used to simulate the unsaturated concrete structures subjected to other aggressive chemicals from de-icing salt. The framework of present model can be extended to simulate the multi-species de-icing salts ingress into non-saturated concrete structures in future work.