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

ABSTRACT

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 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. 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.

Keywords:
Reinforced Concrete; corrosion, time-dependent reliability; Boundary Element Method

Graphical Abstract

1. 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.

2. THE TRANSIENT BEM FORMULATION

2.1 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:

${\nabla }^{2}u+b=0$ (1)

in which $u$ represents the potential and $b$ 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:

• Dirichlet: $u=\stackrel{-}{u}$ in ${\mathrm{\Gamma }}_1$

• Neumann: $q=\stackrel{-}{q}=\frac{du}{d\eta }$ in ${\mathrm{\Gamma }}_2$

where $\stackrel{-}{u}$ describes the prescribed potential values, $\stackrel{-}{q}$ describes the prescribed flux value, ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ represent the boundaries in which potential or flux are prescribed, respectively. The complete body boundary is given by: ${\mathrm{\Gamma }=\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2$. In addition, flux is associated to the potential as follows: $q=\frac{du}{d\eta }$ .Thus, flux represents the directional derivative of $u$with respect to the outward normal vector $\eta$. Figure (1) represents the domain, boundaries and normal vector for an arbitrary body.

Figure 1
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:

${\nabla }^{2}u-\frac{1}{\kappa }\frac{\partial u}{\partial t}=0$ (2)

in which $\kappa$ indicates the domain-related parameter, such as thermal diffusivity or diffusion coefficient, and $t$ 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:

$c\left(\underset{\_}{\xi }\right)u\left(\underset{\_}{\xi },t_F\right)=\kappa \underset{t_0}{\overset{t_F}{{\displaystyle \int }}}\underset{\Gamma }{\overset{}{{\displaystyle \int }}}q\left(\underset{\_}{x},t\right)u^{*}\left(\underset{\_}{\xi },\underset{\_}{x},t_F,t\right)d\Gamma \left(\underset{\_}{x}\right)dt-\kappa \underset{t_0}{\overset{t_F}{{\displaystyle \int }}}\underset{\Gamma }{\overset{}{{\displaystyle \int }}}u\left(\underset{\_}{x},t\right)q^{*}\left(\underset{\_}{\xi },\underset{\_}{x},t_F,t\right)d\Gamma \left(\underset{\_}{x}\right)dt$ (3)

in which $\underset{\_}{\xi }$ indicates the source points, $\underset{\_}{x}$ refers the field points, $t_0$ represents the initial time, $t_F$ represents the observation time, $u^{\mathrm{*}}$ and $q^{\mathrm{*}}$ indicate the fundamental time-dependent solutions for potential and flux, respectively, and $c$ 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):

$u^{*}\left(\underset{\_}{\xi },\underset{\_}{x},t_F,t\right)=\frac{1}{4\pi \kappa \tau }\exp \left(-\frac{r^2}{4\kappa \tau }\right)$ (4)

$q^{*}\left(\underset{\_}{\xi },\underset{\_}{x},t_F,t\right)=\frac{r\frac{\partial r}{\partial \eta }}{8\pi {\kappa }^2{\tau }^2}\exp \left(-\frac{r^2}{4\kappa \tau }\right)$ (5)

where $\tau =t_F-t$, represents the distance between the source $\underset{\_}{\xi }$ and the field $\underset{\_}{x}$ points, and $\frac{\partial r}{\partial \eta }={r,}_k{\eta }_k$.

2.2 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 $\Gamma$ into $Ne$ boundary elements and discretization of the time interval $t_F-t_0$ into $Nt$ time steps enable rewriting Eq. (3) into its discrete form as follows:

$c_i{u}_i^{Nt}=\kappa \sum _{j=1}^{Ne}\sum _{k=1}^{Nt}{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}{\int }_{t_0^k}^{t_f^k}{q}^k{u}^{\mathrm{*}}dtd{\mathrm{\Gamma }}_{\mathrm{j}}-\kappa \sum _{j=1}^{Ne}\sum _{k=1}^{Nt}{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}{\int }_{t_0^k}^{t_f^k}{u}^k{q}^{\mathrm{*}}dtd{\mathrm{\Gamma }}_{\mathrm{j}}$ (6)

where$u_i^{Nt}$ represents the potential value at the time $t^{Nt}$ for a given source point $i$. Equation (6) can be rewritten bearing in mind that fundamental solutions $u^{*}$ and $q^{*}$ are constantly integrated along time. Moreover, one locates source points solely at smooth boundaries in order to ensure enforcement of accurate boundary conditions. Consequently:

$0.5u_i^{Nt}=\kappa \sum _{j=1}^{Ne}{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}\sum _{k=1}^{Nt}{q}^k{\int }_{t_0^k}^{t_f^k}{u}^{\mathrm{*}}dtd{\mathrm{\Gamma }}_{\mathrm{j}}-\kappa \sum _{j=1}^{Ne}{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}\sum _{k=1}^{Nt}{u}^k{\int }_{t_0^k}^{t_f^k}{q}^{\mathrm{*}}dtd{\mathrm{\Gamma }}_{\mathrm{j}}$ (7)

Because the temporal integration is handled analytically, the last equation is further rewritten as follows (Wrobel, 2002):

$0.5u_i^{Nt}=\kappa \sum _{k=1}^{Nt}\left(\sum _{j=1}^{Ne}{q}_j^k{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}{U}_k^{\mathrm{*}}d{\mathrm{\Gamma }}_{\mathrm{j}}-\sum _{j=1}^{Ne}{u}_j^k{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}{Q}_k^{\mathrm{*}}d{\mathrm{\Gamma }}_{\mathrm{j}}\right)$ (8)

The analytical integration scheme along time transforms $u^{*}$ into $U^{*}$and $q^{*}$into $Q^{*}$. These new integral kernels are as follows (Wrobel, 2002):

