Time-Dependent Reliability of Reinforced Concrete Considering Chloride Penetration via Boundary Element Method

STRENGTH DEGRADATION OF STRUCTURAL MATERIALS IS AN INEVITABLE PROCESS, DUE TO DELETERIOUS ACTIONS SUCH AS CORROSION AND FATIGUE. THESE PHENOMENA ARE ALSO TYPICALLY RANDOM, WITH DEGRADATION RATES AND STARTING TIME OF DEGRADATION PROCESS LARGELY UNCERTAIN. IN THE PARTICULAR CASE OF REINFORCED CONCRETE STRUCTURES, CORROSION OF REINFORCING BARS CAUSED BY CHLORIDE IONS IS ONE OF THE MAIN PATHOLOGICAL MANIFESTATIONS. PAST STUDIES ON THE TIME-VARIANT RELIABILITY OF REINFORCED CONCRETE STRUCTURES SUBJECT TO CORROSION HAVE RELIED ON SIMPLIFIED ANALYTICAL MODELS FOR ESTIMATING THE DEPASSIVATION TIME, I.E., THE TIME REQUIRED FOR CHLORIDE CONCENTRATION AT CONCRETE/STEEL INTERFACES TO BECOME CRITICAL. THIS STUDY CONTRIBUTES WITH AN ACCURATE MODELLING OF CHLORIDE DIFFUSION THROUGH CONCRETE USING THE BOUNDARY ELEMENT METHOD, WHICH IS EMPLOYED FOR THE FIRST TIME WITHIN A TIME-VARIANT RELIABILITY FRAMEWORK. CUMULATIVE FAILURE PROBABILITIES ARE EVALUATED IN TIME BY CONSIDERING RANDOM DEPASSIVATION TIMES, RANDOM CORROSION EVOLUTION, AND RANDOM LOAD PROCESSES. THE TIME-VARIANT RELIABILITY PROBLEM IS SOLVED USING MONTE CARLO SIMULATION. AN APPLICATION EXAMPLE IS PRESENTED, DEMONSTRATING THE CAPABILITIES OF THE PROPOSED FRAMEWORK.


