A novel strain-based finite element family for mesh-independent analysis of the tensile failure of reinforced concrete bars

Abstract

The paper presents a novel family of strain-based beam finite elements (FE) for analysis of tensile failure of a reinforced concrete bar (RC bar), with results of the analysis being independent of the applied FE mesh. The composite bar consists of a continuous longitudinal ductile reinforcing bar(s) surrounded by brittle concrete cover, which are considered separately in the model. Longitudinal slip at the contact between the concrete cover and reinforcing bars is allowed, while their relative displacements perpendicular to the axis of the RC bar are prevented. Cracks in concrete cover occur when tensile stress in concrete exceeds its tensile strength. Propagation of partially connected crack, that is, softening of the material at the crack, is described through constitutive law in form of nonlinear relationship between stresses in concrete at the crack and the width of the crack. Each separate crack is considered discretely as a discontinuity in geometry of the element. In the analysis of cracking of concrete, it is commonly assumed that the discrete crack can occur at the nodes of FE only. However, this assumption leads to dependence of the analysis results on the employed FE mesh. The presented family of FE enables occurrence of the crack anywhere along the FE. In order to achieve this, the discrete crack is excluded from equations of FE and additional boundary conditions are introduced at the discontinuity. This approach ensures that the location of the cracks, their number and their propagation are independent of the applied FE mesh. Advantages of the novel family of FE are thoroughly presented in a parametric study which investigates influence of number of FE as well as influence of degrees of interpolation and integration on the cracking of RC bar under tensile loading. Experimental results of tensile tests on the analysed bar are available in literature. It can be concluded that the results obtained with the minimal possible number of novel FE and sufficiently high degree of numerical integration scheme, applied for solving integrals in equations of FE, are considerably more accurate than the results of previous analyses with model of discrete crack at the nodes of FE only.

Introduction

Concrete is a heterogeneous and brittle material, where even very small tensile loading can result in occurrence of cracks. Cracking of concrete significantly influences stiffness, ductility and strength of reinforced concrete (RC) structures; therefore, its consideration is essential for realistic mechanical analysis of such structures. Currently, there are two established advanced models of concrete cracking, both based on the principles of fracture mechanics. Cracking of concrete in the first model is considered through occurrence of discrete cracks (e.g. in Bajc et al., 2013; Dias-da-Costa et al., 2009; Federation internationale du beton, 2013; Yang et al., 2008 and others). When the tensile strength of concrete, f_ct, is exceeded, a discontinuity in the geometry of the element occurs (see Figure 1(a)). The next phase in the process of crack opening is when both parts of concrete cover at the crack are still connected by aggregate grains (Cerioni et al., 2011). This phenomenon of partial connection through aggregate grains is in scientific literature known as aggregate bridging. The constitutive law of crack propagation is defined through mutually related tensile stress in the crack, σ^r, and the width of the crack, r. Several researchers formulated the law as a bilinear relationship (e.g. Federation internationale du beton, 2013; Rabzcuk et al., 2005), while others employed the form of exponential function (e.g. Carpinteri, 1999). Fracture energy of concrete, G_f, is a material parameter introduced in the constitutive law of crack propagation, which describes entire work required for failure of cracked part of concrete cover and equals the area under σ^r-r curve (see Figure 1(a)). Bilinear constitutive law of crack propagation depicted in Figure 1(a) is taken from Rabzcuk et al. (2005). The limit value of the crack width, r = w_crit, reached when σ^r decreases to zero, is calculated through given equation. Here, the parameter α_t is a constant and equals 0.14, while the parameter β_t depends on compressive strength of concrete, f_c. When cracks are considered to be discrete, that is, as discontinuities in the geometry of the elements, the structure is usually analysed with finite element method and appearance of cracks is limited to nodes of the elements. Therefore, in order to achieve sufficiently accurate results of crack propagation analysis, which are independent of finite element mesh, sufficiently dense finite element mesh must be employed. Lately, mesh-free methods, such as free Galerkin method (see e.g. Rabzcuk et al., 2008; Yang et al., 2008), have proven to be a promising alternative to the finite element method.

Figure 1.

Constitutive law of propagation of a discrete and a smeared crack.

The second model considers cracking of concrete in a form of so-called smeared cracks, where the distribution of cracks is not realistic, as presented in Figure 1(b) (see Bažant and Pijaudier-Cabot, 1989; Bažant and Planas, 1997; Federation internationale du beton, 2013; Yang et al., 2008 and others). When the tensile strength of concrete is exceeded, localisation of extensional strains in concrete, ɛ_c, occurs on a band of the element with limited length, L_c, as shown in Figure 1(b).The length of this band of concrete, over which the cracks are smeared, is not arbitrary. Instead, it must be defined as a material parameter and depends on fracture energy of concrete, G_f, and constitutive law of crack propagation. The latter is here defined through tensile stress, σ_c, and localized extensional strain, ${\overline{\epsilon }}_{\text{c}}^{\text{r}}$ , which is notional strain of concrete cover in the band of smeared crack L_c. The fracture energy of concrete G_f is calculated as a product of the length L_c and work required for tensile failure of a volume unit of concrete, γ_f. The latter can be represented by the area under stress-strain curve of tensile loaded concrete at the range of softening of the material. The remaining parameters in Figure 1(b) (ɛ_ct1, ɛ_ct2 and limit strain ɛ_ctu) are inherent to the applied model of crack propagation. Notional extensional strain of concrete ${\overline{\epsilon }}_{\text{c}}^{\text{r}}$ and width of actual crack r is related through well-known equation

$${\overline{\epsilon }}_{\text{c}}^{\text{r}}={\epsilon }_{\text{ct}1}+\frac{r}{L_{\text{c}}}$$

(1)

If it is ensured, that the entire considered structure (or at least parts of the structure, where cracks are most likely to occur) is modelled with finite elements, whose lengths are consistent with characteristic length L_c, it can be assumed, that the results of the analysis are completely independent of finite element mesh (Bažant and Planas, 1997; Federation internationale du beton, 2013). However, this usually requires prediction of crack locations in advance and/or finite element meshes with large number of elements. Consequently, the analyses employing the model of smeared crack are either time-costly or give finite element mesh dependent results.

For the study of occurrence of cracks in concrete cover as well as their propagation, RC bar with tensile load is most convenient. The RC bar is composed of continuous longitudinal reinforcing bar and of enveloping concrete cover. Tensile load can be applied through reinforcing bar, if it continues outside of the concrete cover, or through the concrete cover itself. For the analysis of concrete cracking, it is essential, that cracks in concrete cover appears prior to plastification of reinforcing bars. In order to meet this requirement, the area of cross-section of concrete cover must be sufficiently small (see e.g. Wollrab et al., 1996). The paper presents a novel family of strain-based beam finite elements for analysis of tensile failure of RC bar. The main advantage of the novel finite element is that positions, number and development of discrete cracks are independent of the applied finite element mesh. Furthermore, very good analysis results are obtained with the smallest possible number of finite elements.