$U_k^{\mathrm{*}}={\int }_{t_0^k}^{t_f^k}{u}^{\mathrm{*}}dt=\frac{1}{4\pi \kappa }\left\{E_1\left[\frac{r^2}{4\kappa (t_F-t_0^k)}\right]-E_1\left[\frac{r^2}{4\kappa (t_F-t_f^k)}\right]\right\}$ (9)

$Q_k^{\mathrm{*}}={\int }_{t_0^k}^{t_f^k}{q}^{\mathrm{*}}dt=\frac{1}{2\pi r\kappa }\frac{\partial r}{\partial \eta }\left\{\exp \left[-\frac{r^2}{4\kappa (t_F-t_0^k)}\right]-\exp \left[-\frac{r^2}{4\kappa (t_F-t_f^k)}\right]\right\}$ (10)

in which $E_1$ indicates the exponential-integral function.

As usual in BEM, the integrals illustrated in Eq. (8) lead to the classical influences matrices $H$ and $G$. Then:

$0.5u_i^{Nt}=\sum _{k=1}^{Nt}\left(\sum _{j=1}^{Ne}{G_{ij}^{k}q}_j^k-\sum _{j=1}^{Ne}{\widehat{H_{ij}^k}u}_j^k\right)$ (11)

in which $H$ and $G$ represent the influence matrices, which are defined as follows:

$G_{ij}^k={\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}^{}{U}_k^{\mathrm{*}}d{\mathrm{\Gamma }}_{\mathrm{j}}$ (12)

$\widehat{H_{ij}^k}={\int }_{{\Gamma }_j}^{}{Q}_k^{\mathrm{*}}d{\Gamma }_j$ (13)