INTRODUCTION
The life-time performance of constructed facilities is affected by unavoidable uncertainties in future man-made and natural imposed loads, on the strength and strength degradation of materials, on the quality of construction and on the employed engineering models. A proper way of modelling the impact of uncertainties in loads and degradation mechanisms in the lifetime performance of engineering structures is time-variant reliability analysis.
The performance of constructed facilities can be affected by strength degradation over time. For reinforced concrete structures, in particular, one significant degradation mechanism is corrosion induced by chlorides penetration (Comité Euro-Internacional du Béton, 1992;Gonzalez et al., 1995;Cascudo 1997;Val and Melchers, 1997;Val and Stewart, 2003;Apostolopoulos and Papadakis, 2008;Suo and Stewart, 2009;Zhang et al., 2010). Several studies show that, after depassivation of reinforcement bars due to chloride penetration, structural safety decays sharply (Enright and Frangopol, 1998;Weyers, 1998;Vidal et al., 2007;Dang and François, 2013;Pellizzer et al., 2018).Due to the fast strength reduction, due to corrosion and after depassivation, the largest part of the structures life corresponds to the despassivation time. But this time is also subject to significant uncertainties, due to exposition to random chloride concentrations, variations in concrete cover due to workmanship, and variations in diffusive properties of concrete. Hence, time-variant reliability analysis is a proper tool to investigate the life-time performance of reinforced concrete structures subject to chloride penetration.
The ingress of chloride ions into concrete is a dynamic and nonlinear process. Several studies present analytical formulations to model this phenomenon (Mangat and Molloy, 1994;Vu e Stewart, 2000;Samson et al., 2003;Val et al., 2009;Audenaert et al. 2010;Guzmán et al., 2011). However, analytical solutions apply to very specific domain geometries, and have limitations regarding varying boundary conditions over time. These limitations can be overcome by employing numerical solutions.
In the last years, numerical solutions have been used to address problems like chloride diffusion and reinforcement corrosion. The Finite Element Method (FEM) and the eXtended FEM are two approaches already well explored and reported in the literature (Pan and Lu, 2012;Xiao et al., 2012;Khelifa, 2013;Duddu, 2014;Yoon and Reis, 2017;Al-Alaily et al., 2018). Despite the demonstrated accuracy, these methods require a fine discretization of the time and space domains, which results in very large number of unknowns and, consequently, very large computation time. An efficient alternative, which eliminates spatial domain discretization, is the Boundary Element Method (BEM). Problem equations are written in terms of boundary integrals, with no approximation introduced in the problem domain. Thus, accurate results of diffusive field can be obtained, with one degree less in size of the discretization mesh. In spite of these advantages, only a few studies reported in the literature involve the use of BEM for modeling chloride ion diffusion problems in concrete (Warkus et al., 2006;Piasecka, 2011;Yang et al., 2013;Wang and Chen, 2015). These studies approach the diffusion and strength degradation problems in a deterministic manner.
In recent years, some studies have been done regarding time-variant reliability of reinforced concrete structures subjected to chloride ion penetration and corrosion. Some studies addressed use of inspection and maintenance to maintain reliability indexes over a target value, over the structural lifetime (Mori and Ellingwood, 1994;Enright and Frangopol, 1999;Biondini et al., 2006). Structural reliability studies addressing the interaction between climate variability and the behavior of structures over time should also be mentioned (El Hassan et al., 2010;Stewart et al., 2012;Bastidas-Arteaga and Stewart, 2015). Other studies dealing more specifically with the consequences of uncertainties on the mechanical behavior of structural elements have also been conducted (Thoft-Christensen, 1998;Stewart, 2004;Duprat, 2007;Xiang and Zhao, 2007;Ghosh and Padgett, 2010;Simon et al., 2010).
The studies cited in the last paragraph presented various contributions, but presented some common limitations, like use of approximate reliability methods such as FORM and SORM, considering loads as time-invariant, or using simple analytical diffusion equations. FORM and SORM are known to be inaccurate for problems involving highly nonlinear or multiple limit state equations. Acting live loads are better described as random processes of time, especially in case of bridges and viaducts. In addition, simplified analytical models for chloride diffusion modeling were used in most of these studies. Use of one-dimensional diffusive chloride models in two-dimensional cross-sections of structural elements leads to large differences in reliability, as shown by Sørensen (1996), Val and Trapper (2008) and Bastidas-Arteaga et al. (2011). Finally, boundary conditions like flow or concentration of chloride ions at element surface vary in time, but are assumed time independent in most cited studies.
Addressing the above challenges, the present study proposes a framework for the robust and accurate analysis of time-variant reliability of reinforced concrete structures subjected to corrosion caused by penetration of chlorides. A Monte Carlo based time-dependent reliability analysis is performed, in order to explicitly represent evolution of the process. Chloride ion diffusion through concrete is modelled using the Boundary Element Method (BEM); metal loss due to corrosion after depassivation of reinforcement bars is considered, and loading is modelled as a random process of time. At any time, structural failure can occur due to acting loads becoming larger than instantaneous strength. To the best of the author´s knowledge, this is the first time that BEM is employed in the context of time variant reliability analysis involving concrete strength degradation due to chloride penetration.

