Numerical model of time-dependent diffusion of chlorides in the concrete based on 2D four-node isoparametric element

The enhancement of 2D model for the diffusion of chloride ions into reinforced concrete considering the time-dependent diffusion coefficient is discussed in the article. The non-stationary Finite Element Model is based on the four-node isoparametric element. The algorithm is implemented in MatLab compatible environment and an important part of source code is presented. The results of Finite Element Analysis are compared with the results of 2D analytical diffusion model.


Introduction
One of the most significant problems in case of many bridge decks not only in Czech republic is the corrosion of reinforcing steel because of the ingress of chloride salts applied to melt snow and ice [1][2][3][4]. The chloride ions penetration process is dominantly governed by diffusion in case of bridge decks. Of course, the combined effect of carbonation, chlorides ingress, and mechanical load must be taken into account [5,6]. However, the introduced model considers simply the chloride diffusion, which in some cases may be a major source of deteriorations in the concrete structures.
The article directly follows the development of in house code for the modelling of chloride ion ingress. The original algorithm contains a three-node element that has certain limitations that could cause inaccuracies of a specific input parameter configuration [7]. Also, the possibilities of different protection of the reinforcement or salting-down have been introduced already [8,9]. The finite element model (FEM) is based on the newly derived four-node isoparametric diffusion element that has not been published yet in this form, up to author knowledge. Similar elements are used for solving problems of statics and dynamics, or heat [10,11], but so far they have not been introduced for the process of diffusion of chloride ions in concrete. The following text is a brief overview of the methodology, derivation of the finite element, and an example with comparison of results to analytical solution of the 2D diffusion. It should be noted that the element and the whole model are based on temperature-diffusion analogy. The combination with nonlinear calculation is a suitable tool for the possible occurrence of corrosion of the reinforcement in the concrete structure.