$H_{ij}^k=\left\{\begin{array}{c}\widehat{H_{ij}^k}+0.5ifk=1andi=j\\ \widehat{H_{ij}^k}otherwise\end{array}\right.$ (14)

Therefore, the classical algebraic BEM system of equations is achieved as follows:

$HU=GQ$ (15)

where $H$ and $G$ are $NxN$ matrices, which contain the influence coefficients obtained from Eq. (13) and Eq. (14). $N$ 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 $H$ and $G$ matrices must be performed. As usual in BEM, the unknown values at the boundary are moved to the left side of Eq. (15), whereas the prescribed values are positioned into the right side of this equation. This procedure enable rewriting Eq. (15) as follows:

$AX=F$ (16)

in which $A$ represents the matrix composed of coefficients from $H$ and $G$ associated to unknown values, $X$ indicates the vector containing the unknown values at the boundary and $F$ is the vector obtained from multiplication of the known values and their respective influence coefficients. The solution of Eq. (16) provides the unknown values at the body’s boundary.

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:

$H_{t_f}{U}_{t_f}=G_{t_f}{Q}_{t_f}+S_f$ (17)

Where $S_f$ is assessed as follows:

$S_f=-\sum _{j=1}^{f-1}{H}_j{U}_j+\sum _{j=1}^{f-1}{G}_j{Q}_j$ (18)

Consequently, for the proper solution of Eq. (16), the $S_f$ vector must be added to the $F$ 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 $r\ne 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, $r=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 $c$ 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.

3. 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:

$d\left(t_c\right)=d_0-0.0232i_{corr}{t}_c$ (19)

where $d_0$ is the original diameter of the non-corroded rebar in mm, $i_{corr}$ is the corrosion rate inμA/cm² (microampere per square centimeter) and $t_c$is the elapsed time after rebar depassivation, in years. Once the diameter of the corroded bar $d\left(t_c\right)$ 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 $i_{corr}$ 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:

$i_{corr}\left(t_c\right)=\left[\frac{37.8{\left(1-w/c\right)}^{-1.64}}{cov}\right]0.85t_c^{-0.29}$ (20)

where $i_{corr}\left(t_c\right)$ is the corrosion rate in $\mu A/c{m}^2$ (microampere per square centimeter), $w/c$ is the water/cement ratio, $cov$ is the reinforcement concrete cover thickness in mm and $t_c$ 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.

4. 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 time-variant 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 $M=p+q$ elements $\mathit{X}\left(t,\omega \right)$, where $\omega$ denotes the outcome in the space of outcomes $\mathrm{\Omega }$. The random variables of the problem are a subset of $\mathit{X}\left(t,\omega \right)$, and ordered as $X_j\left(\omega \right),j=1,\dots ,p$. These variables are employed to represent most resistance parameters, such as structural material and geometrical properties. Also, let all $X_k\left(\omega ,t\right),k=p+1,\dots ,p+q$ be another subset of $\mathit{X}\left(t,\omega \right)$, 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 $g(\mathit{d},\mathit{t},\mathit{X}(t,\omega \left)\right)$is used to define the failure domain, for undesirable structural responses, and safety domain, for adequate structural behavior:

$\begin{array}{l}{D}_f\left(t,d\right)=\left\{d,X\left(t,\omega \right):g\left(d,t,X\left(t,\omega \right)\right)\le 0\right\}isthefailuredomain\\ D_s\left(t,d\right)=\left\{d,X\left(t,\omega \right):g\left(d,t,X\left(t,\omega \right)\right)>0\right\}isthesafedomain\end{array}$ (21)

In Eq. (21), vector $\mathit{d}$ groups some deterministic parameters of the problem, which can be employed, for instance, in structural design optimization. The instantaneous probability of failure $P_{fi}$, for a certain limit state, at a certain time $t=\tau$ is then given by:

$P_{fi}\left(\mathit{d},\tau \right)=\mathbb{P}\left(g(\mathit{d},\tau ,\mathit{X}(\tau ,\omega \left)\right)\le 0\right)={\int }_{D_f(\mathit{d},\tau )}^{}f\left(\mathit{x}\right)d\mathit{x}$ (22)

where $\mathbb{P}$ denotes the probability operator, and $f\left(\mathit{x}\right)$ is the joint probability density function of the random variables for a certain limit state and a configuration $\mathit{d}$ at a time $\tau$. In time-dependent problems, it is more convenient to work with the concept of cumulative probability of failure $P_{fc}(t_1,t_2)$. For a certain design configuration $\mathit{d}$ and for a given limit state, the probability of a structure not behaving as intended within the time interval $\left[t_1,t_2\right]$ is given by:

$P_{fc}\left(\mathit{d},t_1,t_2\right)=\mathbb{P}\left[\exists \tau \in \left[t_1,t_2\right]:g\left(\mathit{d},\mathit{t},\mathit{X}\left(t,\omega \right)\right)\le 0\right]$ (23)

Approximate methods have been proposed for the computation of $P_{fc}$ (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.

4.1 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 $\mathit{X}\left(t,\omega \right)$ be a Gaussian random process, with mean $\mu$, standard deviation $\sigma$ and auto-correlation coefficient function ${\rho }_X(t_1,t_2)$. Also, consider that $\mathcal{P}$ time points are selected in the time-series interval $\left[0:T\right]$, with $t_1=0$ and $t_{\mathcal{P}}=T$. The EOLE representation of the stochastic process reads:

$\mathit{X}\left(t,\omega \right)=\mu +\sigma \sum _{i=1}^N\frac{{\xi }_i\left(\omega \right)}{\sqrt{{\lambda }_i}}{\varphi }_i^T{\mathit{C}}_{t,t_i}\left(t\right)$ (24)

In this notation,$\{{\xi }_i\left(\omega \right),i=1,\dots ,\mathcal{P}\}$ are independent standard normal variables and $\{{\varphi }_i^{}{\lambda }_i,i=1,\dots ,N\}$ are the first eigenvectors and eigenvalues of the correlation matrix $\mathit{C}$, sorted in decreasing order, where $C_{ij}={\rho }_X\left(t_i,t_j\right),i,j=${1,...,$\mathcal{P}$}. Since the expression represents a numerical simulation model, the expansion must be truncated to $N\le \mathcal{P}$ terms, where $\mathcal{P}$ is the value that defines the order of expansion. Once all random variables and random processes are sampled, for each studied structural configuration $\mathit{d}$, the limit state equation $g\left(\mathit{d},\mathit{t},\mathit{X}\left(t,\omega \right)\right)$ is evaluated in all points of the discretized time interval $\left[0:T\right]$.

In order to estimate the cumulative failure probability, the simulated values are stored in an array $G$ of dimension $1\times N$, where $N$ is the number of time instants in which the limit state equation is discretized. Each position $i$ of this array corresponds to a time $t_i=(i-1)\mathrm{\Delta }t$ where $\mathrm{\Delta }t=\frac{T}{N-1}$ is the sampling step, so that an uniform discretization is assumed. For each $t_i$, a corresponding counter $k_j$ that is increased whenever $g$ presents the first outcrossing in the interval $[t_i;t_{i+1}]$. The Monte Carlo estimation for the cumulative probability of failure, considering $N_{MC}$ samples is then given by:

$P_{fc}(0,t_i)=\frac{1}{N_{MC}}{k}_j$ (25)

5. 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 $C_{limit}$, the surface chloride concentration $C_s$, the diffusion coefficient in concrete $\kappa$ and the thickness of the concrete cover $CT$. 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 $f_c$ and steel strength $f_y$ and the loads of the problem, hereby generically represented by random process $q\left(t\right)$. 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).

Figure 2
Proposed methodology.

6. 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 $q\left(t\right)$, described by a random process with Gaussian autocorrelation function, with correlation length of 1 year.

Figure 3
Reinforced concrete beam of the case study.

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.

Figure 4
Boundary Element mesh and boundary conditions.

Table 1
Statistical parameters for the diffusion analysis with BEM.

Parameter	Distribution type	Mean	COV
Chloride threshold value $C_{limit}$	Uniform	0.90 kg/m3	0.19
Surface chloride concentration $C_s$	Lognormal	1.15 kg/m3	0.50
Diffusion coefficient in concrete $\kappa$	Lognormal	0. 6742277 cm2/year	0.75
Cover Thickness $\mathrm{C}\mathrm{T}$	Normal	3.00 cm	0.50

Table 2
Statistical parameters for the time-variant reliability analysis.

Parameter	Distribution type	Mean	COV
Concrete Strength $f_c$	Lognormal	30 MPa	0.20
Steel Strength $f_y$	Lognormal	500 MPa	0.20
Load $q\left(t\right)$	Gaussian random process	10 kN/m	0.20

Reliability analysis was carried out over a 50 year time-span, considering ${10}^6$ samples in each of the 500 discrete time-points, so that each correlation length is discretized in 10 time steps. The limit state considered compares resisting bending moment ($M_R$) with load effect bending moment ($M_S$). The resisting bending moment is given by normative recommendations of the Brazilian standard for reinforced concrete design (Brazilian Association of Technical Standards, 2014):

$g\left(\mathit{x},t\right)=M_R\left(\mathit{x},t\right)-M_S\left(\mathit{x},t\right)$ $M_R\left(\mathit{x},t\right)=A_s\left(t\right)f_y\left(d-\frac{\lambda x}{2}\right)$ $M_S\left(\mathit{x},t\right)=\frac{q\left(t\right)l^2}{8}$ (26)

In Eq. (26)$\lambda =0.8$ when $f_c<50\mathrm{M}\mathrm{P}\mathrm{a}$ and $\lambda =0.8-\frac{f_c-50}{400}$when $50\mathrm{M}\mathrm{P}\mathrm{a}\le f_c\le 90\mathrm{M}\mathrm{P}\mathrm{a}$, $x$ is the height of the neutral axis in the cross section, and $A_s$ 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:

$\left\{\begin{array}{c}\frac{x}{d}\le 0.45if{f}_c\le 50\mathrm{M}\mathrm{P}\mathrm{a}\\ \frac{x}{d}\le 0.35if{f}_c>50MPa\end{array}\right.$ (27)

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 ($M_R\left(t\right)$) and load effect ($M_S\left(t\right)$) 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 (5).Cumulative probability of failure results, considering and ignoring corrosion,are shown in Figure (6). In Figure (6a), with no corrosion, failure probabilities increase in time only due to successive load applications. As observed, the effect of depassivation and corrosion is to sharply increase cumulative failure probabilities in time.

Figure 5
One sample of time-evolution of $M_R\left(t\right)$ and $M_S\left(t\right)$.

Figure 6
Results for cumulative probability of failure.

7. 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.

8. 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).

REFERENCES

• Al-Alaily, H.S.; Hassan A.A.A.; Hussein A.A. Use of extended finite element method and statistical analysis for modelling the corrosion-induced cracking in reinforced concrete containing metakaolin. Canadian Journal of Civil Engineering. v.45, 167-178, 2018.
• Andrieu-Renaud, C.; Sudret, B.; Lemaire, M., 2004. The PHI2 method: a way to compute timevariant reliability. Reliab. Eng. Sys. Safety 84, 75-86.
• Apostolopoulos, C.A.; Papadakis, V.G. Consequences of steel corrosion on the ductility properties of reinforcement bar. ConstructionandBuildingMaterials, v. 22, n. 12, p. 2316-2324, 2008.
• Audenaert, K.; Yuan, Q.; De Schutter, G. On the time dependency of the chloride migration coefficient in concrete. Construction and Building Materials v.24, 2010, p. 396-402.
• Bastidas-Arteaga, E.; Chateauneuf, A.; Sánchez-Silva, M.; Bressolette, P.; Schoefs, F. A comprehensive probabilistic model for chloride ingress in unsaturated concrete. Engineering Structures, v. 51, p. 259-266, 2011.
• Bastidas-Arteaga, E.; Stewart, M.G. Damage risks and economic assessment of climate adaptation strategies for design of new concrete structures subject to chloride-induced corrosion. Structural Safety, v.52, p.40- 53, 2015.
• Biondini, F.; Bontempi, F.; Frangopol, D.M.; Malerba, P.G. Probabilistic Service Life Assessment and Maintenance Planning of Concrete Structures. Journal of Structural Engineering, v.132, n.5, p.810-825, 2006.
• Brazilian Association of Technical Standards, 2014. Design of concrete structures: procedure (NBR 6118) (in portuguese: Associação Brasileira de Normas Técnicas, 2014. Projeto de estruturas de concreto: procedimento), Rio de Janeiro, 238 p.
• Carslaw, H.S.; Jaeger, J.C. Conduction of heat in solids, 2nd edn, Clarendon Press, Oxford, 1959.
• Cascudo, O. The control of corrosion in concrete reinforcing bars: inspection and electrochemical techniques. (in portuguese: O controle da corrosão de armaduras em concreto: inspeção e técnicas eletroquímicas.) São Paulo: Pini; Goiania: Editora UFMG, 1997.
• Comité Euro-Internacional du Béton - Durable concrete structures - design guide (CEB), 1992. Bulletin d’information n. 183, London, Thomas Telford.
• Curran, D.; Cross, M.; Lewis, B.A. A preliminary analysis of boundary element methods applied to parabolic partial differential equations. New Developments in Boundary Element Methods, Computational Mechanics Publications, Southampton, 1980.
• Dang, V.H.; François, R. Influence of long-term corrosion in chloride environment on mechanical behaviour of RC beam. Engineering Structures, v.48, p.558-568, 2013.
• Duddu, R. Numerical modeling of corrosion pit propagation using the combined extended finite element and level set method. Computational Mechanics, v. 54, 613-627, 2014.
• Duprat, F. Reliability of RC beams under chloride-ingress. Construction and Building Materials, v.21, p.1605- 1616, 2007.
• El Hassan, J.; Bressolette, P.; Chateauneuf, A.; El Tawil, K. Reliability based assessment of the effect of climatic conditions on the corrosion of RC structures subject to chloride ingress. Engineering Structures, v.32, p. 3279-3287, 2010.
• Enright, M.P.; Frangopol, D.M. Maintenance planning for deteriorating concrete bridges. Journal of Structural Engineering, v.125, n.12, p.1407-1414, 1999
• Enright, M.P.; Frangopol, D.M. Probabilistic analysis of resistance degradation of reinforced concrete bridge beams under corrosion. Engineering Structures, v.20, n.11, p.960-971, 1998
• Ghosh, J.; Padgett, J.E. Aging considerations in the development of time-dependent seismic fragility curves. Journal of Structural Engineering, v.136(12), p.1497-1511, 2010.
• Gonzalez, J.A.; Andrade, C.; Alonso, C.; Feliu, S. Comparison of rates of general corrosion and maximum pitting penetration on concrete embedded steel reinforcement. Cement and Concrete Research, v. 25, 1995, p. 257-264.
• Guzmán, S., Gálvez, J.C., Sancho, J.M. Cover cracking of reinforced concrete due to rebar corrosion induced by chloride penetration. Cement and Concrete Research, v.41, 2011, p. 893-902.
• Khelifa, M-R. New Computational Model Based on Finite Element Method to Quantify Damage Evolution Due to External Sulfate Attack on Self-Compacting Concretes. Computer-Aided Civil and Infrastructure Engineering. v.28, 260-272, 2013.
• Kroetz, H. M.; Moustapha, M.; Beck, A. T.; Sudret, B. A Two-Level Kriging-Based Approach with Active Learning for Solving Time-Variant Risk Optimization Problems. Reliability Engineering & System Safety. 203, 2020, 107033
• Kroetz, H. M.; Tessari, R. K.; Beck, A. T. . Kriging-assisted Risk Optimization with Ranked Weighted Average Simulation. In: 12th International Conference on Structural Safety & Reliability ICOSSAR, 2017a, Viena. Safety, Reliability, Risk, Resilience and Sustainability of Structures and Infrastructure, 2017a. p. 217-226.
• Kroetz, H. M.; Tessari, R.; Beck, A. Performance of global metamodeling techniques in solution of structural reliability problems. Adv. in Eng. Soft. 114, 2017b, 394-404.
• Li, C.; Der Kiureghian, A., 1993. Optimal discretization of random fields. J. Eng. Mech. 119(6), 1136-1154.
• Mangat, P.S.; Molloy, B.T. Prediction of long term chloride concentration in concrete. Materials and Structures, v.27, 1994, p. 338-346.
• Marelli, S.; Sudret, B., 2014. UQLab: a framework for uncertainty quantification in Matlab. In: Vulnerability, uncertainty and risk, proceedings of the 2nd international conference on vulnerability, risk and analysis management (ICVRAM2014), Liverpool, PP 2554-2563
• Mori, Y.; Ellingwood, B.R. Maintaining Reliability of Concrete Structures. I: Role of Inspection/Repair. Journal of Structural Engineering, v.120, n.3, p.824-845, 1994
• Moustapha, M.; Sudret, B.; Bourinet, J.-M.; Guillaume, B., 2016. Quantile-based optimization under uncertainties using adaptive Kriging surrogate models. Struct. Multidisc. Optim. 54(6), 1403-1421.
• Pan, T. ; Lu, Y. Stochastic modeling of reinforced concrete cracking due to nonuniform corrosion: FEM-based cross-scale analysis. Journal of Materials in Civil Engineering, June, Vol.24(6), p.698(9), 2012.
• Pellizzer, G.P.; Leonel, E.D.; Nogueira, C.G. Numerical approach about the effect of the corrosion on the mechanical capacity of the reinforced concrete beams considering material nonlinear models. Ibracon Structures and Materials Journal, v.11, n.7, p.26-51, 2018.
• Piasecka, B.A. Interval boundary element method for 2D transient diffusion problem using the directed interval arithmetic. Engineering Analysis with Boundary Elements, Vol.35(3), pp.259-263, 2011.
• Samson, E.; Marchand, J.; Snyder, K.A. Calculation of ionic diffusion coefficients on the basis of migration test results. Mater Struct, v.36, 2003, p. 156-165.
• Simon, J.; Bracci, J.M.; Gardoni, P. Seismic response and fragility of deteriorated reinforced concrete bridges. Journal of Structural Engineering, v.136(10), p.1273-1281, 2010.
• Sørensen, H.E. Chloridtransportihaerdnetbeton - mekanismer, modellerogmalemetoder (Transport of chloride in hardened concrete - mechanisms, models and test methods). Ph.D thesis, 1996.
• Stewart, M.G. Spatial variability of pitting corrosion and its influence on structural fragility and reliability of RC beams in flexure. Structural Safety, v.26, p.453-470, 2004.
• Stewart, M.G.; Wang, X.; Nguyen, M.N. Climate change adaptation for corrosion control of concrete infrastructure. Structural Safety, v.35, p.29-39, 2012.
• Sudret, B., 2008. Analytical derivation of the outcrossing rate in time-variant reliability problems. Struct. Infrastruct. Eng 4(5), 353-362.
• Sudret, B.; Der Kiureghian, A., 2000. Stochastic finite elements and reliability: a state-of-the-art report. Technical Report UCB/SEMM-2000/08, University of California, Berkeley (173 pages).
• Suo, Q.; Stewart, M.G. Corrosion cracking prediction updating of deteriorating RC structures using inspection information. Reliability Engineering and System Safety, v. 94, p. 1340-1348, 2009.
• Thoft-Christensen, P. Assessment of the reliability profiles for concrete bridges. Engineering Structures, v.20, n.11, p.1004-1009, 1998.
• Thoft-Christensen, P.; Hansen, H.I., 1994. Optimal strategy for maintenance of concrete bridges using expert systems. Proc., ICOSSAR'93, A.A.Balkema, Rotterdam, The Netherlands, p. 939-946.
• Tuutti, K., 1982. Corrosion of steel in concrete Swedish. Cement and concrete research institute Stockholm, 460 p.
• Val, D.V., Chernin, L., Stewart, M.G. Experimental and numerical investigation of corrosion-induced cover cracking in reinforced concrete structures. Journal of Structural Engineering, ASCE, v.135, 2009, p. 376-385.
• Val, D.V.; Melchers, R.E. Reliability of deteriorating RC slab bridges. Journal of Structural Engineering (ASCE), v. 123, n. 12, p. 1638-1644, 1997.
• Val, D.V.; Stewart, M.G. Life-cycle cost analysis of reinforced concrete structures in marine environments. Structural Safety, v. 25, n. 4, p. 343-362, 2003.
• Val, D.V.; Trapper, P.A. Probabilistic evaluation of initiation time of chloride-induced corrosion. Reliability Engineering and System Safety, v.93, p.364-372, 2008.
• Vidal, T.; Castel, A.; François, R. Corrosion process and structural performance of a 17 year old reinforced concrete beam stored in chloride environment. Cement and Concrete Research, v.37, p.1551-1561, 2007.
• Vu, K.A.T.; Stewart, M.G. Structural reliability of concrete bridges including improved chloride-induced corrosion models. Structural Safety, v.22, n.4, 2000, p. 313-333.
• Wang, B.; Chen, X. Fast multipole accelerated dual reciprocity boundary element method for diffusion problems. Annals of Solid and Structural Mechanics, Vol.7(1), pp.45-58, 2015.
• Warkus, J.; Brem, M.; Raupach, M. BEM-models for the propagation period of chloride induced reinforcement corrosion. Materials and Corrosion, v. 57, 636-641, 2006.
• Weyers, R.E. Corrosion service life model. In: Silva-Araya W.P., de Rincon O.T., O’Neill L.P., editors. Repair and rehabilitation of reinforced concrete structures: the state of the art. Reston, VA: American Society of Civil Engineers, 105 p., 1998.
• Wrobel, L.C. The Boundary Element Method Volume 1: Applications in Thermo-Fluids and Acoustics. 1 ed., Wiley, 2002, 468 p.
• Wrobel, L.C.; Brebbia, C.A. A formulation of the boundary element method for axisymmetric transient heat conduction, International Journal of Heat and Mass Transfer, v.24, 843-850, 1981.
• Wrobel, L.C.; Brebbia, C.A. The dual reciprocity boundary element formulation for nonlinear diffusion problems. Computer Methods in Applied Mechanics and Engineering, v.65, 147-164, 1987.
• Xiang, T.; Zhao, R. Reliability evaluation of chloride diffusion in fatigue damaged concrete. Engineering Structures, v.29, p.1539-1547, 2007.
• Xiao, J.; Ying, J.; Shen, L. FEM simulation of chloride diffusion in modeled recycled aggregate concrete, Construction Building Materials. v.29, 12-23, 2012.
• Yang, L.F.; Chen, Z.; Gao, Q.; Ju, J.W. Compensation length of two-dimensional chloride diffusion in concrete using a boundary element model. ActaMechanica, v. 224, 123-137, 2013.
• Yoon, S.; Reis, J.M.L. Simulation of Chloride Ingress through Surface-Coated Concrete during Migration Test Using Finite-Difference and Finite-Element Method. International Journal of Polymer Science, Vol.2017, 12 pages, 2017.
• Zhang, R.; Castel, A.; François, R. Concrete cover cracking with reinforcement corrosion of RC beam during chloride-induced corrosion process. Cementand Concrete Research, v. 40, p. 415-425, 2010.

Appendix A: Singularity Subtraction Method (SSM). Analytical Solutions for Transient Diffusion Case

The singularity subtraction method (SSM) has been applied to regularize the kernel ${\mathrm{u}}^{\mathrm{*}}$. The integration by parts of kernel ${\mathrm{u}}^{\mathrm{*}}$ as a function of time leads to the following:

$U_1^{*}\left(\underset{\_}{\xi },\underset{\_}{x}\right)=\underset{t_0^k}{\overset{t_f^k}{{\displaystyle \int }}}{u}^{*}\left(\underset{\_}{\xi },\underset{\_}{x},t_F,t\right)dt=\frac{1}{4\pi \kappa }\left\{E_1\left[\frac{r^2}{4\kappa \left(t_F-t_0^k\right)}\right]-E_1\left[\frac{r^2}{4\kappa \left(t_F-t_f^k\right)}\right]\right\}=\frac{1}{4\pi \kappa }\left\{E_1\left(a_0^k\right)-E_1\left(a_f^k\right)\right\}$ (A.1)

in which ${\mathrm{E}}_1$ is the exponential-integral function ${\mathrm{E}}_1\left(\mathrm{z}\right)={\int }_{\mathrm{z}}^{\mathrm{\infty }}\frac{{\mathrm{e}}^{-\mathrm{t}}}{\mathrm{t}}\mathrm{d}\mathrm{t}$ and ${\mathrm{a}}_{\mathrm{f}}^{\mathrm{k}}$ and ${\mathrm{a}}_0^{\mathrm{k}}$ are as follows:

${\mathrm{a}}_{\mathrm{f}}^{\mathrm{k}}=\frac{{\mathrm{r}}^2}{4\mathrm{\kappa }\left({\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_{\mathrm{f}}^{\mathrm{k}}\right)}$ (A.2)

${\mathrm{a}}_0^{\mathrm{k}}=\frac{{\mathrm{r}}^2}{4\mathrm{\kappa }\left({\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_0^{\mathrm{k}}\right)}$ (A.3)

The singularity is observed in ${\mathrm{U}}_{\mathrm{k}}^{\mathrm{*}}$ 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:

$U_1^{*}\left(\underset{\_}{\xi },\underset{\_}{x}\right)=\frac{1}{4\pi \kappa }\left[E_1\left(a_0^1\right)-E_1\left(a_f^1\right)\right]$ (A.4)

The variables ${\mathrm{a}}_{\mathrm{f}}^{\mathrm{k}}$ and ${\mathrm{a}}_0^{\mathrm{k}}$ are positive because $\mathrm{\kappa }>0$, $\mathrm{r}\ge 0$ and ${\mathrm{t}}_{\mathrm{F}}\ge {\mathrm{t}}_{\mathrm{f}}^{\mathrm{k}}>{\mathrm{t}}_0^{\mathrm{k}}$. In Eq. (A.2), when $\mathrm{r}$ is nil and ${\mathrm{t}}_{\mathrm{f}}^1\to {\mathrm{t}}_{\mathrm{F}}$, ${\mathrm{a}}_{\mathrm{f}}^1\to +\mathrm{\infty }$. Thus, the exponential-integral function is evaluated at + ∞, which results into a nil value. In equation (A.3) when r is nil, ${\mathrm{a}}_0^1\to 0$, because the denominator (${\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_0^1$) is always greater than zero (${\mathrm{t}}_{\mathrm{F}}\ge {\mathrm{t}}_{\mathrm{f}}^1>{\mathrm{t}}_0^1$). Thus, the exponential-integral function has nil value, which results into the following singularity:

$U_1^{*}\left(\underset{\_}{\xi },\underset{\_}{x}\right)=\frac{1}{4\pi \kappa }\left[E_1\left(0\right)-E_1\left(+\infty \right)\right]=\frac{1}{4\pi \kappa }\left[E_1\left(0\right)\right]$ (A.5)

The singularity arising in this equation is logarithmic for distinct intervals of the exponential-integral function argument, $0\le \mathrm{z}\le 1$and $1<z<\infty$. Thus, the kernel containing ${\mathrm{u}}^{\mathrm{*}}$ can be regularized by the following equation:

$\underset{{\Gamma }_j}{\overset{}{{\displaystyle \int }}}{U}_1^{*}\left(\underset{\_}{\xi },\underset{\_}{x}\right)d{\Gamma }_j=\frac{1}{4\pi \kappa }\underset{{\Gamma }_j}{\overset{}{{\displaystyle \int }}}\left[E_1\left(a_0^1\right)-\ln \left(a_0^1{}^{*}\right)\right]d{\Gamma }_j+\frac{1}{4\pi \kappa }\underset{{\Gamma }_j}{\overset{}{{\displaystyle \int }}}\ln \left(a_0^1{}^{*}\right)d{\Gamma }_j$ (A.6)

where${\mathrm{\xi }}_0$ is the dimensionless coordinate of the source point, $\mathrm{\xi }$the dimensionless coordinate of the field point, ${\mathrm{x}}_{\mathrm{i}}$ the real coordinate of a given point, $\mathrm{\varphi }$ the shape function, ${\mathrm{\varphi }}_{,\mathrm{i}}$ the shape function derivative and the Jacobian is given by: $\mathrm{J}=\sqrt{{\left[{\mathrm{\varphi }}_{,\mathrm{i}}\left({\mathrm{\xi }}_0\right){\mathrm{x}}_{\mathrm{i}}\right]}^2}$, the real distance $\mathrm{r}$ between the source and the field points is evaluated as follows:

$\mathrm{r}={\mathrm{r}}^{\mathrm{*}}=\left|\mathrm{J}\right|\mathrm{\epsilon }=\left|\mathrm{J}\right|\left|\mathrm{\xi }-{\mathrm{\xi }}_0\right|$ (A.7)

where $\mathrm{\epsilon }=\left|\mathrm{\xi }-{\mathrm{\xi }}_0\right|$ is the dimensionless distance between the source and field points. Substituting $\mathrm{d}\mathrm{\Gamma }=\mathrm{J}\mathrm{d}\mathrm{\xi }$ in Eq. (A.6) one has:

$\underset{{\Gamma }_j}{\overset{}{{\displaystyle \int }}}{U}_1^{*}\left(\underset{\_}{\xi },\underset{\_}{x}\right)d{\Gamma }_j=\frac{1}{4\pi \kappa }\underset{-1}{\overset{1}{{\displaystyle \int }}}{E}_1\left(a_0^1\right)\varphi \left(\xi \right)J\left(\xi \right)d\xi -\frac{1}{4\pi \kappa }\underset{-1}{\overset{1}{{\displaystyle \int }}}\ln \left(a_0^1{}^{*}\right)\varphi \left({\xi }_0\right)J\left({\xi }_0\right)d\xi +\frac{1}{4\pi \kappa }\underset{-1}{\overset{1}{{\displaystyle \int }}}\ln \left(a_0^1{}^{*}\right)\varphi \left({\xi }_0\right)J\left({\xi }_0\right)d\xi$ (A.8)

in which:

${{\mathrm{a}}_0^1}^{\mathrm{*}}=\frac{{\left({\mathrm{r}}^{\mathrm{*}}\right)}^2}{4\mathrm{\kappa }({\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_0^1)}=\frac{{\left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)}^2}{4\mathrm{\kappa }({\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_0^1)}$ (A.9)

In Eq. (A.8), the two first integrals on the right side are bounded and no longer singular. Then, they are assessed properly by the Gauss-Legendre quadrature. The last integral is evaluated analytically. The latter term on that equation is named as $\mathrm{I}$:

$\mathrm{I}=\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1}^1\mathrm{l}\mathrm{n}\left({{\mathrm{a}}_0^1}^{\mathrm{*}}\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\xi }$ (A.10)

The change of the integration domain from $\mathrm{d}\mathrm{\xi }$ to $\mathrm{d}\mathrm{\epsilon }$ and consequently the integration intervals lead to the following:

$\mathrm{I}=\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1-{\mathrm{\xi }}_0}^{1-{\mathrm{\xi }}_0}\mathrm{l}\mathrm{n}\left[\frac{{\left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)}^2}{4\mathrm{\kappa }({\mathrm{t}}_{\mathrm{F}}-{\mathrm{t}}_0^1)}\right]\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }$ (A.11)

For simplicity, considering ${\mathrm{t}}_0^1=0$, one obtains:

$\mathrm{I}=\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1-{\mathrm{\xi }}_0}^{1-{\mathrm{\xi }}_0}2\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }-\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1-{\mathrm{\xi }}_0}^{1-{\mathrm{\xi }}_0}\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }$ (A.12)

