An X–FEM Technique for Modeling the FRP Strengthening of Concrete Arches with a Plastic–Damage Model; Numerical and Experimental Investigations

Amir R. Khoei^1,*, Tahmaz Ahmadpour² and Yousef Navidtehrani¹

¹Center of Excellence in Structures and Earthquake Engineering, Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran

²Department of Civil Engineering, School of Science and Engineering, Sharif University of Technology, International Campus, P.O. Box 76417-76655, Kish Island, Iran

E-mail: arkhoei@sharif.edu

${}^{*}$Corresponding Author

Received 25 February 2021; Accepted 11 May 2021; Publication 13 July 2021

Abstract

In this paper, an enriched–FEM method is presented based on the X-FEM technique by applying a damage–plasticity model to investigate the effect of FRP strengthening on the concrete arch. In this manner, the damage strain is lumped into the crack interface while the elastic and plastic strains are employed within the bulk volume of element. The damage stress–strain relation is converted to the traction separation law using an acoustic tensor. The interface between the FRP and concrete is modeled using a cohesive fracture model. The X-FEM technique is applied where the FE mesh is not necessary to be conformed to the fracture geometry, so the regular mesh is utilized independent of the position of the fracture. The accuracy of the proposed plastic-damage model is investigated under the monotonic tension, compression, and cyclic tension loading. Furthermore, the accuracy of the cohesive fracture model is investigated using the experimental data reported for the debonding test. In order to verify the accuracy of the proposed computational algorithm, the numerical results are compared with those of experimental data obtained from two tests conducted on reinforced concrete arches strengthened with FRP. Finally, a parametric study is performed by evaluating the effects of high to span ratio, longitudinal reinforcement ratio, and strengthening method.

Keywords: Concrete arch, FRP retrofitting, X-FEM method, plastic-damage model, Cohesive fracture model.

1 Introduction

The FRP strengthening of concrete is one the most effective retrofitting methods that include the high durability against environmental degradation factors, high tensile strength, simple implementation, low weight, and easy transportation. A number of research works have been performed during last two decades to investigate the structural retrofitting with the FRP; which were focused on the beams, columns, flat slabs, and masonry structures, while a few studies was performed on the behavior of FRP strengthening of concrete arches. Since the arch members are utilized in many structures such as bridges, fluid storage tanks, tunnels, and domes, a study of the methods for retrofitting these structures has gained interest, recently. There are several studies reported in the literature on the masonry arches retrofitted with FRP materials that present a significant improvement in structural performance [1–3]. Dagher et al. [4] studied the bending behavior of concrete-filled tubular FRP arches in bridge structures. The effect of FRP strengthening of concrete arch structures was investigated experimentally by Chen et al. [5], in which the response of concrete arch structures retrofitted by wrapping FRP carbon sheets was studied under the explosive impulses. Hamed et al. [6] studied the performance of the concrete arch retrofitted with externally bonded composite materials, and observed that the maximum load increases about 40%. Zhang et al. [7] performed a set of experimental tests to improve the strength, stiffness and ductility of reinforced concrete arches retrofitted with FRP materials.

Basically, characteristics of the concrete plastic behavior cannot be described by the classical theory of plasticity. Many attempts have been made to demonstrate the behavior of concrete with the theory of plasticity; e.g. Pramono and Willam [8], Feenstra and de Brost [9], Li et al. [10] and Červenka and Papanikolaou [11]. The concrete behavior was also described by the Menetrey-William failure model to model crushing under high confinement, cracks in concrete, and closure due to crushing in different directions. Basically, a constitutive model that is able to predict the behavior of concrete from the beginning of loading to the final failure stage is significant. The damage happens in the concrete due to micro-cracks; in fact, micro-cracks occur because of thermal expansion at the interface between the cement and aggregate. Several studies have been carreied out to model the damage in concrete structures, including: a nonlocal micro-plane model for the fracture and damage by Bažant and Ožbolt [12] and a damage model based on the bounding surface concept by Voyiadjis and Abu-Lebdeh [13]. The combined plastic–damage models have been used by researchers to model the stiffness degradation of concrete, in which the damage variables are assumed according to the plastic deformation in a constitutive model to calibrate parameters with experimental data. In the plastic–damage model, the irreversible plastic behavior was modeled using the concept of plasticity in an effective stress space, while the stiffness degradation was captured using the continuum damage mechanics; see Lubliner et al. [14], Lee and Fenves [15], Faria et al. [16], Ibrahimbegović et al. [17], Salari et al. [18], Gatuingt et al. [19], Yu et al. [20]. Since the behavior of concrete is different in tension and compression, two damage models are utilized for each type of loading; however, some studies only consider one damage variable for different loading conditions [21–23]. However, it is significant to model the concrete damage using two separate damage variables for tension and compression [24, 25].

The plastic–damage models can be generally used to predict the concrete behavior up to the onset of fracture. These approaches are prone to localization and ill-posedness issues that need proper regularization techniques. In order to overcome these drawbacks it is required to introduce discontinuities into the domain. Adaptive finite element technique is one of the most popular approaches used extensively to model fracture problems [26–28]. However, in this method the FE mesh must be confined to the crack body; so the crack growth simulation is time consuming and costly. An alternative technique that alleviates this problem is based on the partition of unity method [29] in which the Extended–FEM (X-FEM) technique is the most popular one [30]. The technique is utilized to enrich the approximation space by adding appropriate enrichment functions through a partition of unity method, where the singularities, high gradients and discontinuities can be achieved [31]. The most research works on the X-FEM with plasticity model is mainly focused on the crack-tip analysis. Broumand and Khoei [32] performed the crack growth simulation in large deformation ductile fracture problem using the Lemaitre damage-plasticity model in the framework of X-FEM. Pañeda et al. [33] proposed a gradient-enhanced numerical framework to enrich the displacement field of the X-FEM with the stress singularity of the gradient-dominated solution. Although there are several research works reported on the X-FEM analysis of plasticity problems to model the ductile fracture behavior, little attention has been paid on modeling the concrete behavior with damage-plastic model within the X-FEM. In this study, a novel technique is presented to model the damage-plastic model for concrete with strong discontinuity via the X-FEM.

In the present paper, the X-FEM method is employed together with a damage–plasticity model to investigate the effect of FRP strengthening on concrete arches. The strain tensor is decomposed into three components, including the elastic, plastic and damage strains, in which the damage strain is employed to model the discontinuity at the crack interface. The procedure to convert the stress – damage strain relation to a traction separation law is presented. The interface between the FRP and concrete is modeled using a cohesive fracture model. Both the experimental and numerical investigations are applied to obtain the effects of FRP strengthening on the behavior of concrete arches. In order to present the capability of the proposed computational model, the experimental data of two experimental tests conducted on the arches strengthened with FRP are utilized. A parametric study is then performed for the arches with internal FRP, external FRP, and both face strengthening at different height to span ratios, longitudinal reinforcement ratios, and strengthening methods. It is shown that the arches with low high to span ratios and the internal FRP strengthening method is the most effective method, while by increasing the high to span ratio, the effectiveness of the external FRP method increases. It is also observed that the main drawback of the internal FRP strengthening method is debonding of FRP from the concrete.

2 The X-FEM Governing Equations

Consider a cracked body $\mathrm{\Omega }$ bounded externally by surface $\mathrm{\Gamma }$. The crack interfaces are ${\mathrm{\Gamma }}_c^{+}$ and ${\mathrm{\Gamma }}_c^{-}$ corresponding to the either sides of the crack. The equation of motion of the body can be expressed as

$$\mathbf{\nabla }\cdot \sigma +\mathbf{b}-\rho \ddot{\mathbf{u}}=0$$

(1)

where $\mathbf{\nabla }$ is the gradient operator, $\sigma$ is the Cauchy stress tensor, and $\mathbf{b}$ is the body force vector. The constitutive relation is defined as $\mathit{\sigma }=\mathbf{D}\mathit{\epsilon }$, where $\mathbf{D}$ is the tangential constitutive matrix of the material. The essential and natural boundary conditions are respectively defined as $\mathbf{u}=\tilde{\mathbf{u}}$ on ${\mathrm{\Gamma }}_u$, $\sigma \cdot {\mathbf{n}}_{\mathrm{\Gamma }}=\overline{\mathbf{t}}$ on ${\mathrm{\Gamma }}_t$, and $\sigma \cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d}=\mathbf{0}$ on ${\mathrm{\Gamma }}_d$, in which ${\mathbf{n}}_{\mathrm{\Gamma }}$ is the outward unit normal vector to the external boundary $\mathrm{\Gamma }$, and ${\mathbf{n}}_{{\mathrm{\Gamma }}_d}$ is the unit normal vector to the discontinuity ${\mathrm{\Gamma }}_d$. In addition, $\tilde{\mathbf{u}}$ is the prescribed displacement at the boundary ${\mathrm{\Gamma }}_u$, and $\overline{\mathbf{t}}$ is the prescribed traction at the boundary ${\mathrm{\Gamma }}_t$.

2.1 The Enrichment of Displacement Field

The X-FEM method was originally introduced by Belytschko and Black [30], and it was then modified by researchers to improve its convergence, such as the enhanced strain method by Borja [34], the cutoff function by Chahine et al. [35], the blending weight function by Tarancon et al. [36] and Fries [37], and the discontinuous Galerkin method by Gracie et al. [38] and Shen and Lew [39]. In order to model the crack discontinuity within the X-FEM, the nodal points of an element cut by the discontinuity are enhanced by the Heaviside function $H(\mathit{x})$. The Heaviside enrichment function is generally employed to model the crack discontinuity because of different displacement fields on each side of the crack as

$$H(\mathit{x})=\{\begin{array}{cc}+1\mathit{\hspace{1em}}\hfill & \phi (\mathit{x})\ge 0\hfill \\ -1\mathit{\hspace{1em}}\hfill & \phi (\mathit{x})<0\hfill \end{array}$$

(2)

where $\phi (\mathit{x})$ is the signed distance function defined using the absolute value of level set function as

$$\phi (\mathit{x})=\min \parallel \mathit{x}-{\mathit{x}}^{*}\parallel \mathrm{sign}((\mathit{x}-{\mathit{x}}^{*})\cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d})$$

(3)

where ${\mathit{x}}^{*}$ is a point on the crack inteface which has the closest distance from the point $\mathit{x}$, and ${\mathbf{n}}_{{\mathrm{\Gamma }}_d}$ is the normal vector to the crack interface at point ${\mathit{x}}^{*}$. Because of independent displacement fields at both sides of the crack, the Heaviside function is a proper function in modeling the crack interface.

Thus, the displacement field can be written for the cracked element as

$\mathbf{u}(\mathit{x},t)={\displaystyle \sum _{I\in \mathcal{𝒩}}}{N}_I(\mathit{x}){\overline{\mathbf{u}}}_I+{\displaystyle \sum _{J\in {\mathcal{𝒩}}^{dis}}}{N}_J(\mathit{x})(H(\mathit{x})-H({\mathit{x}}_J)){\overline{\mathbf{a}}}_J$

(4)