Governing Equations and Integral Representation
The Poisson equation is appropriate for describing several potential problems, like diffusion, torsion, thermal conductivity and conduction of fluids. The approach accounts for time-independent and time-dependent boundary conditions. For the stationary case (time independent), the Poisson's equation is as follows: in which represents the potential and indicates the domain term. In the particular case of b=0, the Poisson's equation leads to the Laplace solution.
The solutions of Eq. (1) require enforcement of the following boundary conditions: where � describes the prescribed potential values, � describes the prescribed flux value, Γ 1 and Γ 2 represent the boundaries in which potential or flux are prescribed, respectively. The complete body boundary is given by: In addition, flux is associated to the potential as follows: = .Thus, flux represents the directional derivative of with respect to the outward normal vector . Figure (1) represents the domain, boundaries and normal vector for an arbitrary body. As previously mentioned, Eq. (1) represents stationary problems. In order to handle transient fields, this equation must be modified by introducing a time dependent term. Therefore, the transient potential problem is governed by the following differential equation: in which indicates the domain-related parameter, such as thermal diffusivity or diffusion coefficient, and is time. Similarly to Eq. (1), solutions of Eq.
(2) are obtained by enforcing boundary conditions. Thus, potential and flux values at the boundary are prescribed at each time step. The diffusion problems addressed in this study employ Eq. (2). Therefore, the potential values indicate chloride concentration along time, whereas flux represents chloride flux along time.
The differential equation (Eq. 2) can be transformed into a boundary integral representation by applying either the finite differences technique, the Laplace transform or the fundamental time dependent solutions (Carslaw and Jaeger, 1959;Curran et al., 1980;Wrobel and Brebbia, 1981;Wrobel and Brebbia, 1987;Wrobel, 2002). In this study, the latter approach was employed. The boundary integral representation for the transient potential problem is obtained by applying initially the Green's second identity in Eq. (2). Then, integration by parts is employed. Finally, the classical limit process is performed, in which the Somigliana's equation is evaluated at the body's boundary. The above-mentioned procedures lead to the following integral representation, in which temporal and spatial integrations are required: in which indicates the source points, refers the field points, 0 represents the initial time, represents the observation time, * and * indicate the fundamental time-dependent solutions for potential and flux, respectively, and is the classical BEM free term. The free-term is the unit for source points positioned at the domain. This parameter is equal to 0.5 for source points positioned at smooth boundary geometries. For the plane case, the time-dependent fundamental solutions are as follows (Wrobel, 2002;Carslaw and Jaeger, 1959) where = − , represents the distance between the source and the field points, and = , .

Algebraic Representation
Equation (3) provides the values of potential and flux at the boundary for a given time step. The solutions for this equation require spatial and temporal integration schemes. The spatial integration is handled by discretizing the entire boundary geometry into boundary elements over which polynomial shape functions approximate geometry and diffusion fields. The temporal integration is performed through the constant approximation. Such approach enables the analytical integration of the time dependent kernels required by the BEM.
Consequently, discretization of the body's boundaries into boundary elements and discretization of the time interval − 0 into time steps enable rewriting Eq. (3) into its discrete form as follows: represents the potential value at the time for a given source point . Equation (6) can be rewritten bearing in mind that fundamental solutions * and * are constantly integrated along time. Moreover, one locates source points solely at smooth boundaries in order to ensure enforcement of accurate boundary conditions. Consequently: Because the temporal integration is handled analytically, the last equation is further rewritten as follows (Wrobel, 2002): The analytical integration scheme along time transforms * into * and * into * . These new integral kernels are as follows (Wrobel, 2002): Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e316 5/17 in which 1 indicates the exponential-integral function.
As usual in BEM, the integrals illustrated in Eq. (8) lead to the classical influences matrices and . Then: in which and represent the influence matrices, which are defined as follows: Therefore, the classical algebraic BEM system of equations is achieved as follows: where and are matrices, which contain the influence coefficients obtained from Eq. (13) and Eq. (14). represents the amount of collocation points into the boundary mesh.
The solution of Eq. (15) requires enforcement of boundary conditions. For this purpose, the classical columns change procedure between and matrices must be performed. As usual in BEM, the unknown values at the boundary are moved to the left side of Eq. (15) Eq. (16) accounts for the fundamental kernels illustrated in Eq. (9) and Eq. (10). Consequently, Eq. (16) is time dependent, i.e., the values at the present time depend on the history of fields values. As a result, a time marching process is required for its proper solution. The constant approach is applied in the present study (Wrobel, 2002). Therefore, the time marching process is as follows: Where is assessed as follows: Consequently, for the proper solution of Eq. (16), the vector must be added to the vector in order to account for the time history effects.
The above-mentioned integral equations have been solved numerically. However, the singular nature of the fundamental kernels requires special attention. The standard Gauss-Legendre quadrature has been applied when the integrated boundary element does not contain the source point. For this case, the kernels are regular and ≠ 0. The kernels become improper for singular elements. In this case, the integrated boundary element contains the source point and at some point of its length, = 0. The singularity-subtraction method (SSM) handles accurately and properly the singular integrals. The SSM handles the improper integrals through a semi-analytical scheme. Then, the application of SSM leads to a regular kernel and to a new singular kernel. Nevertheless, the latter has analytical solution. The regular part is assessed through the standard Gauss-Legendre quadrature, whereas the analytical solution is obtained by the Cauchy principal value. It is worth mentioning that SSM enables the diffusion modelling accounting for high-order boundary elements. In addition, geometrically curved elements may be utilized, which provides better approximation for non-straight geometries. The authors demonstrate the analytical solutions for SSM in transient diffusion problems taking into account the Cauchy principal value in Appendix A.
Finally, Eq. (6) also enables the assessment of internal potential fields. In this case, the free term equals 1 once the source point is positioned into the domain. In addition, for this case, Eq. (6) is evaluated exclusively through the standard Gauss-Legendre quadrature, because the kernels are regular.