In the novel method, cracks occur at so-called ‘excluded’ points anywhere along the length of the finite element, with the only criterion for its occurrence being exceeded tensile strength of concrete. The key criterion for the assessment of suitability of novel finite elements is that the position and number of cracks, together with the rest of physical quantities do not vary significantly with an increase of number of finite elements.

Numerical model of RC bar with discrete crack at ‘excluded’ point

The numerical model is presented on an example of RC bar with only one reinforcing bar with diameter ϕ and area of cross-section A_s. The reinforcing bar is positioned in the centroid axis of the concrete cover with cross-section A_c. A RC bar, cracked due to the application of external tensile load, is schematically presented in Figure 2. Its initial length is L, while L^′ denotes the length of the deformed and cracked bar. Basic assumptions of the model are:

1. the size and the shape of the cross-section of the concrete cover as well as of the reinforcing bar do not change during the deformation; Bernoulli’s hypothesis of planar cross-sections is assumed;

2. the concrete cover and the reinforcing bar are modelled separately; deformations of both parts are limited to longitudinal direction; the influence of the radial pressure exerted by the reinforcing bar on the surrounding concrete is neglected;

3. only the longitudinal slip at the contact between the concrete and reinforcing bar is allowed; the slip law is described through relationship between bond stress τ and relative longitudinal displacement at the contact Δ (see Figure 3(a)); the law is multi-linear with possibility of isotropic hardening;

4. concrete in tension is assumed to be linear elastic material with elastic moduli E_c until cracking occurs; the reinforcing bar is made of bilinear elasto-plastic material with elastic moduli E_s and yield strength f_y;

5. the crack in concrete cover occurs when the stress in concrete reaches its tensile strength f_ct;

6. opening/closing of the crack is described in the same way as it is described in the model with discrete crack, that is, with constitutive law of crack propagation in form of relationship between tensile stress in the crack, σ^r, and width of the crack, r (see Figure 1(a) or Figure 3(b));

7. only one crack per each separate strain finite element can occur.

Figure 2.

Uncracked and cracked RC bar. Locations of characteristic points.

Figure 3.

Constitutive law of (a) bond slip (Comite Euro-International Du Beton and Federation International De la Precontraint, 1993) and (b) crack propagation (Rabzcuk et al., 2005).

Hereinafter, governing equations of a RC bar subjected to tensile load previous to crack occurrence are presented first.

Governing equations of a RC bar previous to cracking

Governing equations of the uncracked RC bar, depicted in Figure 2(a), are taken as in Bajc et al. (2018) and include kinematic, equilibrium and constitutive equation, separately for concrete cover and for reinforcing bar, as well as three constraining equations:

1. Concrete cover

$$\begin{array}{ll}\hfill & u_{\text{c}}^{\prime }-{\epsilon }_{\text{c}}=0\hfill \end{array}$$

(2)

$$\begin{array}{ll}\hfill & N_{\text{c}}^{\prime }+p_{\text{X},\text{c}}=0,\hfill \end{array}$$

(3)

$$\begin{array}{ll}\hfill & N_{\text{c}}=N_{\text{c}}^{\text{c}}=E_{\text{c}}{A}_{\text{c}}{\epsilon }_{\text{c}},\hfill \end{array}$$

(4)

• 2. Reinforcing bar

$$\begin{array}{ll}\hfill & u_{\text{s}}^{\prime }-{\epsilon }_{\text{s}}=0\hfill \end{array}$$

(5)

$$\begin{array}{ll}\hfill & N_{\text{s}}^{\prime }+p_{\text{X},\text{s}}=0,\hfill \end{array}$$

(6)

$$\begin{array}{ll}\hfill & N_{\text{s}}=N_{\text{s}}^{\text{c}}=E_{\text{s}}{A}_{\text{s}}{\epsilon }_{\text{s}}({\sigma }_{\text{s}}=E_{\text{s}}{\epsilon }_{\text{s}}\le f_{\text{y}}),\hfill \end{array}$$

(7)

• 3. Constraining equations

$$\begin{array}{ll}\hfill & \Delta \left(x\right)=u_{\text{s}}\left(x\right)-u_{\text{c}}\left(x\right),\hfill \end{array}$$

(8)

$$\begin{array}{ll}\hfill & p_{\text{X},\text{c}}+p_{\text{X},\text{s}}=0,\hfill \end{array}$$

(9)

$$\begin{array}{ll}\hfill & p_{\text{X},\text{c}}=p_{\text{X},\text{c}}^{\text{c}}=\tau \pi \varphi \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\tau =f(\Delta ).\hfill \end{array}$$

(10)

Herein, prime ()^′ denotes a derivative of certain quantity with respect to material coordinate x. Longitudinal displacements of reference axis of concrete cover and reinforcing bar are denoted with u_c and u_s, respectively, while ɛ_c and ɛ_s are inherent extensional strains of reference axes. N_c and N_s are equilibrium axial forces, while $N_{\text{c}}^{\text{c}}$ and $N_{\text{s}}^{\text{c}}$ are constitutive axial forces.

The constraining equations address longitudinal slip and shear flow at the contact of concrete cover and reinforcing bar. p_X,c and p_X,s denote shear components of contact line load in concrete cover and in reinforcing bar (shear flow), respectively, whereas $p_{\text{X},\text{c}}^{\text{c}}$ stands for constitutively defined shear flow in concrete cover and depends on constitutive law of bond slip, given through relationship τ − Δ.

Along with nine governing equations (2)–(10) for nine unknown functions (namely u_c, u_s, ɛ_c, ɛ_s, N_c, N_s, p_X,c, p_X,s and Δ), the corresponding boundary conditions, either static or kinematic, must be given. Boundary conditions at both ends of the concrete cover are

$$\begin{array}{ll}\hfill & \text{either}\hspace{0.17em}\hspace{0.17em}{S}_{\text{c},1}+N_{\text{c}}\left(0\right)=0\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{u}_{\text{c}}\left(0\right)=u_{\text{c},1}\hfill \end{array}$$

(11)

$$\begin{array}{ll}\hfill & \text{either}\hspace{0.17em}\hspace{0.17em}{S}_{\text{c},2}-N_{\text{c}}\left(L\right)=0\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{u}_{\text{c}}\left(L\right)=u_{\text{c},2}\hfill \end{array}$$

(12)

Boundary conditions at both ends of the reinforcing bar are

$$\begin{array}{ll}\hfill & \text{either}\hspace{0.17em}\hspace{0.17em}{S}_{\text{s},1}+N_{\text{s}}\left(0\right)=0\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{u}_{\text{s}}\left(0\right)=u_{\text{s},1}\hfill \end{array}$$

(13)

$$\begin{array}{ll}\hfill & \text{either}\hspace{0.17em}\hspace{0.17em}{S}_{\text{s},2}-N_{\text{s}}\left(L\right)=0\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{u}_{\text{s}}\left(L\right)=u_{\text{s},2}\hfill \end{array}$$

(14)

