Simulation of Chloride Ingress through Surface-Coated Concrete during Migration Test Using Finite-Difference and Finite-Element Method

Polymer surface coatings are commonly used to enhance the corrosion resistance of concrete structures in saline environments to ionic diffusivity; this diffusivity can be determined by migration tests. This paper presents the simulation of the effects of the surface coatings on migration tests by solving the Nernst-Planck/Poisson model using both finite-difference method and finiteelementmethod.These two numericalmethods were compared in terms of their accuracy and computational speed.The simulation results indicate that the shapes of ionic profiles after migration tests depend on the effectiveness of surface coatings.This is because highly effective surface coatings can cause a high ionic concentration at the interface between coating and concrete.The simulation results were also compared to homogenized cases where a homogenized diffusion coefficient is employed.The result shows that the homogenized diffusion coefficient cannot represent the diffusivity of the surface-coated concrete.


Introduction
The corrosion of the reinforcing steel bars by chloride ingress is a significant problem for the marine concrete structures exposed to seawater.To delay this chloride ingress, the polymer surface coating has been applied to the reinforced concrete.Previous literature clearly shows that this surface coating improves the durability of reinforced concrete (RC) structures [1][2][3].However, the polymer surface coating makes the evaluation of the service life of RC structures difficult, because the diffusivity of the surface-coated concrete as a composite material becomes more complicate than bulk concrete [4].Therefore, the simulation of chloride ingress through surface-treated concrete requires further research.
Migration tests [5,6] provide a convenient way to estimate diffusivity of cement-based systems.The classical mathematical treatment to analyze the experimental data of the migration test solves the movement of single species and approximates the electric field by a constant equal to the electric potential difference between intervals of the sample depth divided by the sample length.Recently, Glasser et al. [7][8][9] used the Nernst-Planck/Poisson (NPP) model to analyze data of the migration tests [5].This approach improves upon the classical model by including multi-ionic interaction, distribution of electric potentials, and evolution of the electric field.The application of the NPP model explains why the diffusion coefficients of identical specimens calculated from the classical model depend on the different boundary conditions.Also, a slight difference in potential profiles leads to a large difference in chloride concentration profiles, so that the assumption of a constant electric field causes significant errors in the distribution of chloride across the entire sample [5].Furthermore, even though most cases of ionic ingress in reinforced concrete structures dominantly occur by diffusion, in aggressive environments the ingress of multiple ions, such as chloride ingress, calcium leaching, and sulfate attack, can occur in concrete simultaneously.For this reason, the NPP model including the migration term can describe the interaction between different ionic species and accurately predict the distribution of multiple ions in reinforced concrete structures.
The NPP model is considerably more complicated than the classical single-species model.For the non-steady-state and surface-coated cases, a numerical scheme is necessary to compute the NPP model.As early as 1995, Kato [10] proposed the finite-difference method (FDM) to solve the NPP model for steady-state cases; however the finite-element method (FEM) has been recently employed preferably for the NPP model [5,6,8].The NPP model was further improved by taking into account unsaturated conditions [11,12] and binding isotherms [13].The analytical solutions of the NPP model for the steady-state response and FEM solutions for the non-steady-state condition were presented by Rani et al. to simulate the migration test.Recently, it was reported that FDM computations achieve 2X to 3X speed-up in parallelization of FDM/FEM computations for large mesh sizes [14].In fact, the surface coating requires increasing the mesh size because the surface coating is a thin section.Therefore, FDM takes advantage of computational speed to simulate the surface-coated concrete.
The present study evaluates the effect of surface coatings on the migration test using these numerical simulations.Since FEM and FDM have different advantages in accuracy and computational speed, it is necessary to discuss the strength and limitations of the methods for simulating the effect of surface coatings.Furthermore, the present study attempts to simulate the effects of the surface coatings by solving the NPP model during migration tests.The simulation results are discussed, and the different shapes of ionic profiles produced by different effectiveness of the surface coatings are compared.In addition, homogenized cases employing a homogenized diffusion coefficient are presented and compared to discretely simulated ones having two diffusion coefficients of surface coating and concrete.