Because of the singular kernel contained in Eq. (A.12), the integral must be evaluated in the Cauchy principal (CPV) sense. Then:

$\mathrm{C}\mathrm{P}\mathrm{V}=\underset{\mathrm{\epsilon }\to 0}{\lim }\left\{\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1-{\mathrm{\xi }}_0}^{-\mathrm{\epsilon }}2\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }+\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{+\mathrm{\epsilon }}^{1-{\mathrm{\xi }}_0}2\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }-\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{-1-{\mathrm{\xi }}_0}^{-\mathrm{\epsilon }}\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }-\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}{\int }_{+\mathrm{\epsilon }}^{1-{\mathrm{\xi }}_0}\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }\right\}\to$ $\mathrm{C}\mathrm{P}\mathrm{V}=\underset{\mathrm{\epsilon }\to 0}{\lim }\left\{\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\left[{\int }_{-1-{\mathrm{\xi }}_0}^{-\mathrm{\epsilon }}2\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)\mathrm{d}\mathrm{\epsilon }+{\int }_{+\mathrm{\epsilon }}^{1-{\mathrm{\xi }}_0}2\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\left|\mathrm{\epsilon }\right|\right)\mathrm{d}\mathrm{\epsilon }-{\int }_{-1-{\mathrm{\xi }}_0}^{-\mathrm{\epsilon }}\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{d}\mathrm{\epsilon }-{\int }_{+\mathrm{\epsilon }}^{1-{\mathrm{\xi }}_0}\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\mathrm{d}\mathrm{\epsilon }\right]\right\}$