CORROSION PROCESS
One of the first qualitative idealizations of the useful life of reinforced concrete structures subjected to chloride ions penetration was proposed by Tuutti (1982). In this model, structural useful life is divided in two phases: initiation and propagation. In the initiation phase, the chloride ions present on the outer surface of the structure move through the micro pores of the concrete. When the concentration of chlorides at the concrete/reinforcement interface reaches a threshold value, reinforcement depassivation occurs, ending the initiation phase. Propagation phase then begins, where the corrosion chemical reactions occur, resulting in deleterious mechanical effects on steel and concrete. Initiation and propagation phases occur in parallel within the same reinforced concrete element, i.e., some reinforcement bars may have already undergone depassivation, while others have not. This is due to the high randomness present in the variables involved in the corrosive process.
In the initiation phase, the main transport mechanism of chloride ions through the concrete's micro pores is diffusion. The diffusion of chloride ions in concrete is often modeled using Fick's laws (Vu and Stewart, 2000;Samson et al., 2003;Val et al., 2009;Guzmán et al., 2011). Two limitations of most differential equation solutions presented in the literature, describing Fick's law, are that they are obtained from semi-infinite domains, and consider constant boundary conditions over time. Since it is extremely important to accurately predict when the corrosion propagation period will begin, in the present work an alternative approach is presented, based on the Boundary Element Method (BEM). Among the advantages of BEM model is the possibility of varying boundary conditions over time, and the possibility of analyzing two-dimensional domains with any contours of geometry. Therefore, the formulation presented here seeks to be a guide for a more reliable modeling of the reality of the chloride diffusion phenomenon.
In the propagation phase, one of the main effects of corrosion is the reduction of reinforcement steel area. Chemical reactions that consume the constituent metal of the rebars produce two basic types of corrosion: uniform and pitting. In this work, uniform corrosion of the reinforcements considered. The model adopted is presented in the work of Thoft-Christensen and Hansen (1994) and Val and Melchers (1997). Based on Faraday's laws, the authors propose an equation to predict the residual diameter of the corroded rebar as follows: where 0 is the original diameter of the non-corroded rebar in mm, is the corrosion rate inμA/cm 2 (microampere per square centimeter) and is the elapsed time after rebar depassivation, in years. Once the diameter of the corroded bar ( ) is known, it is possible to calculate the steel area of each of the steel rebars and hence the total steel area of a given reinforced concrete cross section. The corrosion rate in this work is calculatedusing the empirical equation presented by Vu and Stewart (2000).This equation is based on measurements available in the literature and on the conversion of the oxygen diffusion rate to the corrosion rate by considering the percentage of corrosion products and the molecular equations of corrosion in the cathodic zone. The expression that relates rate of corrosion with elapsed time after the depassivation of the reinforcements is: where ( ) is the corrosion rate in / 2 (microampere per square centimeter), / is the water/cement ratio, is the reinforcement concrete cover thickness in mm and is the elapsed time after the depassivation of the rebar in years. In this study, other deleterious mechanisms of the propagation phase, such as reinforcement stress reduction, modification in the stress transfer mechanism between steel and concrete, and the appearance of internal cracks are not taken into account.