Model Description
In general, the range of pore size in hardened concrete varies from a nanometer to a millimeter.Some pores connect to others, forming continuous paths for ionic transport in a saturated system.These diffusion paths in concrete and a surface coating are tortuous compared to those in free liquid.For this reason, the term "tortuosity" was introduced to account for this complex pore structure [15,16].The diffusion coefficient of each ionic species has the following relation with tortuosity [15]: where  is the tortuosity of the saturated material and    is the diffusion coefficient of the species in free water.Different tortuosity values were used to distinguish between concrete and surface coating to account for their different pore structure.An interfacial resistance may develop when two materials having different pore structures are connected.Zhang et al. [17] found no significant evidence of a strong interfacial resistance between surface treatment and cement-based materials; therefore interfacial resistance was not included in the present mathematical model.The ionic profiles of the simulations are presented in terms of ionic concentration of the pore solution to estimate possibility of chemical reactions at the interfacial region.
Migration tests are commonly performed in saturated conditions and constant temperature.Most studies that simulate the migration test assume that chemical reactions are negligible [7].However, recent research shows that this assumption can overestimate the tortuosity of the material [18,19].The objective of the present study is to show the trend of ionic behavior in the surface-coated concrete during migration tests so that the chemical reactions are not included.
2.1.Nernst-Planck/Poisson (NPP) Model.The NPP model involves a separate mass balance equation for each of the ionic species present.The flux of an ionic species  in solution is given by where   is the concentration,   is the diffusion coefficient,   is the valence number of the species,  is the Faraday constant,  is the perfect gas constant,  is the temperature, and  is the electrical potential.The mass conservation for the species is expressed as Substituting ( 2) into (3) results in the complete Nernst-Planck equation: The Nernst-Planck equation must be solved simultaneously with the Poisson equation that directly relates the electric potential  to the electric charge: where  is the total number of ionic species,  is the absolute permittivity, and  is the fixed charge density.Here, we assume that the fixed charge density is negligible because the fixed charge density is not a relevant parameter for most porous materials [8].

Numerical Schemes
The finite-difference procedure for the surface-coated concrete is newly developed here.The mathematical treatment of the FEM follows the work of Yoon et al. [4], but here the emphasis is on developing a strategy in how to include the coating and have an efficient simulation.
3.1.Finite-Difference Procedure.For many cases of practical relevance, the diffusion coefficient of NPP model is assumed to be a constant value, which can then be extracted from the parenthesis.If the object consists of two different materials, it must be dealt with differently.As the first step, the governing equation of the unidirectional case is expanded as Equation ( 6) is derived by applying the product law of the differential equation.However, for the simulation of surface coatings,   / is problematic because of the discrete (step) function of   () along the -axis.The diffusion coefficient function cannot be directly differentiated by .For this reason, it must be converted to a differentiable form.Since the original step function reveals no change except at the interface, the differentiated function is equal to zero everywhere but the interface must have infinite increase.This is the Dirac-delta function.Instead of the differentiated function (  /), this study used an approximate delta function, which is also known as normal distribution function.
The original step function can be also reexpressed by the integration of the approximate delta function.This simpler expression is similar to the Heaviside step function and suggested in (8).Then ( 7) and ( 8) are plotted in Figure 1 for different values of Depending on different values of , solutions of the equation of the diffusion coefficient may be varied in a complete set.A simple set of the NPP model was analytically solved to compare with the exact solution.This example considers steady-state conditions with no electrical potential.The following assumptions are also made: (1) the dimensionless concentration ( C) at  = 0 is 1; (2) the interface between the two materials is located at the center of the dimensionless depth; (3) the diffusion coefficient of the one material is twice that of the other.Figure 2 and Table 1 show the concentration profiles of the exact solution and the analytical results from the approximate delta function with errors of concentrations for different values of /.As the value of / increases, the result from the approximate delta function is close to the exact solution.
From ( 7) and ( 8), ( 6) can be modified as The partial differential terms must be substituted with difference quotients induced by a Taylor series, and the depth and the time are also divided by each interval  and  to form the lattice.The temporal derivatives can be taken with a Crank-Nicolson theta differencing scheme (i.e.,  = 0.5), and the second spatial derivative is thus From these equations, a linear set of equations can be solved for each time step.The concentration at a new time step depends on the previous time step and can be solved explicitly.First, the Poisson equation is solved by using the initial concentrations of species, and then the gradient potentials are calculated on the -lattice points.The solutions of the Nernst-Planck equation at the next time step are obtained for the species  using the calculated electric potential and (9).A loop is performed on all the ionic species, and then the Poisson equation of the next time step is solved with new concentrations of species.
3.2.Finite-Element Procedure.Yoon et al. [4] proposed using the finite-element method (FEM) to solve the NPP model, suggesting the use of uncoupled and coupled algorithms.In the calculation, both algorithms provided the same results.The presence of surface coating requires large mesh size near the boundaries and interfaces to avoid oscillations.
Because the uncoupled algorithm deals with smaller sized matrices and the coupled one requires a classical Newton-Raphson method to solve the nonlinear set of equations for each time step, the computation time is reduced by using the uncoupled algorithm to simulate the case with the surface coating.The detailed computation of the coupled and uncoupled algorithm is described in Appendix.The Nernst-Planck equation is applied to the weak form of the Galerkin method by using the virtual concentration (  ) over the domain (Ω): By applying the divergence theorem, the weak form is expressed as This weak form is discretized using the explicit scheme ( = 0) and the shape function as where [] is the shape function and [] is the derivative of the shape function.The diffusion coefficient matrix ([]) does not need to be differentiated in this method, because the divergence theorem makes it independent of the spatial derivative.Hence the diffusion coefficient matrix is organized as a diagonal matrix, and the values of the matrix correspond to the diffusion coefficient of each spatial element.To solve the concentration of the next time step, the equation is expressed as where [] is an identity matrix equal in size to the diffusion coefficient matrix.The Poisson equation is also computed in the same way; its discretized weak form is given by The electric potentials with new concentrations can be obtained by This scheme is based on the simulation performed by Samson et al. [8].The difference in this simulation of the surfacecoated concrete is the diffusion coefficient matrix, which depends on the coating and concrete area.