${\mathrm{I}}_1$ (A.13)

To evaluate the terms of Eq. (A.13) properly, they are rewritten in separated form, from ${\mathrm{I}}_1$ to ${\mathrm{I}}_4$ as follows:

$\mathrm{C}\mathrm{P}\mathrm{V}=\underset{\mathrm{\epsilon }\to 0}{\lim }\left\{\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\left[{\mathrm{I}}_1+{\mathrm{I}}_2-({\mathrm{I}}_3+{\mathrm{I}}_4)\right]\right\}$ (A.14)

Thus, evaluating ${\mathrm{I}}_1$one has:

${\mathrm{I}}_1=2{\left[\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)\right]}_{-1-{\mathrm{\xi }}_0}^{-\mathrm{\epsilon }}=2\left\{-\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)+\mathrm{\epsilon }-\left[\left(-1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)(-1-{\mathrm{\xi }}_0)\right|\right)+1+{\mathrm{\xi }}_0\right]\right\}$ (A.15)

Thus, evaluating ${\mathrm{I}}_2$one has:

${\mathrm{I}}_2=2{\left[\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)\right]}_{\mathrm{\epsilon }}^{1-{\mathrm{\xi }}_0}=2\left\{\left(1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(1-{\mathrm{\xi }}_0\right)\right|\right)-1+{\mathrm{\xi }}_0-\left[\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)-\mathrm{\epsilon }\right]\right\}$ (A.16)