where $\mathcal{𝒩}$ consists of all nodal points and ${\mathcal{𝒩}}^{dis}$ is the set of nodes enriched by the Heaviside function. In relation (4), ${\overline{\mathbf{u}}}_I$ are the unknown standard nodal DOFs at node $I$ and ${\overline{\mathbf{a}}}_J$ are the unknown enriched nodal DOFs associated with the Heaviside function at node $J$. Hence, the X-FEM displacement approximation (4) can be rewritten as

$$\mathbf{u}(\mathit{x},t)={\mathbf{N}}^{std}(\mathit{x})\overline{\mathbf{u}}(t)+{\mathbf{N}}^{Hev}(\mathit{x})\overline{\mathbf{a}}(t)$$

(5)

where ${\mathbf{N}}^{std}(\mathit{x})\equiv \mathbf{N}(\mathit{x})$ are the matrix of standard FEM shape functions and ${\mathbf{N}}^{enr}$ are the matrix of enriched shape functions associated with the Heaviside function.

The strain vector corresponding to the X-FEM displacement field (5) can be similarly defined according to the standard and enriched nodal values as

$$\mathit{\epsilon }(\mathit{x},t)={\mathbf{B}}^{std}(\mathbf{x})\overline{\mathbf{u}}(t)+{\mathbf{B}}^{Hev}(\mathbf{x})\overline{\mathbf{a}}(t)$$

(6)

where ${\mathbf{B}}^{std}$and ${\mathbf{B}}^{Hev}$ include the spatial derivatives of the standard and enriched shape functions associated with the Heaviside functions defined as

${\mathbf{B}}_I^{std}$	$=\left[\begin{array}{cc}\hfill \partial N_I\mathrm{/}\partial x\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \partial N_I\mathrm{/}\partial y\hfill \\ \hfill \partial N_I\mathrm{/}\partial y\mathit{\hspace{1em}}\hfill & \hfill \partial N_I\mathrm{/}\partial x\hfill \end{array}\right]$
${\mathbf{B}}_I^{Hev}$	$=\left[\begin{array}{cc}\hfill \partial \left(N_I\left(H(\mathit{x})-H({\mathit{x}}_I)\right)\right)\mathrm{/}\partial x\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \partial \left(N_I\left(H(\mathit{x})-H({\mathit{x}}_I)\right)\right)\mathrm{/}\partial y\hfill \\ \hfill \partial \left(N_I\left(H(\mathit{x})-H({\mathit{x}}_I)\right)\right)\mathrm{/}\partial y\mathit{\hspace{1em}}\hfill & \hfill \partial \left(N_I\left(H(\mathit{x})-H({\mathit{x}}_I)\right)\right)\mathrm{/}\partial x\hfill \end{array}\right]$	(7)

2.2 Discretization of Governing Equations

The weak form of equation (1) can be obtained using the Galerkin discretization technique by integrating the product of equation multiplied by admissible test function over the domain. The test function $\delta \mathbf{u}(\mathit{x},t)$ can be defined according to the displacement field (5) as

$$\delta \mathbf{u}(\mathit{x},t)={\mathbf{N}}^{std}(\mathit{x})\delta \overline{\mathbf{u}}(t)+{\mathbf{N}}^{Hev}(\mathit{x})\delta \overline{\mathbf{a}}(t)$$

(8)

The weak form of equation (1) can be derived by applying the Galerkin method as

$${\int }_{\mathrm{\Omega }}\delta \mathbf{u}(\mathit{x},t)(\mathbf{\nabla }\cdot \sigma +\mathbf{b}-\rho \ddot{\mathbf{u}})\mathit{d}\mathrm{\Omega }=0$$

(9)

The above integral equation can be rewritten by employing the discontinuous Divergence theorem, imposing the natural boundary conditions, and satisfying the traction free boundary condition on the crack interface as

${\displaystyle {\int }_{\mathrm{\Omega }}}\delta \mathbf{u}\cdot \rho \ddot{\mathbf{u}}\mathit{d}\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}\mathbf{\nabla }\delta \mathbf{u}:\mathit{\sigma }d\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_d}}\left[\left[\delta \mathbf{u}\right]\right]\mathit{\sigma }\cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\mathit{d}\mathrm{\Gamma }$
$\mathrm{\hspace{1em}\hspace{1em}}-{\displaystyle {\int }_{\mathrm{\Omega }}}\delta \mathbf{u}\cdot \mathbf{b}\mathit{d}\mathrm{\Omega }-{\displaystyle {\int }_{{\mathrm{\Gamma }}_t}}\delta \mathbf{u}\cdot \overline{\mathbf{t}}\mathit{d}\mathrm{\Gamma }=0$	(10)

where the symbol $\left[\left[\right]\right]$ indicates the jump across the crack interface; this jump can be obtained using the difference of displacement values at two crack surfaces, i.e. $\left[\left[\mathrm{\Xi }\right]\right]={\mathrm{\Xi }}^{+}-{\mathrm{\Xi }}^{-}$. Note that the third integral in the above equation can be removed along ${\mathrm{\Gamma }}_d$ by imposing the traction free boundary condition on the crack interfaces.

Finally, the discretized form of equation (10) can be obtained by removing the third integral as

$$\left(\begin{array}{cc}\hfill {\mathbf{M}}_{uu}\mathit{\hspace{1em}}\hfill & \hfill {\mathbf{M}}_{ua}\hfill \\ \hfill {\mathbf{M}}_{au}\mathit{\hspace{1em}}\hfill & \hfill {\mathbf{M}}_{aa}\hfill \end{array}\right)\left\{\begin{array}{c}\hfill \ddot{\overline{\mathbf{u}}}\hfill \\ \hfill \ddot{\overline{\mathbf{a}}}\hfill \end{array}\right\}+\left(\begin{array}{cc}\hfill {\mathbf{K}}_{uu}\mathit{\hspace{1em}}\hfill & \hfill {\mathbf{K}}_{ua}\hfill \\ \hfill {\mathbf{K}}_{au}\mathit{\hspace{1em}}\hfill & \hfill {\mathbf{K}}_{aa}\hfill \end{array}\right)\left\{\begin{array}{c}\hfill \overline{\mathbf{u}}\hfill \\ \hfill \overline{\mathbf{a}}\hfill \end{array}\right\}-\left\{\begin{array}{c}\hfill {\mathbf{f}}_u^{\mathrm{ext}}\hfill \\ \hfill {\mathbf{f}}_a^{\mathrm{ext}}\hfill \end{array}\right\}=0$$

(11)

where the mass matrix $\mathbf{M}$, stiffness matrix $\mathbf{K}$, and external force vectors ${\mathbf{f}}^{\mathrm{ext}}$ are defined as

${\mathbf{M}}_{\alpha \beta }$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{({\mathbf{N}}^{\alpha })}^T\rho {\mathbf{N}}^{\beta }\mathit{d}\mathrm{\Omega }$
${\mathbf{K}}_{\alpha \beta }$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{({\mathbf{B}}^{\alpha })}^T\mathbf{D}{\mathbf{B}}^{\beta }\mathit{d}\mathrm{\Omega }$
${\mathbf{f}}_{\alpha }^{\mathrm{ext}}$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{({\mathbf{N}}^{\alpha })}^T\rho \mathbf{b}\mathit{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_t}}{({\mathbf{N}}^{\alpha })}^T\overline{\mathbf{t}}\mathit{d}\mathrm{\Gamma }$	(12)

in which $(\alpha ,\beta )\in (std,Hev)$ denote the ‘standard’ and ‘Heaviside’ displacement fields.

3 The Concrete Plastic-Damage Model

In order to describe the behavior of concrete, the plastic-damage model originally presented by Lee and Fenves [15] and then proposed by Nguyen and Houlsby [24, 25] is employed; in which the plastic-damage model is developed on the basis of thermo-dynamical approach. The strain tensor $\epsilon$ is decomposed into the elastic part ${\epsilon }^e$ and the plastic part ${\epsilon }^p$, in which the relation between stress and strain is defined as

$$\sigma =E:(\epsilon -{\epsilon }^p)$$

(13)

where $E$ is the elasticity tensor, $\sigma$ is the stress and, $\epsilon$ and ${\epsilon }^p$ are the total strain and plastic strain, respectively. By mapping the stress into the effective stress space, the plasticity and damage equations can be solved independently. Applying the scalar damage variable $D$, the effective stress can be expressed as

$$\overline{\sigma }=(1-D)E_0:(\epsilon -{\epsilon }^p)$$

(14)

where $E_0$ is the undamaged elastic-stiffness modulus. Considering the non-associated flow rule, the plastic strain rate ${\dot{\epsilon }}^p$ can be defined as

$${\dot{\epsilon }}^p=\dot{\lambda }{\nabla }_{\overline{\sigma }}\mathrm{\Phi }(\overline{\sigma })$$

(15)

where $\dot{\lambda }$ denotes the plastic consistency parameter which is a non-negative function and $\mathrm{\Phi }$ is a scalar plastic potential function. Moreover, the damage variable $k$ is required to represent the damage states variable as

$$\dot{k}=\dot{\lambda }H(\overline{\sigma },k)$$

(16)

3.1 The Damage Model

In order to represent the tensile and compressive damages for concrete material two damage variables are defined. These damage variables have values between zero to one for covering the range from undamaged to completely damaged concrete. In the Barcelona model introduced by Lubliner et al. [14], the uniaxial stress is defined as a function of the plastic strain as

$$\sigma =f_0[(1+a)\exp (-b{\epsilon }^p)-a(\exp (-2{\epsilon }^p))]$$

(17)

where $a$ and $b$ are dimensionless constants and $f_0$ is the initial yield stress. Consider an exponential form for the degradation $D$ as

$$1-D=\exp (-c{\epsilon }^p)$$

(18)

The effective stress can be written as

$$\overline{\sigma }=f_0\left[\left(1+a\right){\left(\exp \left(-b{\epsilon }^p\right)\right)}^{1-\frac{c}{b}}-a{\left(\exp \left(-b{\epsilon }^p\right)\right)}^{2-\frac{c}{b}}\right]$$

(19)

in which the damage variable for the uniaxial loading is denoted by $k$, defined as

$$k=\frac{1}{g}{\int }_0^{{\epsilon }^p}\sigma ({\epsilon }^p)\mathit{d}{\epsilon }^p$$

(20)

where $g$ is defined as

$$g={\int }_0^{\mathrm{\infty }}\sigma ({\epsilon }^p)\mathit{d}{\epsilon }^p$$

(21)

which is the dissipated energy density during the forming of micro-cracking. Substituting equation (17) into (21), the relation between $g$, $a$, and $b$ can be defined as

$$g=\frac{f_0}{b}\left(1+\frac{a}{2}\right)$$

(22)

The uniaxial stress can be defined in term of the damage variable $k$ using (20) and (22) as

$$\sigma =\frac{f_0}{a}\left[\left(1+a\right)\sqrt{1+a\left(2+a\right)k}-1+a\left(2+a\right)k\right]$$

(23)

Moreover, the effective stress can be expressed in terms of damage variable as