Here, S_c,i and S_s,i denote external load at each end of concrete cover and of reinforcing bar, while u_c,i and u_s,i represent displacements at each end. Quantities with lower index i = 1 belong to the end with material coordinate x = 0 and quantities with lower index i = 2 belong to the end with material coordinate x = L. In case that the RC bar is loaded in compliance with Figure 2(a), that is, if the stresses in RC bar are imposed only by load P at the unsupported end of the reinforcing bar (S_s,2 = P), the boundary conditions read as follows

$$N_{\text{c}}\left(0\right)=0,N_{\text{c}}\left(L\right)=0,u_{\text{s}}\left(0\right)=0\hspace{0.17em}\hspace{0.17em}\text{in}\hspace{0.17em}\hspace{0.17em}{N}_{\text{s}}9\left(L\right)=P$$

(15)

RC bar with one crack

A crack in concrete cover occurs when the tensile strength of the concrete, f_ct, is reached. After the appearance of the crack, aggregate grains can still connect both ends of the crack, which is a phenomenon also known as aggregate bridging (see Cerioni et al., 2011; Yang and Chen, 2005; and others). This partial connection of concrete cover with aggregate grains applies only when the crack is sufficiently small (r ≤ w_crit). Crack propagation is defined in the form of relationship between tensile stress in the crack, σ^r, and its width, r, as it was already stated in the list of basic assumptions of the model (see Figure 3(b)).

We assume that the crack occurs in concrete cover at material coordinate $x_{\text{c}}^{\text{r}}$ (see Figure 2). Without any loss of generality an element with only one crack is analysed in this section. The crack is considered as an ‘excluded’ point on the length of the element. A material coordinate ${}^{-}{x}_{\text{c}}^{\text{r}}$ is attributed to concrete cover at the left edge of the crack, while on the right edge of the crack the material coordinate of concrete cover is ${}^{+}{x}_{\text{c}}^{\text{r}}$ , as shown in Figure 2(b). In cracked state, the corresponding points on the reinforcing bar have material coordinates ${}^{-}{x}_{\text{s}}^{\text{r}}$ and ${}^{+}{x}_{\text{s}}^{\text{r}}$ . The following relationships apply

$$\hspace{0.17em}{}^{-}{x}_s^r={}^{-}{x}_c^r-\Delta \left({}^{-}{x}_c^r\right)$$

(16)

$$\begin{array}{ll}\hfill & {}^{+}{x}_{\text{s}}^{\text{r}}={}^{+}{x}_{\text{c}}^{\text{r}}-\Delta \left({}^{+}{x}_{\text{c}}^{\text{r}}\right)\hfill \end{array}$$

(17)

where $\Delta \left({}^{-}{x}_{\text{c}}^{\text{r}}\right)$ and $\Delta \left({}^{+}{x}_{\text{c}}^{\text{r}}\right)$ denote slip of concrete cover at each edge of the crack (see Figure 4(b)). The width of the crack is defined through the difference in longitudinal displacements of left and right concrete cover cross-sections at the crack (see Figure 4(c)).

$$r=u_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)-u_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r}}\right)$$

(18)

Figure 4.

Distribution of characteristic quantities along the length L of the RC bar after appearance of the first crack.

Note, that the characteristic quantities along RC bar after appearance of the first crack in Figure 4 are depicted with respect to the undeformed length of the bar L.

The uncracked concrete cover to the left and right of the crack, that is, at intervals $[0,{}^{-}{x}_{\text{c}}^{\text{r}}]$ and $[\hspace{0.17em}{}^{+}{x}_{\text{c}}^{\text{r}},L]$ , still complies with equations (2)–(4). Similarly, equations (5)–(7) continue to apply along the entire length of the reinforcing bar on interval [0, L]. At intervals $[0,{}^{-}{x}_{\text{c}}^{\text{r}}]$ and $[{}^{+}{x}_{\text{c}}^{\text{r}},L]$ , where uncracked concrete cover and reinforcing bar are in continuous contact, constraining equations (8)–(9) also apply, together with the constitutive law of contact, equation (10). Additionally, the following two boundary conditions at the crack also apply

$$\begin{array}{ll}\hfill & x={}^{-}{x}_{\text{c}}^{\text{r}}:-N_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r}}\right)+N^{\text{r}}=0\hfill \end{array}$$

(19)

$$\begin{array}{ll}\hfill & x={}^{+}{x}_{\text{c}}^{\text{r}}:N_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)-N^{\text{r}}=0\hfill \end{array}$$

(20)

N^r denotes axial force at the crack and depends on constitutive law of crack propagation, which is formulated as relationship σ^r − r

$$N^{\text{r}}=N^{\text{r},\text{c}}={\sigma }^{\text{r}}\left(r\right)A_{\text{c}}$$

(21)

In order to determine 11 unknown quantities of tensile RC bar (u_c, u_s, ɛ_c, ɛ_s, N_c, N_s, p_X,c, p_X,s, Δ, r and N^r), nine governing equations (2)–(10) must be solved, together with kinematic relationship (18) and constitutive equation of the crack (21).

Finite element formulation

The novel family of strain finite elements is derived hereinafter on a basis of modified principle of virtual work, written for entire RC bar with one crack. Note, that the discrete crack is considered as an ‘excluded’ point of the finite element. However, a system of discrete generalized equilibrium equations of finite element without crack is briefly presented first. This system of equations is applicable in the model of RC bar with discrete crack at ‘excluded’ point previous to cracking, and has already been presented in detail by Bajc et al. (2013, 2018).

RC bar previous to cracking

The system of discrete generalized equilibrium equation of uncracked RC bar with tensile load can be found in Bajc et al. (2018). Therefore, the equations themselves are not written down in the present paper, only their most important characteristics are presented. The system consists of 2N + 6 algebraic equations and the same number of unknowns, where 2N + 2 of them are internal degrees of freedom, namely, ɛ_c,n, ɛ_s,n, N_c(0), N_s(0), and four of them are external degrees of freedom, namely, u_c (0), u_s (0), u_c(L) and u_s(L). Here, ɛ_c,n and ɛ_s,n are strain quantities at n^th equidistant points along concrete cover or reinforcing bar (n = 1, 2, …, N). Axial deformation of reference axis of concrete cover, ɛ_c, and of reinforcing bar, ɛ_s, are then approximated with Lagrangian polynomials P _n (x) of N − 1^th order

$${\epsilon }_{\text{c}}\left(x\right)={\displaystyle \sum _{n=1}^N{P}_n}\left(x\right){\epsilon }_{\text{c},n},\hspace{0.17em}\hspace{0.17em}{\epsilon }_{\text{s}}\left(x\right)={\displaystyle \sum _{n=1}^N{P}_n}\left(x\right){\epsilon }_{\text{s},n}$$

(22)

Furthermore, all integrals in equations of finite element are calculated numerically, using Lobatto or Gauss integration scheme.

Modified principle of virtual work for a RC bar with one crack