Simulation Results
A case of the surface coating is solved using numerical and analytical methods (see Table 2).The results are compared to discuss their accuracy and number of required minimum iterations.To compare the numerical results to the exact solutions, the effects of the electric field are excluded so that  the exact solution can be obtained using Fick's law.The exact solution suggested by Crank [20] is as follows: where  is √ 1 / 2 and  is (1 − )/(1 + ).The explicit scheme ( = 0) is used for all numerical methods, and all numerical calculations are performed with 100 elements.As shown in Figure 3, since the analytical solution from the approximate delta function is closed to the exact solution when the value of  is /100, this value is applied to FDM computation.Depending on the chosen shape function of FEM, their results can be different.Thus, both the classical linear and quadratic shape functions are presented.As shown in Figure 3, all numerical methods give good approximations compared to the exact solution.To discuss their accuracy in detail, their errors are compared in Figure 4(a).The degree of error is similar before the interface () between the coated and concrete regions and then increases after that.The FDM shows the highest increase of error at the interface, thereby reducing its total accuracy.After 1.5 mm, when the exact concentrations are small (∼0.001M), the errors in both methods rapidly increase due to the small dominators.Thus, although errors are large after 1.5 mm, differences between the exact and numerical solutions are small in Figure 3. Figure 4(b) shows the number of minimum required time steps of different numerical methods.The computation speed using FDM is three and sixteen times faster than that in FEM, respectively, with linear and quadratic shape functions.Therefore, the FDM improves the computation speed and reduces the accuracy compared to FEM.