$$\overline{\sigma }=\frac{f_0}{a}{\left[\left(1+a\right)-\sqrt{1+a\left(2+a\right)k}\right]}^{\mathrm{\hspace{0.33em}1}-\frac{c}{b}}\left(\sqrt{1+a\left(2+a\right)k}\right)$$

(24)

Hence, $D$ can be defined as

$$D=1-{\left[\frac{1}{a}\left(\left(1+a\right)-\sqrt{1+a\left(2+a\right)k}\right)\right]}^{\frac{c}{b}}$$

(25)

By taking derivative from equation (20), the damage evolution equation for the uniaxial state can be written as

$$\dot{k}=\frac{1}{g}f(k){\dot{\epsilon }}^p$$

(26)

In order to convert the damage from a uniaxial damage evolution to a multi-axial damage evolution, the plastic strain rate is calculated from the following equation

$${\dot{\epsilon }}^p={\delta }_{t}r\left(\overline{\sigma }\right){\dot{\epsilon }}_{\max }^p+{\delta }_c\left(1-r\left(\overline{\sigma }\right)\right){\dot{\epsilon }}_{\min }^p$$

(27)

where $\delta$ is the Kronecker delta, ${\dot{\epsilon }}_{\min }^p$ and ${\dot{\epsilon }}_{\max }^p$ are algebraically the minimum and maximum eigen-values of the plastic strain rate tensor, and $r(\overline{\sigma })$ is a weight factor expressed as