Methodology
Since the diffusion process is governed by the same differential equation as heat transfer, the transient thermal analysis is used for the solution of the mass diffusion problem. Chloride ion diffusion is a non-stationary transport process in which chloride ions penetrate until the steady state in the material is reached.
Thus non-stationary problem is of interest in general cases The task can be described by a system of partial differential equations, which can then be solved by FEM [2,7,[12][13][14] or analytical solution [15,16]. The second Fick's law is used to describe the change in concentration in time and place [15].
The problem of chloride diffusion can be solved by probabilistic assessment (see eg. [1,17]) due to the large scatter of the input parameters. Probabilistic applications (such as [18][19][20][21]) create a demand for large computing capacity. In this case, the model is implemented using the scripting language in the commercial Finite Element Analysis (FEA) package [1,22] and runs relatively slowly. Therefore, creation of the executable code is a process that could speed up the simulation.
The research was based on the code presented in the previous works [7,23], where a three-node isoparametric element was used. This element brings several advantages, such as faster calculation, more accurate boundary conditions for irregular shapes of the structure etc. On this basis, the calculations in the publications [7-9, 24, 25] were presented, where both deterministic and probabilistic procedures were presented.
Conversely, for advanced modelling of a possible crack in a concrete cover, it is desirable to evaluate a four-node element that is less sensitive to aspect ratio and its accuracy may be higher in specific tasks.

Deriving the finite element
The process of derivation and analysis of the finite isoparametric element was prepared on the unitary shape in the coordinate system η, ξ. Using the selected shape functions, the reference element can be displayed on the actual element in the x, z coordinate system. There are several types of isoparametric elements; examples are given in Fig 1. It is a linear triangular element with three nodes (LT3) and bilinear quadrilateral element with four nodes (BC4) [26]. A four-node isoparametric element (BC4) was chosen for further assessment. Four shape functions are required for this element according to the following equation (1): The individual functions N1, N2, N3 and N4 represent the specific corners of the fournode reference element. As an example, the graphical representation of one of the functions is shown in Fig. 2, where the value above the reference point reaches 1 and at the other points the value is 0. The derivation of the functions for the diffusion problem follows. In this case, the chloride ions concentration in the individual nodes are the resulting parameters. The heatdiffusion analogy is used, and the non-stationary heat process is converted on the diffusion problem. At first, it is necessary to express the concentrations C0 by means of approximation functions: where Ni are the functions described in formulas (1) and Ci are the searched concentrations in each node. The expression (2) can be written as: Subsequently, it is possible to use an analogy with a deformation problem that contains geometric equations [27,28] and derivations of the equations: (4) where Dc is the concrete diffusion coefficient [m 2 .s -1 ].
The following vector has been obtained by a combination of the formulas (3) and (4): The individual terms shown in the matrix has been determined by following the derivation in equation (5): (6) Because the members and appear in the matrix but Ni are functions of ξ, η, it is necessary to use the formulas for the derivation of a composite function: which can be written in matrix form as: The parameter J, which occurs in the equation (8) is called Jacobian. It is used for transformation between coordinate systems, and with its help, it is possible to express specific derivatives.
Jacobian J has the form of matrix two by two and its individual members can be expressed by equations: (9) For calculation the original matrix form of Jacobian J is needed: (10) and also, the inverted matrix: (11) where the Jacobian determinant is calculated using Sarrus rule [29]: (12) The original shape functions (1) have been derived by formulas and the matrix of shape functions has been obtained: (13) Equation (8) can be written in a simplified matrix as follows: (14) The last step of finite element analysis is the derivation of the conductivity matrix, which is analogous to the stiffness matrix in classical problems of elasticity theory, and therefore is denoted by K. For nonlinear diffusion problem is needed to derive also the matrix of capacity C. In publication [23], both matrices were derived for a triangular element and had the following form: (15) ( 16) In the case of a bilinear four-node isoparametric element, these equations could be converted to unit coordinates to create the required formulas for conductivity and capacity matrices: In both cases, it is a definite integral that must be calculated in a specific way. Because of the type, it is appropriate to use an approximation with quadrature formulas. Gauss quadrature was chosen for the project [30].
It should be noted that the newly developed algorithm contains all important formulas and expressions. The formulas related to the calculations of Jacobian J (9), (10), (11), (12), the matrix of the shape functions N (13) and the gradients of these functions B (14) have been applicated. Finally, the equations for conductivity matrix K (17) and capacity C (18) enter the calculation.

Comparison with analytical equation
Comparison of the results of chloride concentrations calculated by FEM model and the analytical model is the subject of this section. The 2D analytical formula presented in the paper [31] is used. The geometric and input parameters have been adopted for presented example (see Table 1). Material parameters of two different concrete mixtures have been introduced into the calculation due to the possible sensitivity to input parameters [32]. Table 1. Input parameters of the introduced models.  Table 2 shows the differences of results computed by FEM and the analytical formula. One could note that the results are very close. The differences are less than 1% which is negligible compared with errors coming from other parameters of the task and is influenced by the geometry of the FEM model.

Discussion and conclusions
The article deals with the improvement of the FEM 2D chloride diffusion model. The prepared algorithm is based on a four-node isoparametric finite element. At first, a detailed derivation of the specific element was presented. Then, the example of calculation by using the simple geometric was shown. Two different concrete mixtures as the input parameters were used for consideration of the potential sensitivity of the model. The results show satisfactory accuracy. However, the presented model is relatively simple with respect to bridge deck considers many simplifications. It should be noted that the analysis assumes that the analytical formula gives accurate results irrespective to the material parameters. This new model can bring improvements for modelling with consideration of the cracks in concrete, various protection of reinforcement in the reinforced concrete structures and other problems of engineering.
This article has been prepared as a part of the research project GACR 18-07949S "Probabilistic Modeling of the Durability of Reinforced Concrete Structures Considering Synergic Effect of Carbonation, Chlorides and Mechanical Action" funded by the Czech Science Foundation. Financial support is highly acknowledged.