Discussion of Effects of the Surface Coating on the Migration Test
K. Krabbenhøft and J. Krabbenhøft [21] carried out a simulation study of migration tests in concrete without surface coatings.They applied the NPP model to the migration test using several boundary and initial conditions.In the present study, one of their conditions was used but included 3 mm of surface coating as an additional boundary condition.Table 3 lists the boundary and initial conditions established for the simulations; chemical reactions were neglected.This hypothesis is justified by the fact that non-steady-state migration tests have a short experimental duration (1∼2 days).In addition, neglecting chemical reactions implies that there is no change to the microstructure of the samples during the migration tests [5], which is equivalent to assuming that the diffusivities remain constant.Specimens are assumed to be cured in calcium hydroxide [22] or sodium hydroxide [23] solutions to obtain the saturated condition.Figures 5, 6, and 7, respectively, demonstrate how  concrete / coating ratios of 10, 5, and 2 affect the results of the migration tests using both FDM and FEM.Each value  of  concrete / coating means that surface coatings acting as a barrier are 10, 5, and 2 times more effective than concrete in preventing the penetration of ions.Hence,  concrete / coating of 10, 5, and 2, respectively, represent large, moderate, and small differences in diffusivities of coatings and concrete.The simulation results of the migration tests obtained from FEM and FDM are almost identical in Figures 5, 6, and 7.This is because the test duration is much shorter than the service life of a reinforced concrete structure.The common migration tests were performed between 12 hours and 2 days, but the service life of a reinforced concrete structure is decades.The error by numerical methods increases over time.
Therefore, although FEM and FDM have different accuracy, they can provide similar accuracy when computing short timeframes of the migration tests.In the simulation results, negatively charged ions, such as chloride and hydroxide, tend to move toward a positively charged cell on the right (see Figure 5), whereas positive charged ions, such as sodium and potassium, tend to move in the opposite direction due to the externally applied voltage during the migration test.In the case of  concrete / coating = 10, sodium ions having a positive charge are supplied from the right cell, but their movements to the left are blocked by the surface coating, accumulating at the interface between the surface coating and concrete bulk.Consequently, chloride and hydroxide concentrations at the interface increase over their concentrations at boundaries to maintain the electroneutrality.In the cases where  concrete / coating = 5 and 2, no abnormal behavior where the concentrations of ions at interface exceed their concentrations at boundaries is observed.High concentrations of species increase the possibility of chemical reactions.From these results, it seems clear that substantial differences in diffusivities of surface coatings and concrete bulk increase the potential for unexpected chemical reactions, demonstrating that the  concrete / coating ratio is an important parameter to evaluate chemical reactions and constant diffusivities during the migrations tests.
Both the FDM and FEM provide similar results regardless of the  concrete / coating ratios.In order to compare discrete diffusivities and homogenized effective diffusivity, the diffusion coefficients are integrated over the Representative Elementary Volume (REV) in the homogenization technique [24,25].The volumetric coefficient average is given by where ⟨⟩ is the effective diffusion coefficient.Dashed lines in Figures 5-7 are the results of the simulation using these effective diffusion coefficients.Exceptionally, Figure 7 shows that the homogenized case matches well with the discrete one.This shows that only when the difference in diffusivities of coating and concrete is small, the homogenized diffusion coefficient is reasonable; otherwise it leads the effective diffusivity to be erroneous.It is noteworthy that the solid and dashed lines in Figures 5 and 6 are unmatched.

Conclusions
The present study newly introduced the FDM procedure to solve NPP model of the surface-coated concrete, and further conclusions are as follows: (i) The effects of the surface coatings can be simulated using FEM and FDM.Without considering the effect of the electric field, the comparison between both numerical schemes demonstrates that FDM is characterized by high computational speed and low accuracy, whereas FEM is characterized by low computational speed and high accuracy.However, both FEM and FDM have a similar accuracy in simulations with short timeframes of the migration tests.(ii) Highly effective surface coatings act as a barrier.This barrier action is more effective when the difference in diffusivities of coating and concrete is more substantial.This is because concentrations of all species increase at the interface, meaning increase in the possibility of unexpected chemical reactions.(iii) Instead of two diffusivities of coating and concrete, a homogenized diffusion coefficient results in unmatched diffused depths, except in the case where diffusivities of coating and concrete are similar to each other.Therefore, a homogenized diffusion coefficient cannot represent the diffusivity of the surface-coated concrete.
International Journal of Polymer Science   with a low value of  2 depends on the number of elements.Thus, both algorithms are applicable to any cases.

Figure 2 :
Figure 2: Exact solution and concentration profiles solved from the approximate delta function.

Figure 3 :
Figure 3: Comparison between an exact solution and numerical solutions by computing a coated material.

Figure 4 :
Figure 4: (a) Errors of numerical solutions.(b) Required minimum number of time steps.

Table 1 :
Comparison of solutions from a simple set of NPP model.

Table 2 :
Setup used for the coating effect simulation of Fick's law.

Table 3 :
Boundary and initial conditions for simulated experiments.