TIME-VARIANT RELIABILITY
Structural safety assessment approaches found in most normative codes are based in semi-probabilistic methods, and thus adopted for structural design (Brazilian Association of Technical Standards, 2014). Despite that, structural safety can be explicitly considered and quantified as a probability that the structure won't behave as intended, the probability of failure. Instant probabilities of failure are usually employed in this context. The effect of input random variables is studied as uncertainties are propagated trough structural systems, resulting in stochastic outputs. When stochastic timevariant loads are present, instant probabilities of failure are no longer representative, since structural failure conditions may vary over time. The same is true for resistance degradation scenarios, where accumulated damage or changes in the materials can also affect structural safety, usually diminishing it in time. In this context, a more general approach must be considered, so that structural safety can be addressed as a function of time. In order to do that, consider a set of = + elements ( , ), where denotes the outcome in the space of outcomes Ω. The random variables of the problem are a subset of ( , ), and ordered as ( ), = 1, … , . These variables are employed to represent most resistance parameters, such as structural material and geometrical properties. Also, let all ( , ), = + 1, … , + be another subset of ( , ), gathering stochastic processes, so that time-variant loads are also represented.
A limit state function is defined, representing the behavior of a structural component or system. It denotes safety when assumes values greater than zero and failure otherwise. Thus, the function ( , , ( , ))is used to define the failure domain, for undesirable structural responses, and safety domain, for adequate structural behavior: In Eq. (21), vector groups some deterministic parameters of the problem, which can be employed, for instance, in structural design optimization. The instantaneous probability of failure , for a certain limit state, at a certain time = is then given by: where ℙ denotes the probability operator, and ( ) is the joint probability density function of the random variables for a certain limit state and a configuration at a time . In time-dependent problems, it is more convenient to work with the concept of cumulative probability of failure ( 1 , 2 ). For a certain design configuration and for a given limit state, the probability of a structure not behaving as intended within the time interval [ 1 , 2 ] is given by: Approximate methods have been proposed for the computation of (Andrieu- Renaud et al., 2004, Sudret, 2008. Although efficient, they are not very accurate, particularly when highly non-linear limit states or multiple design points are present. In order to overcome these limitations, more accurate simulation based approaches can be employed, at the cost of considerable computational burden. In this paper, a simulation based approach is adopted, and the corresponding increase in computational cost is diminished by the employment of BEM and normative code-based analytical limit states.

Monte Carlo-based estimation of the cumulative failure probability
When coupled with efficient techniques (Kroetz et al., 2017b), direct simulation becomes feasible in the context of structural optimization (Kroetz et al., 2017a, Kroetz et al., 2020, Moustapha et al. 2016. In order to simulate structural performance from time-dependent limit state equations, stochastic processes are discretized as finite sets of random variables (Sudret and Der Kiureghian, 2000). Hence, the limit state equations that describe the problem are evaluated considering the fixed values sampled from random variables, and the values assumed by the stochastic processes in each instant of the discretized time-series. To represent the stochastic processes, the technique known as "expansion optimal linear estimation" (EOLE) is employed herein, as proposed by Li and Der Kiureghian (1993). In this work, the random process discretization is performed using the software UQlab (Marelli and Sudret, 2014). Let ( , ) be a Gaussian random process, with mean , standard deviation and auto-correlation coefficient function ( 1 , 2 ). Also, consider that time points are selected in the time-series interval [0: ], with 1 = 0 and = . The EOLE representation of the stochastic process reads: In this notation,{ ( ), = 1, … , } are independent standard normal variables and { , = 1, … , } are the first eigenvectors and eigenvalues of the correlation matrix , sorted in decreasing order, where = � , �, , ={1,..., }. Since the expression represents a numerical simulation model, the expansion must be truncated to ≤ terms, where is the value that defines the order of expansion. Once all random variables and random processes are sampled, for each studied structural configuration , the limit state equation � , , ( , )� is evaluated in all points of the discretized time interval [0: ].
In order to estimate the cumulative failure probability, the simulated values are stored in an array of dimension 1 × , where is the number of time instants in which the limit state equation is discretized. Each position of this array corresponds to a time = ( − 1)Δ where Δ = −1 is the sampling step, so that an uniform discretization is assumed. For each , a corresponding counter that is increased whenever presents the first outcrossing in the interval [ ; +1 ]. The Monte Carlo estimation for the cumulative probability of failure, considering samples is then given by:

PROPOSED METHODOLOGY
In order to address the safety of a concrete structure throughout its whole life span, the first step of the proposed methodology is to quantify the time for depassivation of reinforcement bars, thus triggering corrosion processes. The analysis is divided in two independent steps, henceforth addressed as Problem 1 and Problem 2. In Problem 1, uncertainties regarding the chloride penetration process, as detailed in Section 3, are considered as input to the BEM based analysis, performed as explained in Section 2. The input variables to this analysis are the chloride threshold value , the surface chloride concentration , the diffusion coefficient in concrete and the thickness of the concrete cover . This problem is solved 10 6 times, in order to provide a million depassivation times, which are considered as inputs to the second problem. Problem 2 consists in the Monte Carlo based time-dependent reliability analysis. This step also considers uncertainties in the structural parameters: random variables concrete strength and steel strength and the loads of the problem, hereby generically represented by random process ( ). The statistical parameters of the random variables and processes depend on the problem, and are addressed in Section 6. The relationship between Problem 1 and Problem 2 is illustrated in Figure (2). Thus, a time-dependent reliability framework is proposed for degrading reinforced concrete beams. The cumulative probabilities of failure can be accurately computed since a Monte Carlo direct simulation method is employed, as described in Section 4.1. The computational burden is reduced because "Problem 1" is efficiently solved using BEM, hence accurate low-dimension mesh models are employed. The structural behavior considered in the analysis reflects a practice of engineering structural design; hence analytical models are adopted, as suggested by the Brazilian normative code for reinforced concrete beams (Brazilian Association of Technical Standards, 2014).

APPLICATION EXAMPLE
Consider the 4 meter span beam represented in Figure (3). The beam has rectangular cross-section of 40cmX20cm, and 3 longitudinal reinforcement bars with diameter 12.5mm. It is also subject to a time-variant load ( ), described by a random process with Gaussian autocorrelation function, with correlation length of 1 year.  Table 1 shows the random quantities involved in the chloride diffusion problem, and Table 2, the random quantities involved in the structural reliability analysis. The boundary was discretized in 8 quadratic discontinuous and isoparametric boundary elements. The number of Gauss points per element is 8. Figure (4) illustrates adopted boundary element mesh and boundary conditions. The simulation of chloride ion penetration was performed until 200 years, discretized in 50 time steps. Depassivation occurrence time of each rebar was determined by means of interpolation between the successive time steps within which depassivation occurred. A mesh convergence analysis was previously done in order to determine a suitable boundary discretization for the problem.   Reliability analysis was carried out over a 50 year time-span, considering 10 6 samples in each of the 500 discrete timepoints, so that each correlation length is discretized in 10 time steps. The limit state considered compares resisting bending moment ( ) with load effect bending moment ( ). The resisting bending moment is given by normative recommendations of the Brazilian standard for reinforced concrete design (Brazilian Association of Technical Standards, 2014): In Eq. (26) = 0.8 when < 50MPa and = 0.8 − −50 400 when 50MPa ≤ ≤ 90MPa, is the height of the neutral axis in the cross section, and is the remaining area of steel. The steel cross-section area is only reduced by corrosion after depassivation by chlorides. Time-independent corrosion is considered in the first year after depassivation (Eq. 19), and time-dependent corrosion is adopted from this point on (Eq. 20). The structure is also considered to fail when the normative ductility criterion is violated, once again considering NBR6118 specification, which is indicated by: As mentioned, the time-variant reliability problem is solved by generating 10 6 samples of the depassivation time, and same number of samples of resisting ( ( )) and load effect ( ( )) bending moments. One sample realization of the load effect process, and one sample of the strength degradation due to depassivation and corrosion are shown in Figure

CONCLUDING REMARKS
This paper has presented a framework for the robust and accurate analysis of time-variant reliability of reinforced concrete structures subjected to corrosion caused by penetration of chlorides. A Monte Carlo based time-dependent reliability analysis was performed, in order to explicitly represent evolution of the loading and corrosion processes. Chloride ion diffusion through concrete was modelled using the Boundary Element Method (BEM); metal loss due to corrosion after depassivation of reinforcement bars was considered, and loading was modelled as a random process of time. One distinguished feature of the framework proposed herein is the use of BEM in order to accurately model depassivation of reinforcement bars, due to penetration of chlorides. To the best of the author´s knowledge, this is the first time such an accurate BEM model is employed in the context of time-variant reliability analysis.

ACKNOWLEDGMENTS
Valuable comments by the anonymous reviewers are cheerfully acknowledged. The authors acknowledge funding of this research project by CAPES (Brazilian Higher Education Council -Finance Code 001). The fourth author acknowledges funding by CNPq (grant n. 306373/2016-5) and FAPESP (grant n. 2017/01243-5).
The singularity is observed in U k * when k = 1, i.e., in the first time step. In addition, the singularity is herein observed when the source point approaches to the field point, i.e., when r → 0. Thus, when k = 1, Eq. (A.1) becomes: The variables a f k and a 0 k are positive because κ > 0, r ≥ 0 and t F ≥ t f k > t 0 k . In Eq. (A.2), when r is nil and t f 1 → t F , a f 1 → +∞. Thus, the exponential-integral function is evaluated at + ∞, which results into a nil value. In equation (A.3) when r is nil, a 0 1 → 0, because the denominator (t F − t 0 1 ) is always greater than zero (t F ≥ t f 1 > t 0 1 ). Thus, the exponentialintegral function has nil value, which results into the following singularity: The singularity arising in this equation is logarithmic for distinct intervals of the exponential-integral function argument, 0 ≤ z ≤ 1and 1 < < ∞. Thus, the kernel containing u * can be regularized by the following equation: whereξ 0 is the dimensionless coordinate of the source point, ξ the dimensionless coordinate of the field point, x i the real coordinate of a given point, ϕ the shape function, ϕ ,i the shape function derivative and the Jacobian is given by: J = � �ϕ ,i (ξ 0 )x i � 2 , the real distance r between the source and the field points is evaluated as follows: where ε = |ξ − ξ 0 | is the dimensionless distance between the source and field points. Substituting dΓ = Jdξ in Eq. (A.6) one has: ∫ U 1 * �ξ, x� dΓ j Γ j = 1 4πκ ∫ E 1 (a 0 1 )ϕ(ξ)J(ξ)dξ