Principle of virtual work states that virtual work of internal forces equals virtual work of external forces. Here, we begin with modified principle as it is given in Bajc et al. (2013, 2018) for RC bar with discrete cracks at nodes of finite element, however, applicable adaptations for RC bar with discrete crack at ‘excluded’ point (at material coordinates ${}^{-}{x}_{\text{c}}^{\text{r}}$ and ${}^{+}{x}_{\text{c}}^{\text{r}}$ somewhere along the length of finite element) are introduced

$$\begin{array}{ll}\delta W^{\ast }=\hfill & {\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}\left(N_{\text{c}}^{\text{c}}-N_{\text{c}}\right)}\delta {\epsilon }_{\text{c}}\text{d}x+\left(N^{\text{r},\text{c}}-N^{\text{r}}\right)\delta r+{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L\left(N_{\text{c}}^{\text{c}}-N_{\text{c}}\right)}\delta {\epsilon }_{\text{c}}\text{d}x+\hfill \\ \hfill & {\displaystyle {\int }_0^{{}^{-}{x}_{\text{s}}^{\text{r}}}\left(N_{\text{s}}^{\text{c}}-N_{\text{s}}\right)}\delta {\epsilon }_{\text{s}}\text{d}x+{\displaystyle {\int }_r\left(N_{\text{s}}^{\text{c}}-N_{\text{s}}\right)}\delta {\epsilon }_{\text{s}}^{\text{r}}\text{d}x+{\displaystyle {\int }_{{}^{+}{x}_{\text{s}}^{\text{r}}}^L\left(N_{\text{s}}^{\text{c}}-N_{\text{s}}\right)}\delta {\epsilon }_{\text{s}}\text{d}x+\hfill \\ \hfill & \text{ constraining expressions }+\text{ boundary conditions }=0.\hfill \end{array}$$

(23)

Integral boundaries $[0,{}^{-}{x}_{\text{c}}^{\text{r}}]$ and $[0,{}^{-}{x}_{\text{s}}^{\text{r}}]$ specify domains of concrete cover and of reinforcing bar on the left side of the crack, and $[{}^{+}{x}_{\text{c}}^{\text{r}},L]$ in $[{}^{+}{x}_{\text{s}}^{\text{r}},L]$ specify their domains on the right side of the crack. A region of r defines a domain of the crack (i.e. its width). A deformation energy of the reinforcing bar from the domain of the crack is transferred to the remaining two integrals, which concern deformation energy of the reinforcing bar (i.e. left and right from the crack), on the assumption that this simplification does not result in significantly lesser accuracy. The modified principle of virtual work can now be written as

$$\begin{array}{ll}\delta W^{\ast }=\hfill & {\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}(N_{\text{c}}^{\text{c}}-N_{\text{c}})}\delta {\epsilon }_{\text{c}}\text{d}x+(N^{\text{r},\text{c}}-N^{\text{r}})\delta r+{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L(N_{\text{c}}^{\text{c}}-N_{\text{c}})}\delta {\epsilon }_{\text{c}}\text{d}x+\hfill \\ \hfill & {\displaystyle {\int }_0^{{\overline{x}}_{\text{s}}^{\text{r}}}(N_{\text{s}}^{\text{c}}-N_{\text{s}})}\delta {\epsilon }_{\text{s}}\text{d}x+{\displaystyle {\int }_{{\overline{x}}_{\text{s}}^{\text{r}}}^L(N_{\text{s}}^{\text{c}}-N_{\text{s}})}\delta {\epsilon }_{\text{s}}\text{d}x+\text{ constraining expressions }+\hfill \\ \hfill & \text{ boundary conditions }=0,\hfill \end{array}$$

(24)

where

$${\overline{x}}_{\text{s}}^{\text{r}}=\frac{1}{2}\left(+x_{\text{s}}^{\text{r}}+\left(\hspace{0.17em}-x_{\text{s}}^{\text{r}}\right)\right)$$

(25)

Appropriate constraining expressions in equation (24), which ensure compatibility with prescribed boundary displacements at the ends of RC bar as well as at the edges of the crack, are obtained through integration of kinematic equations (2) and (5)

• 1. Constraining expressions regarding left and right end of concrete cover

$$\begin{array}{ll}\hfill & u_{\text{c}}(\hspace{0.17em}-x_{\text{c}}^{\text{r}})=u_{\text{c}}\left(0\right)+{\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}{\epsilon }_{\text{c}}}\text{d}x\hfill \end{array}$$

(26)

$$\begin{array}{ll}\hfill & u_{\text{c}}\left(L\right)=u_{\text{c}}(\hspace{0.17em}+x_{\text{c}}^{\text{r}})+{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L{\epsilon }_{\text{c}}}\text{d}x,\hfill \end{array}$$

(27)

• 2. Constraining expression regarding reinforcing bar

$$u_{\text{s}}\left(L\right)=u_{\text{s}}\left(0\right)+{\displaystyle {\int }_0^{{\overline{x}}_{\text{s}}^{\text{r}}}{\epsilon }_{\text{s}}}\text{d}x+\Delta u_{\text{s}}\left({\overline{x}}_{\text{s}}^{\text{r}}\right)+{\displaystyle {\int }_{{\overline{x}}_{\text{s}}^{\text{r}}}^L{\epsilon }_{\text{s}}}\text{d}x=u_{\text{s}}\left(0\right)+{\displaystyle {\int }_0^L{\epsilon }_{\text{s}}}\text{d}x+\Delta u_{\text{s}}\left({\overline{x}}_{\text{s}}^{\text{r}}\right)$$

(28)

Here, $\Delta u_{\text{s}}\left({\overline{x}}_{\text{s}}^{\text{r}}\right)={\epsilon }_{\text{s}}^{\text{r}}r$ denotes change of displacement of reinforcing bar over the crack.The modified principle of virtual work for cracked RC bar (see equation (24)) includes, besides boundary conditions for concrete cover and reinforcing bar at both ends of RC bar (equations (11)–(14)), also boundary conditions for concrete cover at the crack (equations (19)–(20)). Furthermore, it is assumed that at the point of crack occurrence $x_{\text{c}}^{\text{r}}$ the longitudinal slip between concrete cover and reinforcing bar equals zero, such that $x_{\text{c}}^{\text{r}}=x_{\text{s}}^{\text{r}}={\overline{x}}_{\text{s}}^{\text{r}}$ . According to this assumption, together with inherent constraining expressions equations (26)–(28) as well as with boundary conditions equations (11)–(14) and (19)–(20), modified principle of virtual work can be written as