$$r\left(\overline{\sigma }\right)=\{\begin{array}{cc}0\mathit{\hspace{1em}}\hfill & \mathrm{if}\overline{\sigma }=0\hfill \\ \left(\sum _{i=1}^3⟨{\overline{\sigma }}_i⟩\right)\mathrm{/}\left(\sum _{i=1}^3\left|{\overline{\sigma }}_i\right|\right)\mathit{\hspace{1em}}\hfill & \mathrm{otherwise}\hfill \end{array}$$

(28)

Substituting equation (24) into (13), the evolution equation can be obtained as

$$\dot{k}=h(\overline{\sigma },k):{\dot{\epsilon }}^p$$

(29)

where

$$h(\overline{\sigma },k)=\left[\begin{array}{ccc}\hfill r\left(\overline{\sigma }\right)f_t\left(k_t\right)\mathrm{/}{g}_t\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill (1-r\left(\overline{\sigma }\right))f_c\left(k_c\right)\mathrm{/}{g}_c\hfill \end{array}\right]$$

(30)

Finally, applying the definition of ${\dot{\epsilon }}^p$ in equation (15), $H(\overline{\sigma },k)$ can be obtained as

$$H=h\cdot {\nabla }_{\overline{\sigma }}\mathrm{\Phi }(\overline{\sigma })$$

(31)

3.2 The Plasticity Model and Hardening Function

The definition of yield criterion within the framework of a plasticity model is significant to describe the behavior of concrete under the tensile and compressive loadings [40, 41]. Implementation of similar behavior of tension and compression in concrete leads to an unrealistic plastic deformation. In this study, the yield criterion originally introduced by Lee and Fenves [15] is proposed, which is defined in the effective stress space using an undamaged configuration parameters as

$$F\left(\sigma \right)=\frac{1}{1-\alpha }\left[\alpha I_1+\sqrt{3J_2}+\beta {\sigma }_{\max }\right]-c_c\left(k\right)$$

(32)

where $I_1$ is the first invariant of stress tensor and $J_2$ is the second invariant of the deviatoric stress tensor, respectively. In above, $\alpha$ and $\beta$ are dimensionless constants, in which $\alpha$ depends on the ratio of yield strength under the biaxial and uniaxial compression, defined as

$$\alpha =\frac{(f_{b0}/f_{c0})-1}{2(f_{b0}/f_{c0})-1}$$

(33)

where $f_{b0}$ is the biaxial and $f_{c0}$ is the uniaxial compressive yield stresses. The experimental values of $f_{b0}/f_{c0}$ is in the range of 1.10 and 1.16 that leads to the value of $\alpha$ between 0.08 and 0.12. The parameter $\beta$ is a constant value that is defined as a dimensionless function of the tensile and compressive cohesion parameters $c_t$ and $c_c$ in the Barcelona model as

$$\beta =\frac{c_c(k)}{c_t(k)}(\alpha -1)-(1+\alpha )$$

(34)

The biaxial tensile strength depends on the parameters $\alpha$ and $\beta$ and is always slightly lower than the uniaxial tensile strength.

The flow rule is defined based on the relation between the plastic flow direction and the plastic strain rate. In this study, a non-associative flow rule is required to control the dilatancy in modeling the frictional behavior of material. Hence, a plastic potential function which is the type of Drucker-Prager yield criteria is utilized as

$$\phi =\sqrt{2J_2}+{\alpha }_p{I}_1$$

(35)

where the parameter ${\alpha }_p$ is chosen such that to obtain a proper value of the dilatancy.

3.3 A Combined Plastic-damage Model for Strong Discontinuity

The damage and plasticity models described in preceding sections are combined within the X-FEM framework to describe the damage in the form of strong discontinuity. To this end, the elastic and plastic strain components are employed within the bulk volume of element while the damage strain is lumped into the crack interface. Hence, the additive decomposition of strain is assumed as

$$\epsilon ={\epsilon }^e+{\epsilon }^p+{\epsilon }^d$$

(36)

where ${\epsilon }^e$ is the elastic strain which is recoverable part of the strain in unloading condition, ${\epsilon }^p$ is the plastic strain which is irreversible part of the strain, and ${\epsilon }^d$ is the other recoverable part of strain which is related to the evolution of elastic modulus during the damage state. The incremental displacement field and strain field of the crack body can be written as

$\dot{\mathbf{u}}(x,t)$	$=\underset{\mathrm{regular}(\mathrm{continuous})}{\underbrace{\dot{\overline{\mathbf{u}}}(x,t)}}+\underset{\mathrm{singular}(\mathrm{discontinuous})}{\underbrace{H_{{\mathrm{\Gamma }}^d}(x)\left[\left[\dot{\mathbf{u}}(x,t)\right]\right]}}$	(37)
$\dot{\epsilon }(x,t)$	$=\underset{\mathrm{regular}(\mathrm{bounded})}{\underbrace{{\nabla }^s\dot{\overline{\mathbf{u}}}(x,t)+H_{\mathrm{\Gamma }}(x){\nabla }^s\left[\left[\dot{\mathbf{u}}(x,t)\right]\right]}}+\underset{\mathrm{unbounded}}{\underbrace{{\delta }_{\mathrm{\Gamma }}(x){\left(\left[\left[\dot{\mathbf{u}}(x,t)\right]\right]\otimes {\mathbf{n}}_{\mathrm{\Gamma }}\right)}^s}}$

where the superscript $s$ denotes the symmetric part of the differential operator ${\nabla }^s$, the Dirac-delta function ${\delta }_{\mathrm{\Gamma }}$ is centered on the discontinuity, and $n_{\mathrm{\Gamma }}$ is the unit normal vector to the discontinuity. A regularized type of equation (38) can be defined employing the regularized version of Dirac-Delta function as

$${\delta }_{\mathrm{\Gamma }}(x)=\frac{1}{\mathrm{\ell }}\mu (x)\mathit{\hspace{1em}}\mathrm{with}\mu (x)=\{\begin{array}{cc}1\mathit{\hspace{1em}}\hfill & x\in \mathrm{\Gamma }\hfill \\ 0\mathit{\hspace{1em}}\hfill & x\in \mathrm{\Omega }\hfill \end{array}$$

(39)

where $\mathrm{\Gamma }$ is the discontinuity and $\mathrm{\Omega }$ is the domain of element. As the damage strain is lumped into the crack interface the following equations can be defined as

${\dot{\mathit{\epsilon }}}^e+{\dot{\mathit{\epsilon }}}^p$	$={\nabla }^s\dot{\overline{\mathbf{u}}}(x,t)+H_{\mathrm{\Gamma }}(x){\nabla }^s\left[\left[\dot{\mathbf{u}}(x,t)\right]\right]$	(40)
${\dot{\mathit{\epsilon }}}^d$	$={\displaystyle \frac{1}{\mathrm{\ell }}}\mu (x){\left(\left[\left[\dot{\mathbf{u}}(x,t)\right]\right]\otimes {\mathbf{n}}_{\mathrm{\Gamma }}\right)}^s$	(41)

where $\mathrm{\ell }$ is the perpendicular distance from the crack interface for each Gauss point of the element, as shown in Figure 1.

Figure 1 An element cut by a crack interface; (a) the crack direction and Gauss points, (b) the crack opening displacements.

The stress – strain constitutive relation can be written in the rate form as

$$\dot{\sigma }={\mathbf{C}}^{ep}({\dot{\epsilon }}^e+{\dot{\epsilon }}^p)$$

(42)

where $\dot{\sigma }$, ${\dot{\epsilon }}^e$ and ${\dot{\epsilon }}^p$ are the stress, elastic strain and plastic strain rates, respectively, and ${\mathbf{C}}^{ep}$ is the elastic-plastic tangent operator [42]. Similarly, the stress – strain constitutive relation for the damage regime can be written as

$$\dot{\sigma }={\mathbf{C}}^d{\dot{\epsilon }}^d$$

(43)

where ${\dot{\epsilon }}^d$ is the damage strain rate and ${\mathbf{C}}^d$ is the damage tangent operator. It is necessary to convert the stress – strain constitutive relation into the traction separation law. Substituting equation (41) into (43) and applying some manipulation, the traction separation relation can be obtained as

$\dot{\mathbf{t}}={\mathbf{n}}_{\mathrm{\Gamma }}\dot{\mathit{\sigma }}={\mathbf{n}}_{\mathrm{\Gamma }}{\mathbf{C}}^d\left({\displaystyle \frac{1}{{\mathrm{\ell }}_i}}{\left(\mathrm{\Delta }\left[\left[\mathbf{u}\right]\right]\otimes {\mathbf{n}}_{\mathrm{\Gamma }}\right)}^s\right)={\displaystyle \frac{1}{{\mathrm{\ell }}_i}}{\mathbf{n}}_{\mathrm{\Gamma }}{\mathbf{C}}^d{\left(\mathrm{\Delta }\left[\left[\mathbf{u}\right]\right]\otimes {\mathbf{n}}_{\mathrm{\Gamma }}\right)}^s$

(44)

or

$\dot{\mathbf{t}}={\displaystyle \frac{1}{{\mathrm{\ell }}_i}}\underset{\mathbf{Q}}{\underbrace{\left({\mathbf{n}}_{\mathrm{\Gamma }}{\mathbf{C}}^d{\mathbf{n}}_{\mathrm{\Gamma }}\right)}}\mathrm{\Delta }\left[\left[\mathbf{u}\right]\right]={\displaystyle \frac{1}{{\mathrm{\ell }}_i}}\mathbf{Q}\mathbf{\Delta }\left[\left[\mathbf{u}\right]\right]$

(45)

where $\dot{\mathbf{t}}$ is the traction rate vector of normal and tangential components, $\mathrm{\Delta }\left[\left[\mathbf{u}\right]\right]$ is the incremental jump of displacement, in which its normal ${\left[\left[\mathbf{u}\right]\right]}^n$ and tangential ${\left[\left[\mathbf{u}\right]\right]}^t$ components are shown in Figure 1, and $\mathbf{Q}$ is the so-called acoustic tensor. Hence, equation (45) can be rewritten as

$$\dot{\mathbf{t}}=\mathbf{T}\mathrm{\Delta }\left[\left[\mathbf{u}\right]\right]$$

(46)

where $\mathbf{T}=\frac{1}{{\mathrm{\ell }}_i}\mathbf{Q}$ represent the tangential modulus matrix of the discontinuity.

4 Cohesive Fracture Model for FRP Interface

The cohesive fracture model was originally introduced by Dugdale [43] and Barenblatt [44], which is based on the assumption that the fracture process zone extends along the crack surfaces. The basic assumption of the cohesive fracture model is that the material is capable of transferring stresses while the fracture process zone occurs. Basically, the micro-cracks happen near to the interface and the macro-cracks occur from the collection of micro-cracks [45]. In the fracture process zone, the cohesive constitutive law is defined on the basis of cohesive fracture traction that is given as a function of the separation of two crack surfaces. The cohesive constitutive law includes the tensile strength of the material $f_t$ and the fracture energy $E_f$. In this model, the effective traction and the effective crack separation are defined as follows

$t^e$	$=\sqrt{{(t^n)}^2+{(t^s)}^2}$	(47)
${\mathrm{\Delta }}^e$	$=\sqrt{{({\mathrm{\Delta }}^n)}^2+{({\mathrm{\Delta }}^s)}^2}$	(48)

where $t^n$ and $t^s$ are the normal and tangential tractions, and ${\mathrm{\Delta }}^n$ and ${\mathrm{\Delta }}^s$ are the normal and sliding displacements of the fracture surface, respectively. The damage initiation occurs when the traction or separation reaches the critical value, i.e., $t^e\ge t^0$ or ${\mathrm{\Delta }}^e\ge {\mathrm{\Delta }}^0$ [31].

The constitutive law can be defined between the cohesive traction $t_j$ and crack separation ${\mathrm{\Delta }}_i$ using $t_j=t({\mathrm{\Delta }}_i)$ or ${\dot{t}}_j=T_{ji}{\dot{\mathrm{\Delta }}}_i$, where $T_{ji}$ is the constitutive tangent stiffness tensor of cohesive fracture [44, 46]. In this study, the delamination model proposed by Turon et al. [47] is utilized based on the continuum damage model. The free energy can be defined per unit area of the crack interface as

$$\phi (\mathrm{\Delta },d)=(1-d){\phi }^0({\mathrm{\Delta }}_i)-d{\phi }^0({\delta }_{1i}⟨-{\mathrm{\Delta }}_1⟩)\mathit{\hspace{1em}}(i=1,2)$$

(49)

where $d$ is a scalar damage variable and ${\phi }^0$ is a function of the displacement separation, defined as ${\phi }^0(\mathrm{\Delta })=\frac{1}{2}\mathrm{\Delta }{}_{i}T_{ij}^0{\mathrm{\Delta }}_j$, with $T^0$ denoting the constitutive undamaged tangent stiffness tensor. In above, $⟨\cdot ⟩$ is the MacAuley bracket defined as $⟨\mathrm{\Xi }⟩=\frac{1}{2}(\mathrm{\Xi }+|\mathrm{\Xi }|)$ and ${\delta }_{ij}$ is the Kronecker delta. The negative value of ${\mathrm{\Delta }}_1$ means that the interface is in contact and the damage cannot occur in the normal direction. By taking the derivative from the free energy (50), the constitutive equation for the crack interface can be obtained as

$$t_i=\frac{\partial \phi }{\partial {\mathrm{\Delta }}_i}=(1-d)T_{ij}^0{\mathrm{\Delta }}_j-d{T}_{ij}^0{\delta }_{1i}⟨-{\mathrm{\Delta }}_1⟩$$

(50)

4.1 The Damage Criterion and Damage Evolution Law

The damage criterion can be defined as a function of the displacement separation as [47]

$$\overline{F}({\mathrm{\Delta }}^e,\overline{d})=G({\mathrm{\Delta }}^e)-\overline{d}\le 0$$

(51)

where $\overline{d}$ is the damage threshold and $G$ is a function of damage evolution which is defined in the range of $[0,1]$ as

$$G({\mathrm{\Delta }}^e)=\frac{{\mathrm{\Delta }}^f({\mathrm{\Delta }}^e-{\mathrm{\Delta }}^0)}{{\mathrm{\Delta }}^e({\mathrm{\Delta }}^f-{\mathrm{\Delta }}^0)}$$

(52)

where ${\mathrm{\Delta }}^0$ is the onset displacement jump, ${\mathrm{\Delta }}^e$ is the effective separation and ${\mathrm{\Delta }}^f$ is the final displacement jump that can be determined using the fracture energy $G_f$. Depending on the definition of damage, the damage evolution can be given by $G({\mathrm{\Delta }}^e)=d=1-T/T^0$. If the linear cohesive constitutive relation is considered, $T$ and $T^0$ can be defined as

$T^0$	$={\displaystyle \frac{t^0}{{\mathrm{\Delta }}^0}}$	(53)
$T$	$={\displaystyle \frac{t^e}{{\mathrm{\Delta }}^e}}={\displaystyle \frac{t^0\left({\mathrm{\Delta }}^f-{\mathrm{\Delta }}^e\right)}{{\mathrm{\Delta }}^e({\mathrm{\Delta }}^f-{\mathrm{\Delta }}^0)}}$	(54)

where $t^0$ is the strength of material, as shown in Figure 2. The damage initiates when the effective separation $G({\mathrm{\Delta }}^e)$ exceeds the initial damage threshold. It is required to define an evolution law for the damage model that can be defined as

$$\dot{d}=\dot{\mu }\frac{\partial \overline{F}({\mathrm{\Delta }}^e,\overline{d})}{\partial {\mathrm{\Delta }}^e}=\dot{\mu }\frac{\partial G({\mathrm{\Delta }}^e)}{\partial {\mathrm{\Delta }}^e}$$

(55)

where $\dot{\mu }$ is defined as the damage consistency parameter that should satisfy the loading–unloading conditions, i.e.

$$\dot{\mu }\ge 0,\overline{F}({\mathrm{\Delta }}^e,\overline{d})\le 0,\dot{\mu }F({\mathrm{\Delta }}^e,\overline{d})=0$$

(56)

Figure 2 The bilinear constitutive law for the FRP interface.

In order to derive the constitutive tangent stiffness tensor, the damage model is implemented through a nonlinear approach. Taking the derivative from equation (50) leads to

${\dot{t}}_i={\delta }_{ij}{T}_{ij}^0\left[1-d\left(1+{\delta }_{1j}{\displaystyle \frac{⟨-{\mathrm{\Delta }}_j⟩}{{\mathrm{\Delta }}_j}}\right)\right]{\dot{\mathrm{\Delta }}}_j-{\delta }_{ij}{T}_{ij}^0\left[1+{\delta }_{1j}{\displaystyle \frac{⟨-{\mathrm{\Delta }}_j⟩}{{\mathrm{\Delta }}_j}}\right]{\mathrm{\Delta }}_j\dot{d}$ (57)

and the evolution of the damage variable $d$ can be obtained as

$$\dot{d}=\{\begin{array}{cc}\hfill \dot{G}({\mathrm{\Delta }}^e)=\frac{\partial G({\mathrm{\Delta }}^e)}{\partial {\mathrm{\Delta }}^e}{\dot{\mathrm{\Delta }}}^e\mathit{\hspace{1em}}\hfill & \mathrm{for}\overline{d}<{\mathrm{\Delta }}^e<{\mathrm{\Delta }}^f\hfill \\ \hfill =\frac{{\mathrm{\Delta }}^f{\mathrm{\Delta }}^0}{{\mathrm{\Delta }}^f-{\mathrm{\Delta }}^0}\frac{1}{{({\mathrm{\Delta }}^e)}^2}\mathit{\hspace{1em}}\hfill & \\ \hfill 0\mathit{\hspace{1em}}\hfill & \mathrm{for}\overline{d}>{\mathrm{\Delta }}^e\mathrm{or}{\mathrm{\Delta }}^e>{\mathrm{\Delta }}^f\hfill \end{array}$$

(58)

Finally, the constitutive tangent stiffness tensor $T_{ji}^{\tan }$ in the constitutive law ${\dot{t}}_j=T_{ji}^{\tan }{\dot{\mathrm{\Delta }}}_i$ can be defined as

$$T_{ji}^{\tan }=\{\begin{array}{cc}{\delta }_{ij}{T}_{ij}^0[1-d{H}_j]\mathit{\hspace{1em}}\hfill & \mathrm{for}\overline{d}<{\mathrm{\Delta }}^e<{\mathrm{\Delta }}^f\hfill \\ -T_{ij}^0{H}_j{H}_i\frac{{\mathrm{\Delta }}^f{\mathrm{\Delta }}^0}{{\mathrm{\Delta }}^f-{\mathrm{\Delta }}^0}\frac{1}{{({\mathrm{\Delta }}^e)}^2}\mathit{\hspace{1em}}\hfill & \\ {\delta }_{ij}{T}_{ij}^0\left[1-d{H}_j\right]\mathit{\hspace{1em}}\hfill & \mathrm{for}\overline{d}>{\mathrm{\Delta }}^e\mathrm{or}{\mathrm{\Delta }}^e>{\mathrm{\Delta }}^f\hfill \end{array}$$

(59)

where $H_j=1+{\delta }_{1j}⟨-{\mathrm{\Delta }}_j⟩/{\mathrm{\Delta }}_j$.

5 Crack Initiation and Crack Growth Direction Criteria

In crack growth simulation, there are basically two significant issues that must be taken into account; the crack initiation and the crack propagation direction criteria. In linear elastic fracture mechanics, there are different criteria proposed by researchers for the mixed mode fracture on the basis of the critical stress intensity factor and fracture toughness; these criteria are the minimum strain energy density criterion, the maximum hoop stress criterion, and the maximum energy release rate [31]. However, these criteria are not able to accurately predict the behavior of local non-linearity at the crack-tip region described by the fracture process zone as these criteria are defined based on the whole condition of the structure. In the current study, an approach is proposed based on the damage state that captures the local non-linearity at the fracture process zone. In this way, the crack propagates if the crack segment length ${\mathrm{\ell }}_c$ is damaged. The damage is computed at each Damage Calculation Point (DCP), which is defined at a region of the distant ${\mathrm{\ell }}_c$ from the crack-tip at different angles of ${\theta }_{CP}$ in the range of $-90\le {\theta }_{CP}\le 90$. The value of damage is computed at each DCP by applying a weighted average value of damage at all neighboring Gauss points defined as

$$D_{CP}=\sum _j\left(\frac{D_j}{r_j^2}\right)\mathrm{/}\sum _j\left(\frac{1}{r_j^2}\right)$$

(60)

where $D_{CP}$ is the damage at each DCP and, $D_j$ and $r_j$ are respectively the damage and the distance from each DCP at a nearby Gauss point. The damage value can be obtained by determination of the element containing the DCP, and implementation of the FEM shape functions to damage nodal values. The damage value $D_{CP}$ must be calculated for all DCPs at each loading step, and the computed values must be compared with the critical damage value $D_{cr}$; if $D_{CP}\ge D_{cr}$ the crack propagates in the direction of ${\theta }_{CP}$ at the DCP.

6 Validation of Computational Algorithm

6.1 X-FEM Modeling of a Three-point Bending Beam with an Edge Crack

In order to illustrate the performance of the proposed X-FEM technique, a simply supported beam is modeled, as shown in Figure 3. The beam is subjected to a prescribed displacement at its center on the top edge. The material properties of the beam are chosen as follows; the Young modulus of $E=100\mathrm{MPa}$, the Poisson ratio of $\nu =0.35$, the tensile strength of material $f_t=1.0\mathrm{MPa}$ and the fracture energy $E_f=0.1\mathrm{N}/\mathrm{mm}$. The elements of top edge of the beam are prevented from cracking since the system of equations becomes singular if a crack propagates through the entire beam. The crack commences from the center of the beam at the bottom edge and propagates toward the upper edge. The beam is modeled applying a coarse and a fine mesh with $860$ and $1530$ bilinear elements, respectively. The coarse FE mesh is shown in Figure 3. In Figure 4, the load-displacement curves are plotted for the coarse and fine meshes, and the results are compared with that reported by Wells and Sluys [48]. In Figure 5, the contours of stress in $x-$direction are presented at various displacements where the crack-tip approaches the upper edge of the beam. Finally, the variation of predicted cohesive traction along the crack interface is plotted in Figure 6 with the crack mouth opening displacement (CMOD).

Figure 3 A three-point bending beam with symmetric edge crack; Geometry, boundary conditions and X-FEM mesh.

Figure 4 The load-displacement curves for a three-point bending beam.

Figure 5 The distribution of stress contours in $x-$direction at different steps of crack growth.

Figure 6 The cohesive traction at the crack mouth versus crack mouth opening displacement.

6.2 Modeling the Concrete Behavior with a Plastic-damage Model

In order to verify and validate the performance of the proposed plastic-damage model for the concrete behavior, a single element of $10\times 10{\mathrm{cm}}^2$ is modeled under the monotonic uniaxial compressive and monotonic uniaxial tensile loadings as well as the cyclic uniaxial tensile loading. The material parameters for the numerical modeling of the concrete element are obtained from the experimental test conducted by Gopalaratnam and Shah [49]. In Figure 7(a–b), the strain-stress curves are plotted for the monotonic uniaxial compression and uniaxial tension loadings. A complete agreement can be seen between the predicted results and the experimental stress-strain curves proposed for the compression and tension of the concrete. Also plotted in Figure 7(c) is the strain-stress curve of the tensile cyclic loading; it shows that the damage model can properly describe the behavior of concrete in tensile cyclic loading. It can be seen from this figure that the strain-stress curves obtained from the numerical analysis is in excellent agreement with that of the experimental test. This example clearly demonstrate that the proposed plastic-damage model can be efficiently used to capture the behavior of concrete.

Figure 7 The stress-strain curves of a concrete element with the damage-plasticity model; (a) the monotonic uniaxial compression test, (b) the monotonic uniaxial tension, (c) the cyclic uniaxial tensile test.

6.3 Cohesive Fracture Modeling of the Debonding Test

In order verify the cohesive interface model proposed between the FRP and concrete, the results of the debonding test reported by Au and Büyüköztürk [50] is compared with those obtained using the cohesive interface model. The cohesive model is employed along the interface between the FRP and concrete through the FE analysis. The material properties of the concrete and FRP as well as the inteface between the FRP and concrete are obtained from the experimental test. In Figure 8, the setup of experimental test is presented together with the interface fracture characterization of debonding in FRP plated concrete. In Figure 9, a comparison of the force-displacement curves is performed between the experimental and numerical results. It clearly shows that the cohesive interface model can properly capture the interface behavior between the FRP and concrete. The difference between two diagrams in the softening part of the force-displacement curve can be attributed to the lack of uniformity of the adhesive thickness and the presence of air bubbles.

Figure 8 The debonding test; (a) the setup of experimental test conducted by Au and Büyüköztürk (2006), (b) the interface fracture of FRP strip together with the specimen configuration.

Figure 9 The debonding test; A comparison of the force-displacement curves between the experimental and numerical results.

6.4 Numerical and Experimental Investigations of Concrete Arch

In order to illustrate the performance of the proposed computational algorithm, the two desired concrete arches strengthened with FRP are analyzed numerically and the results are compared with the experimental tests; the first experiment was conducted by Zhang et al. [7] and the second experiment was conducted by the authors at the Strong Floor Laboratory of Sharif University of Technology. In order to present the accuracy of the proposed computational algorithm, the force-displacement curves obtained from the finite element analysis are compared with those reported through the experimental tests. The reinforced concrete arch is modeled using the four-noded bilinear element and the rebar is modeled using a two-noded beam element. The rebars are connected to the concrete using the embedded technique inside the element.

Figure 10 The concrete arch strengthened with internal FRP conducted by Zhang et al. (2015); (a) the concrete arch specimen, (b) the geometry, boundary conditions and the FE mesh of the concrete arch.

In the first case, the concrete arch strengthened with internal FRP, originally conducted by Zhang et al. [7], is modeled using the proposed computational algorithm, as shown in Figure 10. In this figure, the concrete arch specimen is shown together with the geometry, boundary conditions and the FE mesh of the concrete arch. It was observed from the experimental investigation that the cracks are first initiated in the middle of the span and then propagated by increasing the load through the shoulders. It was also observed that a five-hinged arch is formed where the structure becomes unstable. It is interesting to highlight that the five-hinged structure was also detected through the numerical analysis of the concrete arch. In Figure 11, a comparison of the force-displacement curves is performed between the experimental and numerical results; it clearly demonstrates that the proposed computational algorithm can properly capture the behavior of the concrete arch. Obviously, a slight difference can be seen between the numerical analysis and experimental data at the beginning of the loading, which is due to the local splitting of the sample at the time of the sitting of the sample, which cannot be modeled in the numerical model. Furthermore, the numerical analysis demonstrates the higher strength compared to the experiment which can be described due to uncertainties in the properties of the materials utilized in the laboratory and also because of the reduced concrete stress in high strains. According to the results, it can be concluded that the numerical analysis adequately shows a good agreement with those of experiment.

Figure 11 The concrete arch strengthened with internal FRP conducted by Zhang et al. (2015); A comparison of the force-displacement curves between the experimental and numerical results.

Figure 12 The concrete arch strengthened with the internal & external FRP conducted at the Strong Floor Laboratory of Sharif University of Technology; (a) the setup of experimental test, (b) the problem definition, (c) the geometry, boundary conditions and FE mesh of concrete arch.

Figure 13 The concrete arch strengthened with internal & external FRP conducted at the Strong Floor Laboratory of Sharif University of Technology; (a) the crack trajectories obtained from experimental test, (b) the contour of plastic strain obtained from the numerical analysis.

In the second case, the concrete arch strengthened with the internal and external FRP experimentally conducted by the authors at the Strong Floor Laboratory of Sharif University of Technology is modeled, as shown in Figure 12. The data of one test is employed to verify the accuracy of the proposed computational model. The concrete arch specimen has a span length of 110 cm and the high of 35 cm with a rectangular cross-section of $30\times 10{\mathrm{cm}}^2$. The internal and external parts of the arch are retrofitted by a FRP layer with the thickness of 0.35 mm. Three longitudinal bars with a diameter of 8 mm and stirrups with a diameter of 6 mm, with the spacing of 100 mm are used (Figure 12c). All rebars have the strength of 400 MPa. The compressive strength of the concrete varies between 47 and 54 MPa, in which a compressive strength of 48 MPa is considered in the numerical analysis. The parameter of cohesive fracture in the FRP interface is chosen based on the observation of concrete crumbs over the FRP interface after performing the test, so the fracture energy is chosen equal to 146 N/m. It can be observed that the failure of the arch happens in a shear failure mode due to the excessive distance between the stirrups, as shown in Figure 13. The maximum distance of the stirrups should be half of the effective high of the section. Given the section high of 100 mm, the stirrup spacing should be limited to a maximum of 50 mm, however, the distance of stirrups is selected equal to 100 mm in the current experimental test due to practical issues. Because of the tensile stress occurs along the interface between the FRP and concrete at the interface of the arch, debonding happens at this region, as shown in Figure 13. A comparison between Figures 13(a) and 13(b) illustrates that the proposed computational algorithm can be used accurately to capture the crack trajectories as well as debonding of the FRP. In Figure 14, the force-deflection curves are plotted for the experiment and numerical model, which are in a good agreement.

Figure 14 The concrete arch strengthened with internal & external FRP conducted at the Strong Floor Laboratory of Sharif University of Technology; A comparison of the force-displacement curves between the experimental and numerical results.

7 Parametric Studies of FRP Strengthening of Concrete Arches

In this section, the behavior of FRP strengthening of concrete arches is studied numerically for the concrete arch of the span length 110 cm at different highs of $H=10$, $30$, $45$ and $55\mathrm{cm}$ (see Figure 12b) using two percentages of reinforcement of 0.628% and 2.5%, which are retrofitted in three manners, i.e. the Internal FRP (I), External FRP (E), and both External and Internal FRP (E&I). The behavior of the retrofitted arches are compared with the arches without the FRP strengthening. In addition, an additional model is considered in which the debonding of the internal FRP sheet is controlled. A rectangular cross-section of $30\times 10{\mathrm{cm}}^2$is proposed for all concrete arches with 2 cm concrete cover over the rebars. According to the ACI 318-14, the minimum and maximum area of rebars are calculated as follows

${(A_s)}_{\min }$	$=\max ({\displaystyle \frac{0.25\sqrt{f_c^{\prime }}}{f_y}}{b}_{w}d,{\displaystyle \frac{1.4}{f_y}}{b}_{w}d)$
${(A_s)}_{\max }$	$=0.025b_{w}d$	(61)

where ${(A_s)}_{\min }$ and ${(A_s)}_{\max }$ are respectively the minimum and maximum area of the flexural reinforcements, $f_c^{\prime }$ is the specified compressive strength of concrete, $f_y$ is the specified yield strength of the reinforcement, $b_w$ is the web-width of the section, and $d$ is the distance from the extreme compression fiber to the centroid of the longitudinal tensile reinforcements. Accordingly, the minimum allowable rebar percentage is equal to $104{\mathrm{mm}}^2$ and the maximum is $600{\mathrm{mm}}^2$.

7.1 The Concrete Arch with the High of $\mathit{H}$ $=$ 10 cm

In the case of the concrete arch with the high of $H=10\mathrm{cm}$, the behavior of the arch is almost similar to a straight beam due to the low ratio of high to span of the arch. In such case, the crack initiates from the bottom center and propagates to the top center. In Figure 15, the force–deflection curves of the concrete arch are plotted for various retrofitting methods, including: No FRP (N), Internal FRP (I), External FRP (E), and both External and Internal FRP (E&I) at two percentages of reinforcement, i.e. $\rho =0.628\%$ and $2.5\%$. Obviously, the ultimate strength and maximum deformation of the arch increase by increasing the percentages of reinforcement. Accordingly, the ultimate strength of the arch without the FRP increases from $57\mathrm{kN}$ to $110\mathrm{kN}$ with the corresponding deflections from 26 mm to 40 mm when the reinforcements increase from 0.628% to 2.5%; this indicates that the ultimate strength is improved about 100% and the deformation capacity increases about 50%. As mentioned earlier, the crack initiates from the bottom center and propagates to the top center. Before the formation of crack, the stiffness of the arch is identical for the two percentages of reinforcement, however, the stiffness of the arch with higher reinforcement value is more than the arch with lower longitudinal reinforcement when the crack happens. The post cracking stiffness of the arch with $\rho =0.628\%$ and $2.5\%$ are respectively 15% and 32% of the initial stiffness that indicates the effect of longitudinal reinforcement on the post-cracking stiffness.

Figure 15 The force–deflection curves of concrete arch with $H=10\mathrm{cm}$ at different percentages of reinforcement; (a) No FRP (N), (b) Internal FRP (I), (c) External FRP (E), (d) External and Internal FRP (E&I), (e) The effect of FRP for $\rho =0.628\%$, (f) The effect of FRP for $\rho =2.5\%$.

Figure 16 The plastic strain contours of concrete arch with $H=10\mathrm{cm}$ at ultimate strength; (a) No FRP (N) with $\rho =0.628\%$ , (b) No FRP (N) with $\rho =2.5\%$, (c) Internal FRP (I) with $\rho =0.628\%$, (d) Internal FRP (I) with $\rho =2.5\%$, (e) External FRP (E) with $\rho =0.628\%$, (f) External FRP (E) with $\rho =2.5\%$, (g) External and Internal FRP (E&I) with $\rho =0.628\%$, (h) External and Internal FRP (E&I) with $\rho =2.5\%$.

One of the most common failure modes in FRP-strengthened concrete elements is the debonding of the FRP from the surface of concrete. The force-deflection curves plotted in Figure 15(b) clearly present the deboning of the FRP from the concrete with Internal FRP. Moreover, it can be seen that the ultimate strength of the arch is $137\mathrm{kN}$ for $\rho =2.5\%$ with the corresponding deflection of 51 mm, while the ultimate strength is $80\mathrm{kN}$ for $\rho =0.628\%$ with the deflection of 20 mm. It can be observed that increasing the reinforcement from $\rho =0.628\%$ to $2.5\%$ leads to 70% increase in the ultimate strength value. According to Figure 15(b), a sudden drop is obvious, which is due to the debonding of the FRP from the concrete. The main reason of the debonding of the FRP is the presence of tensile stress on the adhesive layer. The higher percentage of reinforcement causes the debonding to occur in higher value of loading. In Figure 16, the plastic strain contours of concrete arch are presented for various retrofitting methods at two percentages of reinforcement. It can be seen that the crack happens at the bottom face and then propagates to the top face; hence, the major failure mode is due to flexure. Similar to the arch without FRP, the reinforcement affects the stiffness of the arch only in the post-cracking stage. Moreover, reduction of the stiffness of the arch with lower percentage of reinforcement is more pronounced. The post-cracking stiffness of the arch with $\rho =0.628\%$ and $2.5\%$ are respectively 16% and 33% of the initial stiffness value. Clearly, the stiffness of the arch decreases significantly when the debonding occurs, however, a fraction of the stiffness is recovered when the load increases. This is due to the fact that increasing the load causes the tensile strength in the FRP after debonding. The strengthening of exterior surface of the arch slightly improves the behavior of concrete arch. In such case, the arch experiences only compression stress on the exterior surface; as a result, it has little effect on the behavior of the arch. According to Figures 15(e) and 15(f), the force-deflection curves of the arches without strengthening are similar to the arches with external strengthening. Furthermore, debonding is more significant in the arches with lower percentage of reinforcement. In addition, the higher reinforcement value leads to higher deformation capacity. It is obvious that the internal FRP strengthening improves the maximum strength of the arch up to 65%. The stiffness reduction can be observed in the post-cracking stage for the arches with Internal FRP strengthening. Figures 16(a–e) and 16(c–g) represent the effect of Internal FRP strengthening on the crack trajectory. Obviously, debonding can be observed in the arch with Internal FRP at the ultimate strength. It can be seen that the arch with both Internal and External FRP demonstrates the higher strength compared to the arch with Internal FRP; this is due to the fact that the External FRP strengthening causes the concrete experiences less tensile stress.

7.2 The Concrete Arch with the High of $\mathit{H}$ $=$ 30 cm

In the case of concrete arch with the high of $H=30\mathrm{cm}$, where the high to span ratio increases, the arch experiences the tensile stress at the top shoulders, hence, the Internal and External FRP strengthening method can be used as an effective scheme. In Figure 17, the force–deflection curves are plotted for the concrete arch with $H=30\mathrm{cm}$ for various retrofitting methods, i.e. No FRP (N), Internal FRP (I), External FRP (E), and both External and Internal FRP (E&I) at two reinforcement values of $\rho =0.628\%$ and $2.5\%$. It is obvious from Figure 17(a) that by increasing the reinforcement ratio in the arch without FRP, the ultimate strength increases about 25%. Moreover, it can be seen from Figure 17(b) that the arch with the reinforcement ratio of 2.5% demonstrates a significant stiffness degradation due to shear cracking. Also, a drop can be observed in both force-deflection curves which is due to deboning of FRP in the arch with Internal FRP (I). However, a partial recovery of the stiffness of the arch can be observed until the stirrups fail at the ultimate stress. In Figure 18, the plastic strain contours of concrete arch are presented for various retrofitting methods with two percentages of reinforcement at the ultimate load. Obviously, the shear cracking occurs from the bottom and propagates to the top of the arch; this is the major difference between the current case and that studied in Section 7.1. The stiffness of the arches is almost similar before occurring the crack, however it shows different behavior in the post-cracking stage, in which the arches with the reinforcement ratio of 2.5% have the greater stiffness value.

Figure 17 The force–deflection curves of concrete arch with $H=30\mathrm{cm}$ at different percentages of reinforcement; (a) No FRP (N), (b) Internal FRP (I), (c) External FRP (E), (d) External and Internal FRP (E&I), (e) The effect of FRP for $\rho =0.628\%$, f) The effect of FRP for $\rho =2.5\%$.

Figure 18 The plastic strain contours of concrete arch with $H=30\mathrm{cm}$ at ultimate strength; (a) No FRP (N) with $\rho =0.628\%$, (b) No FRP (N) with $\rho =2.5\%$, (c) Internal FRP (I) with $\rho =0.628\%$, (d) Internal FRP (I) with $\rho =2.5\%$, (e) External FRP (E) with $\rho =0.628\%$, (f) External FRP (E) with $\rho =2.5\%$, (g) External and Internal FRP (E&I) with $\rho =0.628\%$, (h) External and Internal FRP (E&I) with $\rho =2.5\%$.

According to Figure 17(c), the ultimate strength of the arch with External FRP (E) using $\rho =0.628\%$ is greater than that of the arch with $\rho =2.5\%$; this is due to the exceedance of the reinforcement ratio from its maximum allowable value. As a result, the flexural cracks first occur and then the shear cracks form and propagate from the middle bottom of the arch at $\rho =0.628\%$, as shown in Figure 18(e). However, only the shear cracks initiate and propagate in the arch with the high reinforcement ratio of $\rho =2.5\%$. It can be observed from Figures 18(g) and 18(h) that the main reason of the failure in the arch with External and Internal FRP (E&I) is due to shear cracking in the middle of the arch. According to Figure 17(d), the ultimate strength of the arch with $\rho =0.628\%$ is about 240 kN at the deflection of 24 mm, while its value for the arch with $\rho =2.5\%$ is 190 kN at the deflection of 19 mm; the main reason is because of the occurrence of shear cracks before flexural cracks in the arch.

Finally, it can be highlighted from Figures 17(e) and 17(f) that the highest strength value can be obtained for the arch with External FRP (E) using $\rho =0.628\%$. However, if the debonding of the internal FRP is avoided, the arch with External and Internal FRP (E&I) represents a greater strength value. Moreover, the initial stiffness of all concrete arches are almost similar using $\rho =2.5\%$, where the stiffness reduction can be observed in the post-cracking stage. Moreover, it can be seen from Figure 17(f) that the high strength value belongs to the arch with No FRP (N) and the arch with External FRP (E) due to their brittle behavior. It must be noted that the method of retrofitting has negligible effect on the post-cracking stiffness due to the high reinforcement ratio. Since shear cracking is the main reason of the failure of the arch, the exceedance of longitudinal reinforcement from its maximum allowable value results in a brittle behavior. In order to investigate the effect of deboning of the FRP in the arch, a comparison is performed in Figure 19 between the arch with Internal FRP (I) where the deboning can be happened, the arch with Internal FRP (I) where the deboning is restricted, and the arch with External FRP (E). Obviously, the ultimate strength of the arch with External FRP (E) is almost equal to the arch with Internal FRP (I) where the deboning is restricted.

Figure 19 The force–deflection curves of concrete arch with $H=30\mathrm{cm}$ for the reinforcement of $\rho =0.628\%$; A comparison between the Internal FRP (I) with debonding, Internal FRP (I) without deboning, and External FRP (E).

Figure 20 The force–deflection curves of concrete arch with $H=45\mathrm{cm}$ at different percentages of reinforcement; (a) No FRP (N), (b) Internal FRP (I), (c) External FRP (E), (d) External and Internal FRP (E&I), (e) The effect of FRP for $\rho =0.628\%$, (f) The effect of FRP for $\rho =2.5\%$.

7.3 The Concrete Arch with the High of $\mathit{H}$ $=$ 45 cm

In the next case, the concrete arch is studied for the high of $H=45\mathrm{cm}$. The force–deflection curves are plotted in Figure 20 for various retrofitting methods, i.e. No FRP (N), Internal FRP (I), External FRP (E), and both External and Internal FRP (E&I) at two reinforcement values of $\rho =0.628\%$ and $2.5\%$. In Figure 21, the plastic strain contours are also presented for the concrete arch of various retrofitting methods at two reinforcement ratios. In the arch with No FRP (N), the stiffness of the arch with the reinforcement ratio of 2.5% is slightly greater than that of the arch with $\rho =0.628\%$, as shown in Figure 20(a); moreover, the ultimate strength of the arch is 190 kN at the deflection of 9 mm for $\rho =0.628\%$, and the ultimate strength is 198 kN at the deflection of 13 mm for $\rho =2.5\%$. However, two different mechanisms of the failure can be seen in the two concrete arches, as shown in Figures 21(a) and 21(b). In fact, the arch with $\rho =0.628\%$ represents the occurrence of cracks at the middle of the arch and the formation of plastic hinge in the shoulder of the arch; while the arch with $\rho =2.5\%$ demonstrates the shear cracking in the arch.

Figure 21 The plastic strain contours of concrete arch with $H=45\mathrm{cm}$ at ultimate strength; (a) No FRP (N) with $\rho =0.628\%$ , (b) No FRP (N) with $\rho =2.5\%$, (c) Internal FRP (I) with $\rho =0.628\%$, (d) Internal FRP (I) with $\rho =2.5\%$, (e) External FRP (E) with $\rho =0.628\%$, (f) External FRP (E) with $\rho =2.5\%$, (g) External and Internal FRP (E&I) with $\rho =0.628\%$, (h) External and Internal FRP (E&I) with $\rho =2.5\%$.

In Figure 20(b), the force–deflection curves are plotted for the arch with Internal FRP. Obviously, because of the tensile stress at the internal surface, debonding occurs between the FRP and concrete in the arch with $\rho =0.628\%$ that causes a sudden drop in the force–deflection curve. The main failure of this arch is its instability due to formation of the hinges at the shoulders of the arch, resulting in a five-hinged arch. It can be observed that the ultimate strength of the arch with $\rho =0.628\%$ is 200 kN at the deflection of 8 mm, while the ultimate strength for $\rho =2.5\%$ is 190 kN at the deflection of 12 mm. Moreover, it can be observed that the ultimate strength of the arch increases by applying the External FRP for the arch with $\rho =0.628\%$, as shown in Figures 20(c) and 20(d), in which the ultimate strength of the arch with External FRP is 275 kN at the deflection of 20 mm and the ultimate strength is 280 kN at the deflection of 23 mm for the External and Internal FRP. Clearly, cracks occur at the middle and shoulders of the arch that leads to the failure of a five-hinged arch, as shown in Figures 21(e) and 21(g). The force-deflection curves plotted in Figure 20(e) clearly illustrate that the External FRP is an efficient manner for the arch with the high of $H=45\mathrm{cm}$ at the lower percentage of reinforcement. It shows that the external strengthening improves the ultimate strength and deflection of the arch about 50%. Moreover, it can be seen from Figures 21(a), 21(c), 21(e) and 21(g) that the failure of the arch with $\rho =0.628\%$ happens when a five-hinged arch occurs, while the failure of the arch with $\rho =2.5\%$ happens when the shearing cracks take place, as shown in Figures 21(b), 21(d), 21(f) and 21(h).

Figure 22 The force–deflection curves of concrete arch with $H=55\mathrm{cm}$ at different percentages of reinforcement; (a) No FRP (N), (b) Internal FRP (I), (c) External FRP (E), (d) External and Internal FRP (E&I), (e) The effect of FRP for $\rho =0.628\%$, f) The effect of FRP for $\rho =2.5\%$.

Figure 23 The plastic strain contours of concrete arch with $H=55\mathrm{cm}$ at ultimate strength; (a) No FRP (N) with $\rho =0.628\%$ , (b) No FRP (N) with $\rho =2.5\%$, (c) Internal FRP (I) with $\rho =0.628\%$, (d) Internal FRP (I) with $\rho =2.5\%$, (e) External FRP (E) with $\rho =0.628\%$, (f) External FRP (E) with $\rho =2.5\%$, (g) External and Internal FRP (E&I) with $\rho =0.628\%$, (h) External and Internal FRP (E&I) with $\rho =2.5\%$.

7.4 The Concrete Arch with the High of $\mathit{H}$ $=$ 55 cm

The last case study refers to a semi-circular concrete arch with the high of $H=55\mathrm{cm}$ that demonstrates the behavior of the half-circle arch. In Figure 22, the force–deflection curves are plotted for various retrofitting methods using two reinforcement values. Also presented in Figure 23 are the corresponding plastic strain contours at the ultimate strength of the arch. It can be seen from the force-deflection curves that the arch with No FRP using $\rho =2.5\%$ has the ultimate strength of 190 kN at the deflection of 13 mm, while the ultimate strength of the arch with $\rho =0.628\%$ is 140 kN at the deflection of 8 mm. The force–deflection curves of Figure 22(b) demonstrate that debonding occurs between the FRP and concrete for the arch with Internal FRP. It can be observed from Figures 22(c) and 22(d) that the External FRP increases the ultimate strength of the arch using $\rho =0.628\%$. Obviously, the ultimate strength of the arch with External FRP is 270 kN at the deflection of 13 mm and its value for the arch with External and Internal FRP is 245 kN at the deflection of 15 mm. It can be seen from Figure 23 that the main failure of the arch with External FRP is due to shearing cracks as the External FRP prevents the formation of flexural cracks at the shoulders of the arch. In such case, shearing cracks occur at the middle of the arch where debonding happens between the FRP and concrete. It can be highlighted that the arch with $\rho =2.5\%$ illustrates more brittle behavior in comparison with the arch with $\rho =0.628\%$, so shearing cracks happen at lower loading level. However, the failure of the arch with No FRP and the arch with Internal FRP is due to the formation of flexural cracks at the middle of the arch.

Figure 24 The force–deflection curves of concrete arch with different high; (a) No FRP (N) with $\rho =0.628\%$, (b) No FRP (N) with $\rho =2.5\%$, (c) Internal FRP (I) with $\rho =0.628\%$, (d) Internal FRP (I) with $\rho =2.5\%$, (e) External FRP (E) with $\rho =0.628\%$, (f) External FRP (E) with $\rho =2.5\%$, (g) External and Internal FRP (E&I) with $\rho =0.628\%$, (h) External and Internal FRP (E&I) with $\rho =2.5\%$.

Figure 25 The stress evolutions along the FRP interface for a half of the concrete arch with different high before debonding; (a) the normal stress along the Internal FRP (I) interface, (b) the normal stress along the External FRP (E) interface, (c) the shear stress along the Internal FRP (I) interface, (d) the shear stress along the External FRP (E) interface.

A comparison of force–deflection curves is performed in Figure 24 between the preceding case studies for the concrete arches with the high of $H=10$, $30$, $45$ and $55\mathrm{cm}$. Clearly, increasing the high of the arch increases the stiffness of the concrete arch for both values of reinforcement. However, the ultimate strength of the arch with the high of $55\mathrm{cm}$ is less than the arches of $30$ and $45\mathrm{cm}$ using the reinforcement of $\rho =0.628\%$. It is because of the geometry of the semi-circular arch that leads to the formation and propagation of flexural cracks at the shoulders of the arch in a lower loading level. In addition, the arch of $30\mathrm{cm}$ with Internal FRP has the highest value of ultimate strength. Also, the ultimate strength of the arch with External FRP increases as the high of the arch increases using the reinforcement of $\rho =0.628\%$, in which the arches of $45$ and $55\mathrm{cm}$ have the highest ultimate strength. Moreover, the failure of the arches with the high of $10$ and $30\mathrm{cm}$ is due to formation of shearing cracks in the middle of the arch, while the failure of the arches with the high of $45$ and $55\mathrm{cm}$ is due to formation of cracks in the middle and shoulders of the arch that leads to its instability and failure where a five-hinged arch occurs. It can also be observed from the force-deflection curves of the arch with $\rho =2.5\%$ that the ultimate strengths of the arches with the high of $45$ and $55\mathrm{cm}$ are almost similar, while the arch of $H=30\mathrm{cm}$ represents the higher value of the ultimate strength. Furthermore, the failure of the arches with $\rho =2.5\%$ is due to shearing cracks because of the high value of reinforcement, except for the arch of $10\mathrm{cm}$ that demonstrates a beam-type behavior. Finally, the evolutions of normal and shear stresses are plotted in Figure 25 for a half of the arch with Internal FRP and External FRP interfaces at different high before debonding. Obviously, the arch of External FRP with low ratio of the high to span experiences the high stress values at its shoulders, while the arch of Internal FRP with high ratio of the high to span experiences the high values of stresses.

8 Conclusion

In the present paper, an enriched–FEM method was presented based on the X-FEM technique by applying a damage–plasticity model to investigate the effect of FRP strengthening on concrete arches, and employing a cohesive fracture model to capture the interface behavior between the FRP and concrete. The validity of the proposed X-FEM method was presented using a three-point bending beam with an edge crack. In addition, the accuracy of the plastic-damage model was studied under the monotonic tension, compression, and cyclic tension loading, and the accuracy of the cohesive fracture model was investigated using the experimental data reported for the debonding test. In order to verify the accuracy of the proposed computational algorithm, the numerical simulation results of FRP strengthening of concrete arches were compared with two sets of experimental test; the first experimental test was reported in the literature and the second experiment was conducted by the authors at the Strong Floor Laboratory of Sharif University of Technology. Finally, a parametric study was performed by evaluating the effects of high to span ratio, longitudinal reinforcement ratio, and strengthening method. To this end, the behavior of FRP strengthening of concrete arches was studied numerically for the concrete arch of the span length 110 cm at different highs of $H=10$, 30, 45 and 55 cm using two percentages of reinforcement of 0.628% and 2.5% The force–deflection curves of the concrete arch were plotted for various retrofitting methods, including: No FRP (N), Internal FRP (I), External FRP (E), and both External and Internal FRP (E&I) at two percentages of reinforcement. It was shown that for the arch with low ratio of the high to span, the Internal FRP is the most effective method. Furthermore, increasing the ratio of reinforcement decreases the efficiency of FRP strengthening. On the other hand, for the arch with the geometry of a semi-circular shape, the External FRP can be an appropriate retrofitting method. Moreover, the mechanism of failure in the arch can be transferred from the flexural failure into the shearing failure as the high to span ratio increases. In addition, a more brittle behavior was observed in the arch with higher ratio of the high to span, where the high value of the reinforcement was utilized.

References

[1] Tao, Y., Stratford, T.J., Chen, J.F. (2011) Behaviour of a masonry arch bridge repaired using fibre-reinforced polymer composites. Engineering Structures 33, 1594–1606.

[2] Alecci, V., Misseri, G., Rovero, L., et al. (2016) Experimental investigation on masonry arches strengthened with PBO-FRCM composite. Composites Part B: Engineering 100, 228–239.

[3] Martin, T., Taylor, S., Robinson, D., Cleland, D. (2019) Finite element modelling of FRP strengthened restrained concrete slabs. Engineering Structures 187, 101–119.

[4] Dagher, H.J., Bannon, D.J., Davids, W.G., Lopez-Anido, R.A. (2012) Bending behavior of concrete-filled tubular FRP arches for bridge structures. Construction and Building Materials 37, 432–439.

[5] Chen, H., Zhou, J., Fan, H., et al. (2014) Dynamic responses of buried arch structure subjected to subsurface localized impulsive loading: Experimental study. International Journal of Impact Engineering 65, 89–101.

[6] Hamed, E., Chang, Z.T., Rabinovitch, O. (2015) Strengthening of reinforced concrete arches with externally bonded composite materials: Testing and analysis. Journal of Composites for Construction 19, 04014031.

[7] Zhang, X., Wang, P., Jiang, M., et al. (2015) CFRP strengthening reinforced concrete arches: Strengthening methods and experimental studies. Composite Structures 131, 852–867.

[8] Pramono, E., Willam, K. (1989) Fracture energy-based plasticity formulation of plain concrete. Journal of Engineering Mechanics ASCE 115, 1183–1204.

[9] Feenstra, P.H., De Borst, R. (1996) A composite plasticity model for concrete. International Journal of Solids and Structures 33, 707–730.

[10] Li, Y.F., Lin, C.T., Sung, Y.Y. (2003) A constitutive model for concrete confined with carbon fiber reinforced plastics. Mechanics of Materials 35, 603–619.

[11] Červenka, J., Papanikolaou, V.K. (2008) Three dimensional combined fracture–plastic material model for concrete. International Journal of Plasticity 24, 2192–2220.

[12] Bažant, Z.P., Ožbolt, J. (1990) Nonlocal microplane model for fracture, damage, and size effect in structures. Journal of Engineering Mechanics 116, 2485–2505.

[13] Voyiadjis, G.Z., Abu-Lebdeh, T.M. (1993) Damage model for concrete using bounding surface concept. Journal of Engineering Mechanics 119, 1865–1885.

[14] Lubliner, J., Oliver, J., Oller, S., Oñate, E. (1989) A plastic-damage model for concrete. International Journal of Solids and Structures 25, 299–326.

[15] Lee, J., Fenves, G.L. (1998) Plastic-damage model for cyclic loading of concrete structures. Journal of Engineering Mechanics 124, 892–900.

[16] Faria, R., Oliver, J., Cervera, M. (1998) A strain-based plastic viscous-damage model for massive concrete structures. International Journal of Solids and Structures 35, 1533–1558.

[17] Ibrahimbegović, A., Markovič, D., Gatuingt, F. (2003) Constitutive model of coupled damage-plasticity and its finite element implementation. European Journal of Computational Mechanics 12, 381–405.

[18] Salari, M.R., Saeb, S., Willam, K.J., et al. (2004) A coupled elastoplastic damage model for geomaterials. Computer Methods in Applied Mechanics and Engineering 193, 2625–2643.

[19] Gatuingt, F., Desmorat, R., Chambart, M., et al. (2008) Anisotropic 3D delay-damage model to simulate concrete structures. European Journal of Computational Mechanics 17, 749–760.

[20] Yu, T., Teng, J.G., Wong, Y.L., Dong, S.L. (2010) Finite element modeling of confined concrete-II: Plastic-damage model. Engineering Structures 32, 680–691.

[21] Grassl, P., Jirásek, M. (2006) Damage-plastic model for concrete failure. International Journal of Solids and Structures 43, 7166–7196.

[22] Nguyen, G.D., Korsunsky, A.M. (2006) Damage-plasticity modelling of concrete: calibration of parameters using separation of fracture energy. International Journal of Fracture 139: 325–332.

[23] Hwang, Y.K., Bolander, J.E., Lim, Y.M. (2016) Simulation of concrete tensile failure under high loading rates using three-dimensional irregular lattice models. Mechanics of Materials 101, 136–146.

[24] Nguyen, G.D., Houlsby, G.T. (2008) A coupled damage–plasticity model for concrete based on thermo-dynamic principles: Part I: model formulation and parameter identification. International Journal for Numerical and Analytical Methods in Geomechanics 32, 353–389.

[25] Nguyen, G.D., Houlsby, G.T. (2008) A coupled damage–plasticity model for concrete based on thermo-dynamic principles: Part II: non-local regularization and numerical implementation. International Journal for Numerical and Analytical Methods in Geomechanics 32, 391–413.

[26] Moslemi, H., Khoei, A.R. (2009) 3D adaptive finite element modeling of non-planar curved crack growth using the weighted superconvergent patch recovery method. Engineering Fracture Mechanics 76, 1703–1728.

[27] Khoei, A.R., Moslemi, H., Sharifi, M. (2012) Three-dimensional cohesive fracture modeling of non-planar crack growth using adaptive FE technique. International Journal of Solids and Structures 49, 2334–2348.

[28] Khoei, A.R., Eghbalian, M., Moslemi, H., Azadi, H. (2013) Crack growth modelling via 3D automatic adaptive mesh refinement based on modified–SPR technique. Applied Mathematical Modeling 37, 357–383.

[29] Babuska, I., Melenk, J.M. (1997) The partition of unity method. International Journal for Numerical Methods in Engineering 40, 727–758.

[30] Belytschko, T., Black, T. (1999) Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 45, 601–620.

[31] Khoei, A.R. (2015) Extended Finite Element Method: Theory and Applications. John Wiley.

[32] Broumand, P., Khoei, A.R. (2013) The extended finite element method for large deformation ductile fracture problems with a non-local damage-plasticity model. Engineering Fracture Mechanics 112, 97–125.

[33] Pañeda, E.M., Natarajan, S., Bordas, S. (2017) Gradient plasticity crack tip characterization by means of the extended finite element method. Computational Mechanics 59, 831–842.

[34] Borja, R.I. (2008) Assumed enhanced strain and the extended finite element methods: A unification of concepts. Computer Methods in Applied Mechanics and Engineering 197, 2789–2803.

[35] Chahine, E., Laborde, P., Renard, Y. (2008) Crack tip enrichment in the XFEM using a cutoff function. International Journal for Numerical Methods in Engineering 75, 629–646.

[36] Tarancon, J.E., Vercher, A., Giner, E., Fuenmayor, F.J. (2009) Enhanced blending elements for XFEM applied to linear elastic fracture mechanics. International Journal for Numerical Methods in Engineering 77, 126–148.

[37] Fries, T.P. (2008) A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering 75: 503–532.

[38] Gracie, R., Wang, H.W., Belytschko, T. (2008) Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods. International Journal for Numerical Methods in Engineering 74, 1645–1669.

[39] Shen, Y.X., Lew, A. (2010) An optimally convergent discontinuous Galerkin-based extended finite element method for fracture mechanics. International Journal for Numerical Methods in Engineering 82, 716–755.

[40] Khalfallah, S., Charif, A., Naimi, M. (2004) Nonlinear analysis of reinforced concrete structures using a new constitutive model. European Journal of Computational Mechanics 13, 841–856.

[41] Bertolesi, E., Milani, G., Fedele, R. (2016) Fast and reliable non-linear heterogeneous FE approach for the analysis of FRP-reinforced masonry arches. Composites Part B: Engineering 88, 189–200.

[42] Ziad, N. (2008) Elasto-plastic and damage modeling of reinforced concrete. PhD thesis, Louisiana State University.

[43] Dugdale, D.S. (1960) Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids 8, 100–104.

[44] Barenblatt, G.I. (1962) The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics 7, 55–129.

[45] Shah, S.G., Chandra Kishen, J.M. (2010) Nonlinear fracture properties of concrete–concrete interfaces. Mechanics of Materials 42, 916–931.

[46] Khoei, A.R., Moslemi, H., Ardakany, K.M., Barani, O.R., Azadi, H. (2009) Modeling of cohesive crack growth using an adaptive mesh refinement via the modified–SPR technique. International Journal of Fracture 159, 21–41.

[47] Turon, A., Camanho, P.P., Costa, J., Dávila, C.G. (2006) A damage model for the simulation of delamination in advanced composites under variable-mode loading, Mechanics of Materials 38, 1072–1089.

[48] Wells, G.N., Sluys, L.J. (2001) A new method for modelling cohesive cracks using finite elements. International Journal for Numerical Methods in Engineering 50, 2667–2682.

[49] Gopalaratnam, V., Shah, S.P. (1985) Softening response of plain concrete in direct tension. ACI Materials Journal 82: 310–323.

[50] Au, C., Büyüköztürk, O. (2006) Peel and shear fracture characterization of debonding in FRP plated concrete affected by moisture. Journal of Composites for Construction 10, 35–47.

Biographies

Amir R. Khoei received his Ph.D. in Civil Engineering from the University of Wales Swansea in UK in 1998. He is currently a professor in the Civil Engineering Department at Sharif University of Technology. He is a member of the editorial board in the journals of Finite Elements in Analysis and Design and European Journal of Computational Mechanics. He has been selected several times as a distinguished professor at the Sharif University of Technology, and as a distinguished professor by the Ministry of Science, Research and Technology in 2008. He is the silver medal winner of Khwarizmi International Award organized by the Iranian Research Organization for Science and Technology.

Tahmaz Ahmadpour received his BSc degree in Civil Engineering from Iran University of Science and Technology, Tehran, Iran in 2001. He received his MSc degree in Structural Engineering from Sahand University of Technology, Tabriz, Iran in 2004. He is currently a PhD student in Civil Engineering at the International Campus of the Sharif University of Technology. His main research interests are Computational Fracture and Damage Mechanics for Brittle material and Extended Finite Element Method (X-FEM).

Yousef Navidtehrani received his BSc degree in Civil Engineering from Iran University of Science and Technology, Tehran, Iran in 2014. He received his MSc degree in Structural Engineering from Sharif University of Technology, Tehran, Iran in 2017. He is currently working as a researcher in the University of Oviedo, Spain cooperating with Sharif University of Technology, Iran and Imperial College London, UK. His main research interests are Computational Plasticity, Computational Fracture and Damage Mechanics, Phase Field Modeling, Material Characterization and Multiscale Modeling.

European Journal of Computational Mechanics, Vol. 30_1, 1–50.
doi: 10.13052/ejcm1779-7179.3011
© 2021 River Publishers

Abstract

1 Introduction

2 The X-FEM Governing Equations

2.1 The Enrichment of Displacement Field

2.2 Discretization of Governing Equations

3 The Concrete Plastic-Damage Model

3.1 The Damage Model

3.2 The Plasticity Model and Hardening Function

3.3 A Combined Plastic-damage Model for Strong Discontinuity

4 Cohesive Fracture Model for FRP Interface

4.1 The Damage Criterion and Damage Evolution Law

5 Crack Initiation and Crack Growth Direction Criteria

6 Validation of Computational Algorithm

6.1 X-FEM Modeling of a Three-point Bending Beam with an Edge Crack

6.2 Modeling the Concrete Behavior with a Plastic-damage Model

6.3 Cohesive Fracture Modeling of the Debonding Test

6.4 Numerical and Experimental Investigations of Concrete Arch

7 Parametric Studies of FRP Strengthening of Concrete Arches

7.1 The Concrete Arch with the High of $\mathit{H}$ $=$ 10 cm

7.2 The Concrete Arch with the High of $\mathit{H}$ $=$ 30 cm

7.3 The Concrete Arch with the High of $\mathit{H}$ $=$ 45 cm

7.4 The Concrete Arch with the High of $\mathit{H}$ $=$ 55 cm

8 Conclusion

References

Biographies