Then, evaluating ${\mathrm{I}}_3$one has:

${\mathrm{I}}_3=\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)[-\mathrm{\epsilon }+1+{\mathrm{\xi }}_0]$ (A.17)

Then, evaluating ${\mathrm{I}}_4$ one has:

${\mathrm{I}}_4=\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)[1-{\mathrm{\xi }}_0-\mathrm{\epsilon }]$ (A.18)

Adding ${\mathrm{I}}_1$ and ${\mathrm{I}}_2$ one obtains:

${\mathrm{I}}_1+{\mathrm{I}}_2=2\left\{-\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)+\mathrm{\epsilon }+\left(1+{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(-1-{\mathrm{\xi }}_0\right)\right|\right)-1-{\mathrm{\xi }}_0+\left(1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(1-{\mathrm{\xi }}_0\right)\right|\right)-1+{\mathrm{\xi }}_0-\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)+\mathrm{\epsilon }\right\}\to$ ${\mathrm{I}}_1+{\mathrm{I}}_2=2\left\{-2\mathrm{\epsilon }\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right|\right)+2\mathrm{\epsilon }-2+\left(1+{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(-1-{\mathrm{\xi }}_0\right)\right|\right)+\left(1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(1-{\mathrm{\xi }}_0\right)\right|\right)\right\}$ (A.19)