$$\begin{array}{ll}\delta W^{\ast }=\hfill & {\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}\left(E_{\text{c}}{A}_{\text{c}}{\epsilon }_{\text{c}}-N_{\text{c}}\right)}\delta {\epsilon }_{\text{c}}\text{d}x+\left({\sigma }^{\text{r}}\left(r\right)A_{\text{c}}-N^{\text{r}}\right)\delta r+{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L\left(E_{\text{c}}{A}_{\text{c}}{\epsilon }_{\text{c}}-N_{\text{c}}\right)}\delta {\epsilon }_{\text{c}}\text{d}x+\hfill \\ \hfill & {\displaystyle {\int }_0^L\left(E_{\text{s}}{A}_{\text{s}}{\epsilon }_{\text{s}}-N_{\text{s}}\right)}\delta {\epsilon }_{\text{s}}\text{d}x+\left[u_{\text{c}}\left(\hspace{0.17em}-x_{\text{c}}^{\text{r}}\right)-u_{\text{c}}\left(0\right)-{\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}{\epsilon }_{\text{c}}}\text{d}x\right]\delta N_{\text{c}}\left(0\right)+\hfill \\ \hfill & \left[u_{\text{c}}\left(L\right)-u_{\text{c}}\left(\hspace{0.17em}+x_{\text{c}}^{\text{r}}\right)-{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L{\epsilon }_{\text{c}}}\text{d}x\right]\delta N_{\text{c}}\left(\hspace{0.17em}+\hspace{0.17em}{x}_{\text{c}}^{\text{r}}\right)+\left[u_{\text{s}}\left(L\right)-u_{\text{s}}\left(0\right)-{\displaystyle {\int }_0^L{\epsilon }_{\text{s}}}\text{d}x-{\epsilon }_{\text{s}}^{\text{r}}r\right]\delta N_{\text{s}}\left(0\right)-\hfill \\ \hfill & \left[S_{\text{c},1}+N_{\text{c}}\left(0\right)\right]\delta u_{\text{c}}\left(0\right)-\left[S_{\text{c},2}-N_{\text{c}}\left(L\right)\right]\delta u_{\text{c}}\left(L\right)-\hfill \\ \hfill & \left[N^{\text{r}}-N_{\text{c}}\left(\hspace{0.17em}-x_{\text{c}}^{\text{r}}\right)\right]\delta u_{\text{c}}\left(\hspace{0.17em}-x_{\text{c}}^{\text{r}}\right)-\left[-N^{\text{r}}+N_{\text{c}}\left(\hspace{0.17em}+x_{\text{c}}^{\text{r}}\right)\right]\delta u_{\text{c}}\left(\hspace{0.17em}+x_{\text{c}}^{\text{r}}\right)-\hfill \\ \hfill & \left[S_{\text{s},1}+N_{\text{s}}\left(0\right)\right]\delta u_{\text{s}}\left(0\right)-\left[S_{\text{s},2}-N_{\text{s}}\left(L\right)\right]\delta u_{\text{s}}\left(L\right)=0\hfill \end{array}$$

(29)

Equation (29) is written with consideration of relations (4), (7) and (21) in place of constitutive quantities $N_{\text{c}}^{\text{c}}$ , $N_{\text{s}}^{\text{c}}$ and N^r,c.

The only two unknown functions in the expression for δW* are axial deformations of concrete cover and of reinforcing bar (ɛ_c and ɛ_s). The remaining unknown quantities in equation (29) are represented by their boundary values only. Axial deformations of reference axis of concrete cover, ɛ_c, and of reinforcing bar, ɛ_s, are interpolated with Lagrangian polynomials P _n (x) (n = 1, …, N) of N − 1^th order. The same formulation as per finite element prior to crack occurrence applies (see equation (22)). Expressions in equation (22) are used in derivation of expressions for variations of strain quantities δɛ_c and δɛ_s, which are then inserted in equation (29). The expression at variation δr is assumed to be identically satisfied, hence, N^r = σ^r(r)A_c. Furthermore, the expressions at independent variations δɛ_c,n, δɛ_s,n, δN_c(0), δN_s(0), $\delta N_{\text{c}}{(}^{+}{x}_{\text{c}}^{\text{r}})$ , δu_c (0), δu_s (0), $\delta u_{\text{c}}{(}^{-}{x}_{\text{c}}^{\text{r}})$ , $\delta u_{\text{c}}{(}^{+}{x}_{\text{c}}^{\text{r}})$ , δu_c(L) and δu_s(L) must equal zero. Accordingly, a system of discrete generalized equilibrium equations for a finite element with one crack is obtained.

$$\begin{array}{ll}\hfill & g_n={\displaystyle {\int }_0^L\left(E_{\text{c}}{A}_{\text{c}}{\epsilon }_{\text{c}}-N_{\text{c}}\right)}{P}_n\text{d}x=0,\left(n=\mathrm{1,2},\dots ,N\right)\hfill \end{array}$$

(30)

$$\begin{array}{ll}\hfill & g_{N+n}={\displaystyle {\int }_0^L\left(E_{\text{s}}{A}_{\text{s}}{\epsilon }_{\text{s}}-N_{\text{s}}\right)}{P}_n\text{d}x=0,\left(n=\mathrm{1,2},\dots ,N\right)\hfill \end{array}$$

(31)

$$\begin{array}{ll}\hfill & g_{2N+1}=u_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r}}\right)-u_{\text{c}}\left(0\right)-{\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}{\displaystyle \sum _{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{c},n}\text{d}x=0\hfill \end{array}$$

(32)

$$\begin{array}{ll}\hfill & g_{2N+2}=u_{\text{s}}\left(L\right)-u_{\text{s}}\left(0\right)-{\displaystyle {\int }_0^L{\displaystyle \sum _{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{s},n}\text{d}x-{\displaystyle \sum _{n=1}^N{P}_n}\left({\overline{x}}_{\text{s}}^{\text{r}}\right){\epsilon }_{\text{s},n}r=0\hfill \end{array}$$

(33)

$$\begin{array}{ll}\hfill & g_{2N+3}=u_{\text{c}}\left(L\right)-u_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)-{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^L{\displaystyle \sum _{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{c},n}\text{d}x=0\hfill \end{array}$$

(34)

$$\begin{array}{ll}\hfill & g_{2N+4}=-S_{\text{c},1}-N_{\text{c}}\left(0\right)=0\hfill \end{array}$$

(35)

$$\begin{array}{ll}\hfill & g_{2N+5}=-S_{\text{s},1}-N_{\text{s}}\left(0\right)=0\hfill \end{array}$$

(36)

$$\begin{array}{ll}\hfill & g_{2N+6}=-{\sigma }^{\text{r}}\left(r\right)A_{\text{c}}+N_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r}}\right)=0\hfill \end{array}$$

(37)

$$\begin{array}{ll}\hfill & g_{2N+7}={\sigma }^{\text{r}}\left(r\right)A_{\text{c}}-N_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)=0\hfill \end{array}$$

(38)

$$\begin{array}{ll}\hfill & g_{2N+8}=-S_{\text{c},2}+N_{\text{c}}\left(L\right)=0\hfill \end{array}$$

(39)

$$\begin{array}{ll}\hfill & g_{2N+9}=-S_{\text{s},2}+N_{\text{s}}\left(L\right)=0\hfill \end{array}$$

(40)

For given external load, the equations (30)–(40) represent a system of 2N + 9 algebraic equations for the same number of unknowns. 2N + 3 of the unknowns are internal degrees of freedom, namely, ɛ_c,n, ɛ_s,n, N_c(0), N_s(0), $N_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)$ , while the remaining six unknowns include external degrees of freedom of the finite element, namely, u_c (0), u_s (0), $u_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r}}\right)$ , $u_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r}}\right)$ , u_c(L) and u_s(L). Functions u_c(x), u_s(x), N_c(x) and N_s(x) are calculated with the following equations.

$$\begin{array}{ll}\hfill & u_{\text{c}}\left(x\right)=\{\begin{array}{l}{u}_{\text{c}}\left(0\right)+{\displaystyle {\int }_0^x{\displaystyle {\sum }_{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{c},n}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{<}^{-}{x}_{\text{c}}^{\text{r}}\hfill \\ u_{\text{c}}\left(0\right)+{\displaystyle {\int }_0^x{\displaystyle {\sum }_{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{c},n}\text{d}x+r,\hspace{1em}\text{if}\hspace{0.17em}x{\ge }^{+}{x}_{\text{c}}^{\text{r}}\hfill \end{array}\hfill \end{array}$$

(41)

$$\begin{array}{ll}\hfill & u_{\text{s}}\left(x\right)=\{\begin{array}{l}{u}_{\text{s}}\left(0\right)+{\displaystyle {\int }_0^x{\displaystyle {\sum }_{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{s},n}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{<}^{-}{x}_{\text{c}}^{\text{r}}\hfill \\ u_{\text{s}}\left(0\right)+{\displaystyle {\int }_0^x{\displaystyle {\sum }_{n=1}^N{P}_n}}\left(x\right){\epsilon }_{\text{s},n}\text{d}x+{\displaystyle {\sum }_{n=1}^N{P}_n}\left({\overline{x}}_{\text{s}}^{\text{r}}\right){\epsilon }_{\text{s},n}·r,\hspace{1em}\text{if}\hspace{0.17em}x{\ge }^{+}{x}_{\text{c}}^{\text{r}}\hfill \end{array}\hfill \end{array}$$

(42)

$$\begin{array}{ll}\hfill & N_{\text{c}}\left(x\right)=\{\begin{array}{l}{N}_{\text{c}}\left(0\right)-{\displaystyle {\int }_0^x{p}_{\text{X},\text{c}}}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{<}^{-}{x}_{\text{c}}^{\text{r}}\hfill \\ N_{\text{c}}\left(0\right)-{\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}{p}_{\text{X},\text{c}}}\text{d}x-{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^x{p}_{\text{X},\text{c}}}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{\ge }^{+}{x}_{\text{c}}^{\text{r}}\hfill \end{array}\hfill \end{array}$$

(43)

$$\begin{array}{ll}\hfill & N_{\text{s}}\left(x\right)=\{\begin{array}{l}{N}_{\text{s}}\left(0\right)-{\displaystyle {\int }_0^x{p}_{\text{X},\text{s}}}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{<}^{-}{x}_{\text{c}}^{\text{r}}\hfill \\ N_{\text{s}}\left(0\right)-{\displaystyle {\int }_0^{{}^{-}{x}_{\text{c}}^{\text{r}}}{p}_{\text{X},\text{s}}}\text{d}x-{\displaystyle {\int }_{{}^{+}{x}_{\text{c}}^{\text{r}}}^x{p}_{\text{X},\text{s}}}\text{d}x,\hspace{1em}\text{if}\hspace{0.17em}x{\ge }^{+}{x}_{\text{c}}^{\text{r}}\hfill \end{array}\hfill \end{array}$$

(44)

Again, all integrals in equations of finite elements are evaluated numerically with the Lobatto or Gauss integration scheme.

To sum up, the occurrence of the crack in finite element leads to increase of the system of discrete generalized equilibrium equations from 2N + 6 to 2N + 9 degrees of freedom. At the same time, the forms for interpolation of deformations ɛ_c and ɛ_s, as well as the integration schemes, remain unchanged. Since the crack can occur in any of finite elements along the RC bar, with the only criterion for its occurrence being reached tensile strength of concrete cover, we denote all finite elements with unanimous denotation $E_{p-t}^{\text{crack}}$ , regardless of the presence of the crack. Here, index p stands for a degree of interpolation polynomial for axial deformations of concrete cover ɛ_c and reinforcing bar ɛ_s, and index t defines a prescribed degree of integration scheme.

Solving algorithm

A RC bar is now modelled with arbitrary finite element mesh. Symbol EL stands for the number of finite elements, and symbol NOD = EL + 1 for a number of nodes. In the finite element, where the crack appears during the numerical analysis, the equations for cracked finite element (as presented in previous sections of this paper) are applied. In the remaining finite elements, the equations for un-cracked finite element (as in Bajc et al. (2013; 2018)) are used. All equations are then combined into the equation of system through well established procedures of numerical theory of structures

$$\mathbf{G}\left(\mathbf{x},\lambda \right)=\mathbf{R}\left(\mathbf{x}\right)-\lambda \mathbf{P}=0$$

(45)

x in equation (45) represents vector of nodal unknowns. Parts of vector x are: (1) strain quantities along the concrete cover and reinforcing bar of separate finite elements, ${\epsilon }_{\text{c},n}^{el}$ and ${\epsilon }_{\text{s},n}^{el}$ (el = 1, …, EL); (2) static quantities at the beginning of each element, $N_{\text{c}}^{el}\left(0\right)$ and $N_{\text{s}}^{el}\left(0\right)$ ; (3) static quantities at edges of cracks, $N_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r},k}\right)$ (with k ≤ EL); (4) longitudinal displacements in nodes, $u_{\text{c}}^{nod}$ and $u_{\text{s}}^{nod}$ (where nod = 1, …, NOD) and (v) longitudinal displacements of concrete cover at the edges of cracks, $u_{\text{c}}\left({}^{-}{x}_{\text{c}}^{\text{r},k}\right)$ and $u_{\text{c}}\left({}^{+}{x}_{\text{c}}^{\text{r},k}\right)$ . λ is a loading factor of the system. With static and kinematic boundary conditions being satisfied, nonlinear system (45) can be solved by applying Newton incremental-iterative method.

Appropriate computer programme for analysis of cracked RC bar was developed in software environment Matlab (The MathWorks, Inc. (2016)). Elements of tangent stiffness matrix of finite element both with and without crack were derived analytically and their exactness was checked with software tool Mathematica (Wolfram Research, Inc. (2017)).

Parametric study on accuracy of the introduced strain finite element

Suitability and accuracy of the novel family of strain finite elements are hereinafter presented on example of numerical analysis of tensile failure of RC bar. Characteristics of the bar as well as applied boundary conditions are consistent with experimental data on tensile test from literature (Wollrab et al. (1996)). The same RC bar has already been analysed by Bajc et al. (2018) with their model of discrete cracks in nodes of finite elements.

Geometrical, material and loading data of the analysed reinforced concrete bar

A RC bar with rectangular cross-section b/h = 12.7/5.08 cm and with initial length L = 63.5 cm is analysed. The bar is reinforced with seven reinforcing bars with diameter 6 mm, which are evenly distributed along the middle of the cross-section, as depicted in Figure 5. Experimental results of tensile test were presented by Wollrab et al. (1996). The load on the tested bar was applied to concrete cover through steel fixtures at both ends of the bar. Two notches, each 2 mm wide and 10 mm deep were made at x = L/2 along the shorter sides of the bar, so as to impose the assumed location of the crack. Above each notch, an extensometer with a gauge length l_e = 1.27 cm was installed, in order to measure relative displacements $r^{\ast }=l_{\text{e}}^{\prime }-l_{\text{e}}$ at the crack. Nevertheless, several cracks appeared at various locations along the RC bar during the tensile test, therefore, elongation of the bar $\Delta u_{\text{c}}=l_{\text{g}}^{\prime }-l_{\text{g}}$ was also measured on a longer scale with extensometer of gauge length l_g = 25.4 cm. Relevant geometrical and material data are shown in Figure 5.Mechanical properties of the contact between the reinforcing bar and the concrete cover are not given in Wollrab et al. (1996). Thus, contact law depicted in Figure 3(a), which is a modification of constitutive law of bond slip as defined in Comite Euro-International Du Beton and Federation International De la Precontraint (1993), is applied in our numerical analyses. The following values of the parameters apply in the model: τ₁ = 0.4, τ_u = 0.8, τ₂ = 0.33 kN/cm², Δ₁ = 0.03, Δ₂ = 1.0, Δ₃ = 3.0 and Δ₄ = 10 mm. Constitutive law of crack propagation is applied in accordance with recommendations by Rabzcuk et al. (2005), as shown in Figure 3(b), and is the same as has already been used by Bajc et al. (2018). The parameters are used with the following values: f_ct = 0.319 kN/cm², G_f = 80 N/m, α_t = 0.14, β_t = 0.2485 and w_crit = 0.129 mm.

Figure 5.

Geometrical, material and loading data of the analysed RC bar (Wollrab et al., 1996).

Analysis of concrete cracking of Wollrab RC bar using novel finite element with discrete crack at ‘excluded’ point

For the analysis of cracked RC bar with model of discrete crack (Bajc et al., 2018), where cracks can occur at nodes of finite element mesh only, 48 finite elements were used for modelling of the central part of the bar (i.e. between the steel fixtures), as presented in Figure 6(a). For the present numerical analysis, the central part of the bar is modelled with nine finite elements of type $E_{4-5}^{\text{crack}}$ . The number of finite elements has been chosen the same as the number of calculated cracks in stabilized state (i.e. 9, see Figure 6(b)). This way it is possible to detect all cracks with current formulation of the novel finite element presented in previous sections, where only one crack per element can occur. Additionally, Figure 6(b) shows the order of appearance of cracks; starting from the first two identical cracks, which are denoted with r₁ and appear next to the steel fixtures, towards the last two identical cracks r₅.

Figure 6.

Comparison of the results: (a) numerical model and finite element mesh, (b) computed distribution of cracks in stable state.

Numerically and experimentally determined elongations of the RC bar in the middle of the bar over the length l_g = 25.4 cm in dependence of tensile load P are compared in Figure 7(a). The results of both numerical models are almost identical, which has been expected, since both models employ the same equations but differ importantly in the applied method of calculation. Furthermore, Figure 7(b) depicts comparison of calculated and measured relative displacements of RC bar at the location of the notch in the middle of the bar. Again, very similar results of both numerical models can be observed.The model with discrete crack (Bajc et al., 2018) requires relatively large number of finite elements, depending on the number of cracks and their positions. On the other hand, the novel family of strain finite elements for analysis of tensile failure of RC bar, where the crack can occur at ‘excluded’ point anywhere along the length of the finite element, gives almost identical results with the smallest possible number of the elements (i.e. one element per one crack). The main advantage of the novel family of finite elements, which is independence of numerical results of the prescribed finite element mesh, is thoroughly presented hereinafter.

Figure 7.

Comparison of the numerical results with experimental results from literature: (a) load–displacement curve on the length l_g, (b) load–displacement curve on the length l_e.

On accuracy of novel family of finite elements

In the parametric study on accuracy of novel family of strain finite elements, the same RC bar as in previous section is analysed. The central part of the bar between the steel fixtures is here modelled with various numbers of finite elements (5, 9, 15, 30, 60 and 100) with equal lengths and of various element types $E_{p-t}^{\text{crack}}$ , as shown in Figure 8(a)–(f). Furthermore, the central part of the bar is also modelled with 15 finite elements of type $E_{4-5}^{\text{crack}}$ with various lengths, as depicted in Figure 8(g). In the parametric study, numerical integration with Lobatto integration scheme is applied for evaluation of all integrals in equations of finite elements.

Figure 8.

Considered finite element meshes.

Influence of number of finite elements

First, state of the bar immediately after the occurrence of the crack r₂ at the location of the notch in the middle of the RC bar is analysed (see Figure 6(b)). Influence of number and type of the finite elements on elongation of concrete cover of RC bar (Δu_c = L^′ − L) is presented in Figure 9(a), while Figure 9(b) depicts influence of number and type of the finite elements on the width of the crack r₂. The chosen reference values are numerical results obtained with 100 finite elements of type $E_{8-9}^{\text{crack}}$ . It can be noted that the accuracy of the results increases with the increase of number of finite elements. Very accurate results are obtained with 30 (or more) finite elements, where the error is smaller than 0.015%. The only exception is the finite element of type $E_{2-3}^{\text{crack}}$ , that is, the element with the lowest degrees of interpolation polynomial and Lobatto integration method. The accuracy of the results can mainly be improved with the use of higher degree of numerical integration t, while the degree of the interpolation polynomial has lesser influence on the accuracy of the results. This observation is clearly visible at smaller numbers of finite elements. Furthermore, results obtained with 15 finite elements of type $E_{4-5}^{\text{crack}}$ with various lengths (see Figure 8(g)) are almost identical to the results obtained with the same number of finite elements of the same type and with equal lengths. Accuracy of the results obtained with the analysis of RC bar with the model of discrete crack according to (Bajc et al., 2018), where 48 finite elements are employed (denoted as $48E_{4-5}^{\text{discr}.}$ ), is additionally depicted in Figure 9 for comparison.

State of the bar immediately after the occurrence of the last crack r₅, that is, when the stabilized state of cracks is established (see Figure 6(b)), is analysed next. As before, influence of number and type of the finite elements on elongation of concrete cover of RC bar is presented in Figure 10(a). On the other hand Figure 10(b), now depicts influence of number and type of the finite elements on the location of the crack r₅. Again, the chosen reference values are numerical results obtained with 100 finite elements of type $E_{8-9}^{\text{crack}}$ . The smallest number of finite elements, with which the last cracks r₅ (8th and 9th crack) can be observed is nine. Results with finite element of type $E_{2-3}^{\text{crack}}$ are not presented due to some convergence problems during calculation, which appeared after occurrence of multiple cracks as a result of to low accuracy of $E_{2-3}^{\text{crack}}$ element in combination with Lobatto integration scheme. Otherwise, observations are similar as at analysing the state of the bar after occurrence of r₂: (1) increase of number of finite elements generally results in increased accuracy of the results and (2) higher degree of numerical integration t is reflected in much higher accuracy of the results than higher degree of interpolation polynomial p, as can be seen at smaller numbers of finite elements. However, note that the accuracy of the elongation of the concrete cover obtained with nine finite elements exceeds the accuracy at 15 finite elements (see Figure 10(a)). This is a consequence of occurrence of certain number of cracks and different integration processes over cracked and uncracked finite elements; with finite element mesh with nine elements one crack occurs in each of the elements. Integrals over all finite elements are thus evaluated separately for part of the element on the left side of the crack and part of the element on the right side of the crack, which leads to high accuracy of the integration process. On the other hand, with nine cracks and 15 finite elements in the finite element mesh, nine elements in stabilized state are cracked and six remain uncracked. The integration in uncracked finite elements is calculated over their undivided length and consequently less accurately, which leads to lesser accuracy of the finale results. The lesser accuracy of integration over the entire undivided length of the uncracked finite elements becomes irrelevant with the use of sufficiently high number of (consequently shorter) finite elements, as can be seen from results for 30 or more finite elements. On the basis of results depicted in Figure 10, it can also be confirmed that modelling of central part of the RC bar with finite elements of different lengths has little influence on the accuracy of the results (see denotation $15E_{4-5}^{\text{crack}}$ ). Again, accuracy of the results obtained with the analysis of RC bar with the model of discrete crack according to (Bajc et al., 2018), denoted as $48E_{4-5}^{\text{discr}.}$ , is additionally depicted in Figure 10 for comparison. As already mentioned, the model of discrete crack enables the occurrence of cracks at the nodes of finite elements only. Consequently, both depicted quantities, elongation of the concrete cover and even more prominently the location of the crack r₅, obtained with $48E_{4-5}^{\text{discr}.}$ finite elements are calculated with notably lesser accuracy.

Figure 9.

Influence of different number of finite elements on: (a) elongation of concrete part of RC bar and (b) crack width r₂ at x = L/2; immediately after appearance of crack r₂.

Figure 10.

Influence of different number of finite elements on: (a) elongation of concrete part of RC bar and (b) relative position of crack r₅; immediately after appearance of the last crack r₅.

Influence of numerical integration

The second part of parametric study on accuracy of the novel family of strain finite elements is dedicated to analysis of influence of the degree and the scheme of numerical integration over the length of finite elements. As before, state of the RC bar after the appearance of the crack r₂ at the notch in the middle of the bar is analysed first. Figure 11(a) depicts influence of degree of numerical integration on the elongation of the concrete cover (Δu_c = L^′ − L), while in Figure 11(b), its influence on the width of the crack r₂ is presented. The results are shown for five different finite element meshes, namely, with 5, 9, 15, 30 and 60 finite elements of equal length, and for two integration schemes, namely, Lobatto and Gauss integration schemes, with indices L and G, respectively. Similar as in previous sections the chosen reference values are numerical results obtained with 100 finite elements of type $E_{8-9}^{\text{crack}}$ . The results show that the increase of degree of numerical integration t results in increase of accuracy of the results for all applied finite element meshes. Furthermore, it can be noted that the Gauss integration scheme returns more accurate results.

Figure 11.

Influence of order of numerical integration on: (a) elongation of concrete part of RC bar and (b) crack width r₂ at x = L/2; immediately after appearance of crack r₂.

Finally, the state of the bar after the occurrence of the last crack r₅ is analysed. Figure 12(a) shows influence of degree of numerical integration on the elongation of the concrete cover and Figure 12(b) depicts its influence on the location of the crack r₅. Again, the chosen reference values are numerical results obtained with 100 finite elements of type $E_{8-9}^{\text{crack}}$ . The increase of the accuracy of the results as a consequence of increased degree of numerical integration t is confirmed for all analysed finite element meshes (with 9, 15, 30 and 60 finite elements). However, higher accuracy of the results due to the use of Gauss integration scheme is in stabilized state less prominent than immediately after occurrence of crack r₂ (see Figure 12(b)).

Figure 12.

Influence of order of numerical integration on: (a) elongation of concrete part of RC bar and (b) relative position of crack r₅; immediately after appearance of the last crack r₅.

Conclusions

The paper presented a novel family of strain-based beam finite elements for analysis of tensile failure of a composite bar. The latter consists of longitudinal reinforcement surrounded by brittle concrete cover. Cracking of the cover occurs when tensile stress exceeds tensile strength of concrete. In the present model, a crack is considered as a discontinuity in a geometry of the element, while reinforcing bar bridges the crack continuously. Material softening of concrete in tension has also been considered at the crack, due to partial connection of concrete cover. The advantage of presented elements is that the crack can occur anywhere along the finite element and not necessarily at its nodes. Furthermore, the results of analysis are independent on the applied finite element mesh. In order to achieve that independence, the crack is excluded from equations of finite element and additional boundary conditions at the location of the crack are introduced. Parametric study aimed at analysis of crack occurrence and crack propagation in tensile loaded reinforced concrete bar and its results were compared to experimental results from literature. It has been shown that:

1. analysis with the lowest possible number of the novel beam finite elements $\left(9E_{4-5}^{\text{crack}}\right)$ where the crack occurs anywhere along the length of the element, gives practically equal results as the analysis with 48 finite elements $\left(48E_{4-5}^{\text{discr}.}\right)$ , where discrete crack occurs at the nodes of the finite element mesh only,

2. higher number of novel finite elements $\left(E_{p-t}^{\text{crack}}\right)$ in a finite element mesh leads to increased accuracy of the results,

3. in analysis with lower number of novel finite elements $\left(E_{p-t}^{\text{crack}}\right)$ , a degree of numerical integration t (applied in solving integrals in finite element equations) has more influence on the accuracy of the results than higher degree of interpolation polynomial p (for interpolation of strain quantities along the element),

4. results of analysis with finite element mesh of 15 FE with different lengths are perfectly comparable with results of analysis with same number of FE with equal lengths; the conclusion can thus be made that the use of finite element mesh with FE with different lengths for the analysis has no significant influence on the accuracy of the results.

The presented model enables efficient analysis of crack occurrence and crack propagation in RC elements with centric tensile load (such as tensile RC hangers or geotechnical anchors). As such, the model represents conceptual basis for development of the model for analysis of cracking, softening and tensile failure of brittle materials subjected to combined axial and bending load. So far, such problems in beam theory have been approached and more or less efficiently solved through the concept of plastic hinges or through the model of smeared crack (crack band element). An upgrading of the present model for possibility of occurrence of several cracks in one finite element will also be part of our future research.