Adding ${\mathrm{I}}_3$ and ${\mathrm{I}}_4$ one obtains:

${\mathrm{I}}_3+{\mathrm{I}}_4=\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\left[-\mathrm{\epsilon }+1+{\mathrm{\xi }}_0+1-{\mathrm{\xi }}_0-\mathrm{\epsilon }\right]=\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\left[-2\mathrm{\epsilon }+2\right]$ (A.20)

By applying the L'Hôpital’s theorem one obtains:

$\underset{\mathrm{\epsilon }\to 0}{\mathrm{l}\mathrm{i}\mathrm{m}}\left\{\mathrm{\epsilon }\ln \left(\mathrm{J}\left({\mathrm{\xi }}_0\right)\mathrm{\epsilon }\right)\right\}=0$ (A.21)

Therefore, the limit of Eq. (A.14) is achieved as follows:

$\mathrm{C}\mathrm{P}\mathrm{V}=\frac{1}{4\mathrm{\pi }\mathrm{\kappa }}\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\left\{-4+2\left(1+{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(-1-{\mathrm{\xi }}_0\right)\right|\right)+2\left(1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(1-{\mathrm{\xi }}_0\right)\right|\right)-2\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\right\}$ (A.22)

Consequently:

$\mathrm{C}\mathrm{P}\mathrm{V}=\frac{1}{2\mathrm{\pi }\mathrm{\kappa }}\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\left\{-2+\left(1+{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(-1-{\mathrm{\xi }}_0\right)\right|\right)+\left(1-{\mathrm{\xi }}_0\right)\ln \left(\left|\mathrm{J}\left({\mathrm{\xi }}_0\right)\left(1-{\mathrm{\xi }}_0\right)\right|\right)-\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\right\}$ (A.23)

It is worth mentioning that Eq. (A.23) is valid for discontinuous elements, i.e., when the source points are not positioned at the element ends ${\mathrm{\xi }}_0=\pm 1$. When continuous elements are utilized, the finite part of Eq. (A.23) is evaluated, leading to the following:

$\mathrm{C}\mathrm{P}\mathrm{V}=\frac{1}{2\mathrm{\pi }\mathrm{\kappa }}\mathrm{\varphi }\left({\mathrm{\xi }}_0\right)\mathrm{J}\left({\mathrm{\xi }}_0\right)\left\{-2+2\ln \left(\left|2\mathrm{J}\left({\mathrm{\xi }}_0\right)\right|\right)-\ln \left(4\mathrm{\kappa }{\mathrm{t}}_{\mathrm{F}}\right)\right\}$ (A.24)

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  11 Dec 2020
• Date of issue
  2020

History

• Received
  02 Dec 2019
• Reviewed
  24 June 2020
• Accepted
  31 Aug 2020
• Published
  02 Sept 2020