# Stochastic Lattice Modelling of the Force-Displacement and Cracking Behaviour of the Steel Reinforced Tie

Authors: Darius Zabulionis; Ona Lukoševičienė; Rimantas Kačianauskas; Liudas Tumonis; Romualdas Kliukas
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The stochastic modelling of the microcracking and the force-displacement behaviour of the tensile steel reinforced tie using the lattice model is presented in the current article. The three-dimension problem of the modelling of the tie is reduced to the two-dimensional so as the main stiffness parameters of the concrete and the reinforcement of the two-dimensional model would be the same as for the three-dimensional. The concrete and steel obey the Hook law. All elastic constants, as well as dimensions of the tie, were assumed as the deterministic quantities except for the critical concrete tensile strains which were treated as a two-dimensional stationary uncorrelated truncated Gaussian random field. The discrete element approach and the explicit integration scheme have been used for the modelling. The estimations of the main parameters of the force-displacement behaviour stochastic process and other statistical indexes were obtained using 72 realization of the force-displacement behaviour of a chosen model. Extra two stochastic realizations of the two different models, as well as three deterministic models, were modelled to compare stochastic and deterministic behaviour of the force-displacement behaviour. The analysis showed that the force-displacement behaviour of the tie under tensile force cannot be treated as a Gaussian stochastic process when the p value is 0.05 at the small displacements and within the interval when the cracking of the concrete is very intensive. However, at the bigger displacements, when the cracking becomes less intensive, the tensile force can be treated as a Gaussian random variable.


---

## Body

## 1. Introduction

The cracking and the force-displacement behaviour of the reinforced concrete structures is still under intensive consideration for several decades, see Hegemier, Murakami and Hageman [1], Balevičius and Augonis [2], Choi and Cheung [3], Yankelevsky, Jabareen and Abutbul [4], Stramandinoli and Rovere [5], Morelli et al. [6]. As majority phenomenon in nature, the behaviour of the reinforced tie subject to the tensile force is a stochastic phenomenon depending on various random factors: moduli of elasticity, Poison’s ratio, the compressive and tensile strength of the concrete and steel, and the bond between the concrete and the steel rod, etc. The compressive strength is the most investigated random phenomenon of concrete. Other random properties such as the tensile strength, moduli of elasticity, and fracture strains are much less investigated. The knowledge on the stochastic force-displacement behaviour of a steel reinforced tie is even less understood due to lack of the investigations on this subject. There are several theoretical approaches to model the force-displacement behaviour of the tie: by solving the bond-slip behaviour model analytically, see for example Zabulionis et al. [7]; Hegemier, Murakami and Hageman [1], or numerically, see for example Pamin and de Borst [8]. The last approach suggests a wide choice of the numerical methodologies used in the solid mechanics. The disadvantage of the analytical or semianalytical approaches is that the cracking of the concrete can be taken into account very roughly. The numerical approaches allow us to take into account the heterogeneity of the concrete and to model the fracture process. Direct modelling of the mechanical fracture behaviour of the heterogeneous solids and particularly concrete by using the stochastic finite element method is a widely accepted technique [9, 10]. In general, each particular concrete damage model is problem-dependent and presents a combination of various details concerning the field variables described by the constitutive models, the mechanisms of material degradation, and the removing of the lattice elements in accordance with the failure criteria. Smeared or discrete cracks, splitting, breaking, or removing the material are just illustrative samples of the failure modes. The extended finite element method, the interelement crack method, and the element deletion method are reviewed by Song et al. [11]. The usage of the finite element method in the stochastic analysis can be restricted due to the complexity of the method, limited time, and computation resources. Therefore, various alternatives are applied to model fracture of the reinforced concrete members. One of them is the so-called lattice model, which can be used for the modelling of the behaviour of reinforced solids; see Breysse, Fokwa, and Drahy [12]. Historically, lattice model was introduced for homogeneous materials; for instance, see Hrennikoff [13], and later this method was applied for cement-based materials and damage problems; for instance, see Sawamoto et al. [14], Bolander and Hikosaka [15], Kozicki and Tejchman [16], and Lilliu, G. & van Mier [17]. A review of the hybrid techniques combining lattice approach with the discrete particle model, molecular dynamics, and PeriDynamic is given in Pan et al. [18]. It was observed here that despite a large variety, the spring stiffness is still the main concern in many newly developed models. It should be noticed that in these and majority of other works the connecting element, i.e., a rod or a spring, of the lattice model represents a particular amount of homogenous material say filler or binder. Thus, the connecting element represents only the filler or binder. According to other approaches, heterogeneous material is modelled as a web of connecting elements that undertake only axial force and represent particular volume consisting of different materials, said filler and binder, for instance, see Potyondy and Cundall [19], Rojek et al. [20], Pilkavičius, Kačianauskas and Norkus [21], Zabulionis et al. [22]. According to this approach, the elastic constants of the connecting element, mainly only axial stiffness, are calculated by taking into account the interaction of the particles via the matrix and by taking into account different mechanical properties of the materials that constitute the representative volume of the connecting element of the lattice model. Another way to obtain the required elastic constants of the connecting element is to use methodologies developed for homogeneous media; for instance, see Ostoja-Starzewski [23]. However, for this approach, the effective properties of the heterogeneous material must be known, and these properties may be calculated by analytical methodologies: Hashin and Shtrikman [24], Mori and Tanaka [25], method of conditional moments Nazarenko et al. [26]. Or by using the numerical homogenization methodologies, see Pellegrino, Galvanetto and Schrefler [27], Węglewski et al. [28]. An advantage of the lattice method is the ability to model fracture process of the materials, see Kosteski, D’Ambra and Iturrioz [29], Braun and Fernández-Sáez [30].The bond between concrete and reinforcement as well as cracking of the concrete has a crucial influence on the bearing capacity and the deformation behaviour of the reinforced concrete structures. It is well known that such phenomena as deflections and cracking of the reinforced concrete structures have a big variability that has to be taken into account in practical usage of the structures to provide the required reliability and safety. Despite the fact that there is a huge amount of the experimental and theoretical investigations of the force-displacement and tension stiffening behaviour of the reinforced concrete ties, there is very few works of literature in which the force-displacement behaviour of the reinforcement tie is analysed as a stochastic process. Therefore, the identification of the force-displacement stochastic process and the estimation of the parameters of the process are actual issues.In the present article, the stochastic modelling of the two-dimensional computational model of the steel reinforced tie is presented when the critical concrete tensile strain is a truncated normal random variable while other properties of the tie, concrete and steel, that is moduli of elasticity, Poison’s ratios as well as the dimensions of the two-dimensional computational model are deterministic quantities. The behaviour of the tie tensile force dependent on the displacement of the reinforcement is the principal object of the investigation. The investigation has been conducted by modelling six identical models of the tie except for the critical concrete tensile strains: three deterministic models were modelled with different deterministic (fixed) critical concrete tensile strains, and three stochastic models were modelled with spatially distributed critical concrete tensile strains. Stochastic properties of one model were estimated using 72 realizations. Number of realizations of remaining stochastic models is one. The novelty of the article is the identification and characterisation of the force-displacement stochastic process of the reinforced tie, also the estimation of the parameters of the stochastic process such as the range, coefficient of correlations and so on.


## 2. Modelling Object

### 2.1. Geometry and Properties of the Reinforced Tie

The geometry of the tie is shown in Figure1. The tie is rectangular with the steel reinforcement placed in the middle of the tie. The initial length of the tie (the initial length of the steel reinforcing rod and concrete part of the tie)L0 = 56.80 mm. The height and the width of the model are of the same size:H=b = 67.9 mm, Figure 1(a). The diameter of the steel reinforcement is adopted as 14 mm.Figure 1
Geometry of the steel: general view of the tie (a), a cross-section of the tie (b).
(a) (b)The left face of the end of the tie is fixed. The displacementu is applied to the free end of the reinforcement from the right, Figure 1(a). The concrete comprises the hardened cement paste and the aggregates assumed as the spherical particles of diameter 0.35 mm embedded in the paste. It is assumed that the moduli of elasticity of the steel reinforcement and concrete are the following: Es=210 GPa, Ec=9.954 GPa. The Poison coefficients are νc=0.2 for concrete and νs=0.33 for the steel reinforcement, respectively.


### 2.2. Computational Model of the Tie

To apply the planar equilateral triangular elastic spring lattice method for the tie, the 3D tie model is reduced to the 2D model consisting of the concrete and the steel reinforcement. The geometry of the 2D model is shown in Figure2. The thickness of the 2D model equals the diameter of the aggregate, i.e., δ = 2Rp=0.35 mm, Figure 2(b). The initial length of the 2D reinforcement and concrete parts are the same as of 3D model L0=56.80 mm and the height of the modelH= 67.9 mm, Figure 2(a). The ribs of the steel reinforcement are trapezoidal with height tr=1.73 mm, and the length cII=6.00 mm, Figure 2(d). The space between the ribs is equal to cI=8.2 mm. It is assumed that the displacements of the left plane of the tie is restrained fully, that is, for the left plane the displacements ux=uy=0.Figure 2
Geometry, boundary conditions, and action of the 2D model of the tie: a scheme of the 2D model (a), a section of the 2D model (b), a geometry of the computational model of the tie (c), and a model of the steel part (d).


### 2.3. Geometry of the Computational Model

For the simulation of the tie, the half of the tie is modelled by using the lattice methodology. The geometry of the model used for the simulation is shown in Figure2(c). The simulation model consists of one part of the concrete and the height of the steel reinforcement equals half of the tie. The thickness of the simulation model, the reinforcement, and the concrete part is the same as for the reduced 2D tie model, δ = 2Rp=0.35 mm. The displacements of the left plane of the simulation model is restrained fully, that is, for the left plane the displacements ux=uy=0 and also the displacements of the bottom plane of the steel reinforcement is restrained in directiony, that is, uy=0, Figure 2(c).It is assumed that the concrete of the 2D model of the tie consists of the spherical aggregate particles of diameter 0.35 mm embedded in the hardened cement paste, Figure3. In the plane, the particles form equilateral triangular cells. The distance of the centres of the two any vicinal particles is 0.4 mm. There is only one layer of the particles in the concrete part of the tie. The volume ratio of the particles (the ratio of the volume of the particles to the volume of the entire concrete) is 0.45815.Figure 3
Structure of the concrete part of the 2D model.The cross-section areas and moduli of elasticity of the 2D model are obtained from the following equations:Es,eff·As,eff = Es·As, and Ec,eff·Ac,eff = Ec·Ac, where As,eff, Ac,eff, Es,eff, and Ec,eff are effective cross-section areas and effective moduli of elasticity of the 2D model of the concrete part and 2D reinforcement, respectively.The effective cross-section areaAs,eff of the 2-reinforcement can be obtained from the condition that the stiffness of the steel part of the 2D model(1)Krnf,2D=As,effEs,effL0=14/Krnf,2D,I+2/Krnf,2D,KII+6/Krnf,2D,IIIwhere Krnf,2D,i, i∈{I,II,III}, is the stiffness of the segmenti of the 2D reinforcement; see Figure 2(d). For the rectangular segments Krnf,2D,i = Es,effδhi/ci, i∈{I,II}, where hi and ci are height and length of the segmenti of the 2D reinforcement. The stiffness Krnf,2D,III of the trapezium segmentIII(2)Krnf,2D,III=Es,eff∫0cIIIdz/As,effz=δEs,eff∫0cIIIdz/hz=δEs,eff∫0cIIIdz/hs+2λz=δEs,efflog⁡hs+2λc1/2λ-log⁡hs/2λwhereλ=tr/cIII = 1.73/1 = 1.73, tr=1.73 mm is the height of the rib, Figure 2(d). From (1) we obtain that the effective cross-section area of the 2D reinforcement As,eff=4.754·10-6 m2. Then, the effective height of the 2D reinforcement h∗s,eff = As,eff/δ = 13.583 mm. However, duo to short length of the modelled tie the plane section hypothesis does not hold and these analytically obtained quantities are not exact. Using the program made for the modelling of the tie, it was determined that the effective height hs,eff = 12.93 mm and the effective area of the cross-section As,eff=4.525419·10-6 m2 and these values were adopted for the modelling of the 2D model of the tie. On the bases of the equations Es,eff·As,eff = Es·As, and Ec,eff·Ac,eff = Ec·Ac it was obtained that Es,eff = 7143.423 GPa and Ec,eff = 2333.284 GPa, where Ac,eff=δH–h∗s,eff=1.901·10-5 m2, and h∗s,eff = 13.583 mm was obtained using (2).


### 2.4. Discretisation of the 2D Model

The concrete and steel parts of the 2D tie model are discretised by the equilateral triangular elastic spring lattices (ETESL). A spring connects the centres of two vicinal aggregate particles. The springs are pin-connected. The length of the leg of a triangular cell is 0.4 mm; that is, the same as the distance between two vicinal aggregate particles. The principal scheme of the discretisation of the 2D model of the ties is shown in Figure4.Figure 4
The discrete model of the computational 2D model of the tie: a plane view with the boundary conditions (a); a fragment of the model in the vicinity of the rib (b); the one-dimensional spring element combining two particles acting via the interface (c) and (d), and fragment of the lattice (e).The final discrete model consists of 42264 spring elements of equal length, 8734 of them belonging to the 2D reinforcement, and 33530 elements belonging to the concrete.The lattice method is an approximation of the continuum, therefore, has disadvantages. Most important of them are the following: the parameters of the lattice such as the lattice geometry and size, the axial or bending stiffness, the critical tensile strain or stress of the failure criterion, the constitutive relation of the beam or rod have to be determined depending on the continuum, also depending on the chosen criterion of proximity between modelled and the required results of the target continuum. The criterion can be energy stored in the volume of the continuum, displacements or deformations, strength and so on. The next disadvantage is that the Poisson ratio of the continuum approximated by the triangular lattice always equals 1/3. To overcome this problem, the special type of lattices, for example with angular springs, is used; see Ostoja-Starzewski [23]. Another disadvantage is a huge computational time required to perform a large scale analysis. There are also some issues concerning the element removal techniques. Two main techniques “event-by-event” and “event-driven” are used; see Pan et al. [18]. The former technique demands a big computational time while the latter can be inaccurate if the loading increment is too big and too many elements have to be removed; see Pan et al. [18].


### 2.5. Stiffness Properties of the Spring of the Computational Model

According to Hrennikoff [13], the required spring stiffness can be calculated according to the following equation Ks,i = √3/2Eeff,iδ,i∈{c,s}. This formula can be obtained using the well-known formula of the required cross-section area of the spring of the ETESL As = √3/2Lδ given in Hrennikoff [13], where it is assumed that the moduli of elasticity of the spring and the approximated body are the same. The Poison ratio of the body approximated by the ETESL is νi = 1/3,i∈c,s, Hrennikoff [13]. It should be noted that the equation Ks = √3/2Eδ is obtained when the springs are connected rigidly at the nodes. In the present study, the springs are pin-connected at the nodes. However, the given formulas also are valid for the pin-connected ETESL. These formulas were verified numerically using BRIGADE Plus 5.1 program. For the present 2D computational model the following stiffnesses of the springs of ETESL were obtained: Ks,s = 2165.235 MN/m for the steel reinforcement and Ks,c = 707.239 MN/m for the concrete.The fracture of the concrete is modelled by breaking the spring when the strain of the spring exceeds the critical strain, i.e., whenεc>εcrit, where εcrit is the critical tensile strain, and εc is the strain of the spring that represents concrete. According to experimental investigation by Rossello and Elices [31] the cracks always try to avoid the strong aggregates and follows the weaker aggregate-matrix interface or cement matrix. According to Tschegg, Rotter, Roelfstra et al. [32], the first cracks occur in aggregate-matrix interfaces in modeI. Therefore, in the present simulation it is assumed that a spring breaks through the cement paste. The element-size dependent failure criterion elaborated by Tavarez and Plesha [33] reflecting opening fracture mode was applied for the critical strain εcrit. The normal critical tensile strain may be calculated as follows, Tavarez and Plesha [33]:(3)εcrit=πKIc3-νc4Ec6πRp1-νcwhere KIc, νc, and Ec are the fracture toughness, Poisson‘s ratio, and modulus of elasticity of the concrete of the 3D tie model. The fracture toughness of the concrete KIc depends on concrete strength and modulus of elasticity Ec. Usually, KIc increases with increasing concrete strength and modulus of elasticity. In the present research, the following values of the fracture toughness are adopted: KIc∈0.47,0.27,0.067 MPa√m. If Ec=9.954 GPa, then, according to (3) we have εcrit,1 = 2.260·10-3 when KIc=0.47 MPa√m; εcrit,2 = 1.298·10-3 when KIc=0.27 MPa√m; and εcrit,3 = 0.322·10-3 when KIc=0.067 MPa√m. It should be noticed that the triangular lattice of rods without strain softening may lead to too brittle failure of concrete, see Vervuurt, Van Mier and Schlangen [34]. To avoid this shortcoming, the strain softening can be introduced [34–36]. Other approaches to determine the critical tensile strain εcrit are elaborated, for example, dependent on the fracture energy criterion; see Birck et al. [36].


## 2.1. Geometry and Properties of the Reinforced Tie

The geometry of the tie is shown in Figure1. The tie is rectangular with the steel reinforcement placed in the middle of the tie. The initial length of the tie (the initial length of the steel reinforcing rod and concrete part of the tie)L0 = 56.80 mm. The height and the width of the model are of the same size:H=b = 67.9 mm, Figure 1(a). The diameter of the steel reinforcement is adopted as 14 mm.Figure 1
Geometry of the steel: general view of the tie (a), a cross-section of the tie (b).
(a) (b)The left face of the end of the tie is fixed. The displacementu is applied to the free end of the reinforcement from the right, Figure 1(a). The concrete comprises the hardened cement paste and the aggregates assumed as the spherical particles of diameter 0.35 mm embedded in the paste. It is assumed that the moduli of elasticity of the steel reinforcement and concrete are the following: Es=210 GPa, Ec=9.954 GPa. The Poison coefficients are νc=0.2 for concrete and νs=0.33 for the steel reinforcement, respectively.


## 2.2. Computational Model of the Tie

To apply the planar equilateral triangular elastic spring lattice method for the tie, the 3D tie model is reduced to the 2D model consisting of the concrete and the steel reinforcement. The geometry of the 2D model is shown in Figure2. The thickness of the 2D model equals the diameter of the aggregate, i.e., δ = 2Rp=0.35 mm, Figure 2(b). The initial length of the 2D reinforcement and concrete parts are the same as of 3D model L0=56.80 mm and the height of the modelH= 67.9 mm, Figure 2(a). The ribs of the steel reinforcement are trapezoidal with height tr=1.73 mm, and the length cII=6.00 mm, Figure 2(d). The space between the ribs is equal to cI=8.2 mm. It is assumed that the displacements of the left plane of the tie is restrained fully, that is, for the left plane the displacements ux=uy=0.Figure 2
Geometry, boundary conditions, and action of the 2D model of the tie: a scheme of the 2D model (a), a section of the 2D model (b), a geometry of the computational model of the tie (c), and a model of the steel part (d).


## 2.3. Geometry of the Computational Model

For the simulation of the tie, the half of the tie is modelled by using the lattice methodology. The geometry of the model used for the simulation is shown in Figure2(c). The simulation model consists of one part of the concrete and the height of the steel reinforcement equals half of the tie. The thickness of the simulation model, the reinforcement, and the concrete part is the same as for the reduced 2D tie model, δ = 2Rp=0.35 mm. The displacements of the left plane of the simulation model is restrained fully, that is, for the left plane the displacements ux=uy=0 and also the displacements of the bottom plane of the steel reinforcement is restrained in directiony, that is, uy=0, Figure 2(c).It is assumed that the concrete of the 2D model of the tie consists of the spherical aggregate particles of diameter 0.35 mm embedded in the hardened cement paste, Figure3. In the plane, the particles form equilateral triangular cells. The distance of the centres of the two any vicinal particles is 0.4 mm. There is only one layer of the particles in the concrete part of the tie. The volume ratio of the particles (the ratio of the volume of the particles to the volume of the entire concrete) is 0.45815.Figure 3
Structure of the concrete part of the 2D model.The cross-section areas and moduli of elasticity of the 2D model are obtained from the following equations:Es,eff·As,eff = Es·As, and Ec,eff·Ac,eff = Ec·Ac, where As,eff, Ac,eff, Es,eff, and Ec,eff are effective cross-section areas and effective moduli of elasticity of the 2D model of the concrete part and 2D reinforcement, respectively.The effective cross-section areaAs,eff of the 2-reinforcement can be obtained from the condition that the stiffness of the steel part of the 2D model(1)Krnf,2D=As,effEs,effL0=14/Krnf,2D,I+2/Krnf,2D,KII+6/Krnf,2D,IIIwhere Krnf,2D,i, i∈{I,II,III}, is the stiffness of the segmenti of the 2D reinforcement; see Figure 2(d). For the rectangular segments Krnf,2D,i = Es,effδhi/ci, i∈{I,II}, where hi and ci are height and length of the segmenti of the 2D reinforcement. The stiffness Krnf,2D,III of the trapezium segmentIII(2)Krnf,2D,III=Es,eff∫0cIIIdz/As,effz=δEs,eff∫0cIIIdz/hz=δEs,eff∫0cIIIdz/hs+2λz=δEs,efflog⁡hs+2λc1/2λ-log⁡hs/2λwhereλ=tr/cIII = 1.73/1 = 1.73, tr=1.73 mm is the height of the rib, Figure 2(d). From (1) we obtain that the effective cross-section area of the 2D reinforcement As,eff=4.754·10-6 m2. Then, the effective height of the 2D reinforcement h∗s,eff = As,eff/δ = 13.583 mm. However, duo to short length of the modelled tie the plane section hypothesis does not hold and these analytically obtained quantities are not exact. Using the program made for the modelling of the tie, it was determined that the effective height hs,eff = 12.93 mm and the effective area of the cross-section As,eff=4.525419·10-6 m2 and these values were adopted for the modelling of the 2D model of the tie. On the bases of the equations Es,eff·As,eff = Es·As, and Ec,eff·Ac,eff = Ec·Ac it was obtained that Es,eff = 7143.423 GPa and Ec,eff = 2333.284 GPa, where Ac,eff=δH–h∗s,eff=1.901·10-5 m2, and h∗s,eff = 13.583 mm was obtained using (2).


## 2.4. Discretisation of the 2D Model

The concrete and steel parts of the 2D tie model are discretised by the equilateral triangular elastic spring lattices (ETESL). A spring connects the centres of two vicinal aggregate particles. The springs are pin-connected. The length of the leg of a triangular cell is 0.4 mm; that is, the same as the distance between two vicinal aggregate particles. The principal scheme of the discretisation of the 2D model of the ties is shown in Figure4.Figure 4
The discrete model of the computational 2D model of the tie: a plane view with the boundary conditions (a); a fragment of the model in the vicinity of the rib (b); the one-dimensional spring element combining two particles acting via the interface (c) and (d), and fragment of the lattice (e).The final discrete model consists of 42264 spring elements of equal length, 8734 of them belonging to the 2D reinforcement, and 33530 elements belonging to the concrete.The lattice method is an approximation of the continuum, therefore, has disadvantages. Most important of them are the following: the parameters of the lattice such as the lattice geometry and size, the axial or bending stiffness, the critical tensile strain or stress of the failure criterion, the constitutive relation of the beam or rod have to be determined depending on the continuum, also depending on the chosen criterion of proximity between modelled and the required results of the target continuum. The criterion can be energy stored in the volume of the continuum, displacements or deformations, strength and so on. The next disadvantage is that the Poisson ratio of the continuum approximated by the triangular lattice always equals 1/3. To overcome this problem, the special type of lattices, for example with angular springs, is used; see Ostoja-Starzewski [23]. Another disadvantage is a huge computational time required to perform a large scale analysis. There are also some issues concerning the element removal techniques. Two main techniques “event-by-event” and “event-driven” are used; see Pan et al. [18]. The former technique demands a big computational time while the latter can be inaccurate if the loading increment is too big and too many elements have to be removed; see Pan et al. [18].


## 2.5. Stiffness Properties of the Spring of the Computational Model

According to Hrennikoff [13], the required spring stiffness can be calculated according to the following equation Ks,i = √3/2Eeff,iδ,i∈{c,s}. This formula can be obtained using the well-known formula of the required cross-section area of the spring of the ETESL As = √3/2Lδ given in Hrennikoff [13], where it is assumed that the moduli of elasticity of the spring and the approximated body are the same. The Poison ratio of the body approximated by the ETESL is νi = 1/3,i∈c,s, Hrennikoff [13]. It should be noted that the equation Ks = √3/2Eδ is obtained when the springs are connected rigidly at the nodes. In the present study, the springs are pin-connected at the nodes. However, the given formulas also are valid for the pin-connected ETESL. These formulas were verified numerically using BRIGADE Plus 5.1 program. For the present 2D computational model the following stiffnesses of the springs of ETESL were obtained: Ks,s = 2165.235 MN/m for the steel reinforcement and Ks,c = 707.239 MN/m for the concrete.The fracture of the concrete is modelled by breaking the spring when the strain of the spring exceeds the critical strain, i.e., whenεc>εcrit, where εcrit is the critical tensile strain, and εc is the strain of the spring that represents concrete. According to experimental investigation by Rossello and Elices [31] the cracks always try to avoid the strong aggregates and follows the weaker aggregate-matrix interface or cement matrix. According to Tschegg, Rotter, Roelfstra et al. [32], the first cracks occur in aggregate-matrix interfaces in modeI. Therefore, in the present simulation it is assumed that a spring breaks through the cement paste. The element-size dependent failure criterion elaborated by Tavarez and Plesha [33] reflecting opening fracture mode was applied for the critical strain εcrit. The normal critical tensile strain may be calculated as follows, Tavarez and Plesha [33]:(3)εcrit=πKIc3-νc4Ec6πRp1-νcwhere KIc, νc, and Ec are the fracture toughness, Poisson‘s ratio, and modulus of elasticity of the concrete of the 3D tie model. The fracture toughness of the concrete KIc depends on concrete strength and modulus of elasticity Ec. Usually, KIc increases with increasing concrete strength and modulus of elasticity. In the present research, the following values of the fracture toughness are adopted: KIc∈0.47,0.27,0.067 MPa√m. If Ec=9.954 GPa, then, according to (3) we have εcrit,1 = 2.260·10-3 when KIc=0.47 MPa√m; εcrit,2 = 1.298·10-3 when KIc=0.27 MPa√m; and εcrit,3 = 0.322·10-3 when KIc=0.067 MPa√m. It should be noticed that the triangular lattice of rods without strain softening may lead to too brittle failure of concrete, see Vervuurt, Van Mier and Schlangen [34]. To avoid this shortcoming, the strain softening can be introduced [34–36]. Other approaches to determine the critical tensile strain εcrit are elaborated, for example, dependent on the fracture energy criterion; see Birck et al. [36].


## 3. Investigated Models

Six kinds of the tie of the same properties except for the critical concrete strainsεcrit were modelled to investigate the force tensile behaviour of the tie. That is, the geometry of the 2D computation model, Figures 2 and 3, the discretisation and lattice, Figure 4, and the stiffness properties of the lattice of the concrete and steel parts of the computational 2D model are the same for all six specimens. Only the critical tensile strains εcrit differ for each examined model. Three deterministic models, DM1, DM2, and DM3, were modelled with constant, i.e., deterministic, values of the critical tensile strains of the concrete εcrit∈2.27·10-3,1.298·10-3,0.324·10-3 and three random models, RM1, RM2, and RM3, were modelled with spatially variable properties of the critical strains εcrit of the concrete.For the stochastic modelling of the cracking of the concrete and the tie behaviour the critical concrete tensile strainsεcrit are modelled as two-dimensional stationary uncorrelated truncated Gaussian random field Xϵ,crit,ix:x∈S⊂R2, wherex is coordinate vector of a spring of the concrete part of the 2D tie model,S is the region of the concrete part of the 2D tie model, and Xϵ,crit,ix is the truncated normal random variable. The mean E(Xϵ,crit,ix)=με,crit,i, x∈S and covariation cov⁡(Xϵ,crit,ix1,Xϵ,crit,ix2)=0 for all x1,x2∈S,x1≠x2. In general, the concrete properties are spatially correlated. However, when the conditions of the casting of a tie are constant, the fresh concrete mass is homogeneous and volume of the tie is not large; it is reasonable to assume that the spatial variability of the properties of the concrete is uncorrelated and only random deviations of the concrete properties can be taken into account.The random variableXϵ,crit,ix (hereafter Xε,crit,i) is the truncated normal r.v. attaining values within interval [με,crit,i–3sdε,crit,i,με,crit,i+3sdε,crit,i] = [με,crit,i(1-3Cvar(Yε,crit,i)), με,crit,i(1+3Cvar(Yε,crit,i))], where με,crit,i, sdε,crit,i, and Cvar(Yε,crit,i) are mean, standard deviation and coefficient of variation of the untruncated r.v. Yε,crit,i ~ N(με,crit,i,sdε,crit,i): M(Yε,crit,i) = M(Xε,crit,i) = με,crit,i = εcrit,i, √D(Yε,crit,i) = sdε,crit,i andCvar(Yε,crit,i) =sdε,crit,i /με,crit,i, i∈1,2,3. Three cases of the behaviour of the tie were considered when the critical tensile strain truncated normal r.v. Xε,crit,i∈[με,crit,i-3sdε,crit,i,με,crit,i+3sdε,crit,i], i∈1,2,3, where με,crit,1 = εcrit,1 = 2.260·10-3, sdε,crit,1 = 4.080·10-4; με,crit,2 = εcrit,2 = 1.298·10-3, sdε,crit,2 = 2.336·10-4, and με,crit,3 = εcrit,3 = 0.322·10-3, sdε,crit,3 = 5.796·10-5; Cvar(Yε,crit,i) = 0.18 for all r.v. Yε,crit,i. The standard deviations √D(Xε,crit,i) of the truncated r.v. Xε,crit,i are the following: √D(Xε,crit,1) = 4.0128·10-4, √D(Xε,crit,2) = 2.305·10-4, and √DXε,crit,3 = 5.717·10-5. The coefficient of variation of the truncated normal random variables of the critical tensile concrete strain is Cvar(Xε,crit,i) = 0.1776, i∈1,2,3.For the sake of the convenience, the elastic constants of the reinforcement, concrete, the effective elastic constants of the 2D model, and the stiffness of the springs of the 2D model are summarised in Table1.Table 1
Summary of properties of the deterministic and random models of the computational models of the ties.
Designation of the specimen Moduli of elasticity for the 3D model in GPa Effective moduli of elasticity for the 2D model in GPa Stiffness of the springKs for the lattice model in MN/m Critical tensile strainsεcrit of the concrete E s for the steel reinforcement E c for of the concrete E s , e f f for 2D reinforcement E c , e f f for concrete K s , s for the 2D reinforcement K s , c for the concrete Mean valuesμε,crit Coefficient of variationCvar (Yε,crit) 3D tie model 2D model and computational model  Deterministic models  DM1 210 9.954 7143.423 2333.284 2165.235 707.239 2.260·10−3 0 DM2 1.298·10−3 DM3 0.322·10−3  Spatially variable model  RM1 210 9.954 7143.423 2333.284 2165.235 707.239 2.260·10−3 0.18 RM2 1.298·10−3 RM3 0.322·10−3


## 4. Simulation Methodology

The discrete element methodology, Cundall [37], has been applied to solve the problem of the deformation of the tie. The coordinates and velocities of each particle are determined by the integrating the dynamic equation mid2xi(t)/dt2 = Fi(t), wherei∈Ip = 1,...,N, Ip is an indexing set of the nodes,N is number of the nodes, mi is node mass, mi=2.54·10-4 kg for the concrete nodes and mi=9.063·10-4 kg for the steel nodes, and xi=xi,yiT and Fi=Fx,i,Fy,iT are the vectors of the coordinates and force respectively. The total force Fi=∑j∈Ip∖{i}Fij, where Fij=nijKs,αΔlij – ηvij(t) is the vector of the force acting between nodesi andj; Ks,α, α∈s,c is the stiffness of a spring representing steel, in case of Ks,s, or concrete, in case of Ks,c, Ks,α is given in Table 1; Δlij = √(lij,x2 + lij,y2) is the elongation or shortening of the spring connecting nodesi andj; nij = xij/xij is the unit vector, where xij=xij,yijT is vector connecting nodesi andj, and xij = √(xij2 + yij2) is the Euclidean norm or the length of the vector,η = 2c√(miKs,i) is the damping coefficient whose the damping ratioc = 0.005, vij(t) = dxij(t)/dt is the relative velocity vector. The damping has been introduced so as to make Euler’s explicit integration scheme more stable: vi(tn+1) = vi(tn) + Fi(tn)/miΔt, and xi(tn+1) = xi(tn) + vi(tn)Δt, where Δt=8·10-8 s is the time step. The external displacementu has been imposed with the constant velocitydu(t)/dt = 1 mm/s in the directionx to the free end of the reinforcement.The simulation program has been written in C++11. The geometry of the model and nodes coordinates for the computation program have been created by the finite element program ANSYS 12. One realization of the stochastic process of the force-displacement behaviour of a tie took at least 7 hours using the one core of the computer whose parameters were as follows: RAM 4GB, CPU AMD Phenom II X4 955 3.2 GHz 4 cores.


## 5. Simulation Results and Discussions

In this section, the results of the stochastic analysis of the tensile force and the number of the broken springs of the concrete of the tie dependent on the displacement of the reinforcement are presented. The statistical analysis was performed with the program “RStudio” [38], and using the following R packages: “stats” [38], “corrplot” [39], and “e1071” [40]. The coefficients of variations, standard deviations, range, autocorrelation or coefficients of correlation as well as the hypotheses about the normality of the distribution of the force r.v. of the 72 realizations dependent on the reinforcement displacement of the model RM2 were considered. Also, the crack patterns of the deterministic DMi and random RMi, i∈1,2,…,3, specimens are considered.


### 5.1. Force-Elongation Behaviour of the Steel Reinforcement

The influence of the random critical concrete tensile strainεcrit on the force-displacement and fracture behaviour of the tie is presented in the current section. Also, a comparison of the force-displacement and fracture behaviour of the ties with the stochastic and deterministic critical concrete tensile strains εcrit is presented in this section. 72 realizations of the specimen RM2 and one realization for each specimen RM1 and RM3 were calculated.The dependence of the tensile axial forcesFDM2 and estimated average force F-RM2 of the models DM2 and RM2 on the displacement u∈0,0.1 mm of the free end of the steel reinforcement is depicted in Figure 5. The estimated average force F-RM2u=1/n∑i=1nFRM2,iu,n = 72, where FRM2,i is realizationi of the axial forceF of the specimen RM2. In this figure, the supremum Fsup is calculated by assuming that εcrit=∞ and the infimum Finf (u) = (AsˑEs)ˑu/L0 = (As,efffˑEs,eff)ˑu/L0.Figure 5
Tensile axial force-displacement relationships of the tie: deterministicFDM2, estimated average force F-RM2 of the model RM2 as well as infimum Finf and supremum Fsup.As we can see from Figure5, the estimated average tensile axial force F-RM2 is smooth in comparison to the FDM2. Also, from Figure 5 we can see that FDM2 and F-RM2 approach to Finf as the displacementu increases and FDM2 and F-RM2 approach to Fsup as the displacementu decreases. It is clear that Finf<F<Fsup, where F∈FDM2,F-RM2. The FDM2 and F-RM2 forces almost coincide when the displacementu > 0.016 mm andu < 0.07 mm. However, when the displacementu∈(0.016, 0.07), the force FDM2 may differ from F-RM2 significantly.The dependencies of the axial tensile forceFRM2 and the number of broken springs Nbr,RM2 of the three different realizations of the specimen RM2 on the displacementu are shown in Figure 6. Also, in this figure the following dependencies on the displacementu are shown: the minimum FRM2,min(u) = min⁡FRM2,iu,i∈IRM2 and the maximum FRM2,max(u) = max⁡{FRM2,i(u),i∈IRM2} of the axial tensile force, as well as the minimum Nbr,RM2,min(u) = min⁡Nbr,RM2,iu,i∈IRM2 and the maximum Nbr,RM2,max(u) = max⁡{Nbr,RM2,i(u) of the number of the broken springs, where FRM2,i and Nbr,RM2,i are realizationi of the r.vs. XF,RM2 and XN,br,RM2 respectively, i∈IRM2, IRM2 = {1,2,…,72} is an indexing set of the realizations of the specimen RM2. In Figure 6, the vertical dotted lines denote the displacementsu at which the estimated standard deviation and coefficient of variation of r.v. XN,br,RM2 attain local maximums.Figure 6
Dependencies of the tensile forceFRM2  (a) and number of broken springs Nbr,RM2  (b) on the displacement of the steel reinforcementu of model RM2 obtained by the stochastic analysis when the critical tensile strain mean με,crit,2 = 1.298·10-3 and the standard deviation √D(Xε,crit,2) = 2.305·10-4.Figure6 shows that the variation of the force FRM2 and the number of broken springs Nbr,RM2 is small at the small displacementsu. When the spring, which represents concrete tensile strain εc, is very small, then the probabilityPr(Xε,crit ≤ εc) of the breaking of the spring is also very small. Therefore, the number Nbr,RM2 of broken springs is also very small. The tie is almost monolithic. The density of the probabilityPr(Xε,crit ≤ εc) at the small displacementu is also small. As a consequence, the tensile force FRM2 relationship on the displacementu practically does not depend on the critical concrete tensile strains. Therefore, at the beginning of the deformation, we have that the relationship between the tensile forceF and the displacementu is almost deterministic rather that Finf(u) ≤F(u) ≤ Fsup(u) and Finf(u) → Fsup(u) asu → 0. With increasing the displacementu, the intensity of the breaking of the springs increases and the variability of the tensile forceF increases too. It can be clearly seen from Figures 6(a) and 6(b) whenu∈ [0.03, 0.07] mm. In this interval, the increase of the number of the broken springs is the biggest in comparison to the other intervalsu. When the displacement increases, i.e., whenu > 0.06 mm, the probabilityPr(Xε,crit ≤ εc) of the spring breaking is big; however, the density of the probabilityPr(Xε,crit ≤ εc) is decreased and is small. In addition, the tensile forceF(u) approaches to the Finf(u) asu increases since the damaged bond between the reinforcement and the concrete decreases the stiffness of the tie. These reasons lead to the decreased variability of the tensile forceF with increasing the displacementu. This statement is clearly illustrated in Figures 7, 8, 9, and 10.Figure 7
Coefficient of variationCvarF,RM2 of the force r.v. XF,RM2.Figure 8
Estimated coefficient of variationCvarN,br,RM2 of the broken spring number r.v. XN,br,RM2.Figure 9
The differencesFRM2,max - FRM2,min and the estimated standard deviation sF,RM2 of the force r.v. XF,RM2.Figure 10
The differencesNbr,RM2,max – Nbr,RM2,min and the estimated standard deviation sN,br,RM2 of the number of broken springs r.v. XN,br,RM2.The dependencies of the estimated coefficients of variationCvarF,RM2 = sF,RM2/mF,RM2 and CvarN,br,RM2 = sN,br,RM2/mN,br,RM2 of the force r.v. XF,RM2 and the number of the broken spring r.v. XN,br, respectively, are shown in Figures 7 and 8, where sF,RM2, sN,br,RM2, mF,RM2, and mN,br,RM2 are estimations of the standard deviations and the means of the force r.v. XF,RM2 and the number of the broken spring r.v. XN,br,RM2, respectively, mj,RM2=1/n∑i=1njRM2,i, s2j,RM2=1/n-1∑i=1njRM2,i-mj,RM22, j∈F,Nbr,n = 72. In Figures 6, 7, and 8, the first four vertical dotted lines show the displacementsu∈0.016397, 0.035185, 0.04509, 0.060122 mm at which CvarF,RM2 attains local maximums.As we can see, from Figures7 and 8, the maximums of CvarF,RM2 are very close to the maximums of the estimated coefficients of variations CvarN,br,RM2 of the number of the broken springs. From Figures 7 and 8 we also can see that CvarF,RM2 in general decreases while CvarN,br,RM2 is almost constant whenu∈[0.07, 0.35] mm. The same statement is valid for the estimation of the standard deviation sF,RM2 of r.v. XF,RM2, Figure 9.The dependencies of the differencesFRM2,max – FRM2,min and Nbr,RM2,max – Nbr,RM2,min of the axial tensile force FRM2 and the number of broken springs Nbr,RM2 of the specimen RM2 on the displacementu are shown in Figures 9 and 10. From Figure 10 we can see that, in general, sN,brRM2 increases whenu∈[0, 0.35] mm, while the difference Nbr,RM2,max – Nbr,RM2,min does not have the unambiguous tendency. However, we cannot say that the difference Nbr,RM2,max – Nbr,RM2,min decreases with increasing the displacementu.The histograms and the kernel density estimations of the realizationsFRM2,i(u),i∈IRM2, of the force r.v. XF,RM2(u), atu∈{0.016397,0.035185,0.04509,0.060122} mm at which the estimations CvarF,RM2 and sN,br,RM2 attain local maximums, are shown in Figure 11. As we can see, the histograms and the estimations of the densities are not symmetrical with respect to the mean. The coefficient of skewness varies between -0.0172 and 0.4306. Also, the histograms and the kernel density estimations show that r.v. XF,RM2 is neither a unimodal nor a Gaussian. For example, thepvalues,p, of the Shapiro-Wilk test for the hypothesis on the normality of the models of realizations of the r.v. XF,RM2(u) are the following: p=5.98·10-7 foru = 0.016397 mm, p=2.77·10-9 foru = 0.03518 mm, p=8.6·10-3 foru = 0.04509 mm, andp = 0.131 foru = 0.060122 mm. Only whenu = 0.060122 the hypothesis about the normality of r.v. XF,RM2 can be accepted. For the first three cases, Figures 11(a), 11(b), and 11(c), thep << α = 0.05, where α is significance level. Only for the fourth case whenu = 0.06012, Figure 11(d),p > α = 0.05 and the hypothesis on normality of r.v. XF,RM2 can be accepted.Figure 11
Histograms of the tensile forceFRM2 and kernel density estimations of 72 realizations of the force random variable XF,RM2 whenu∈ {0.016397,0.035185,0.04509,0.060122} mm.
(a) (b) (c) (d)The histograms and kernel density estimations of the realizationsFRM2,i(u),i∈IRM2, of the force r.v. XF,RM2(u), at u∈0.10,0.1503,0.20,0.25 mm, are shown in Figure 12. As we can see from these figures, the histograms and the estimations of the densities also are not symmetrical with respect to the estimated means. The skewness varies between -0.339 and 1.180. In the present case r.v. XF,RM2 can be treated as a unimodal and Gaussian. Thepvalues of the Shapiro-Wilk test for the hypothesis on the normality of the models of the realizations of the r.v. XF,RM2(u) are the following:p = 0.953 foru = 0.10 mm, p=4.795·10-5 foru = 0.1503 mm,p = 0.5378 foru = 0.200 mm,p = 0.153 foru = 0.250 mm. Only whenu = 0.1503 mm, Figure 12(b), the hypothesis of the normality of r.v. XF,RM2 must be rejected at the significance level α = 0.05; for the rest cases the hypotheses about normality can be accepted; see Figures 12(a), 12(c), and 12(d).Figure 12
Histograms of the tensile forceFRM2 and kernel density estimations of 72 realizations of the force random variable XF,RM2 when u∈{0.1,01503,020,0.250} mm.
(a) (b) (c) (d)The analysis of the normality of the realizations of the r.v.XF,RM2(u), where u∈{0.01+0.0025k:k∈{0,…,36}}={0.01,…,0.1} mm, showedpvalue of the Shapiro-Wilk test,p < α = 0.05, whenu ≤ 0.05 mm. Therefore, it is a reason to claim that r.v. XF,RM2(u) is not a Gaussian whenu ≤ 0.05 mm. Whenu > 0.05 mm, then in majority of cases, at 15 points of the remaining 20,p ≥ 0.05. Therefore, whenu > 0.05 mm, there is a reason to treat r.v. XF,RM2(u) as a Gaussian r.v. Also, on the basis of the performed analysis it is reasonable to claim that the stochastic process {XF,RM2(u):u∈[0,0.05]mm} is not Gaussian, while the stochastic process {XF,RM2(u):u∈(0.05,0.35]mm} may be treated as a Gaussian. Combining these conclusions and Figure 6(b) we can also state that the stochastic process {XF,RM2(u):u∈[0,0.05]mm} is not a Gaussian within the interval in which the new concrete cracks form most intensive. Whenu > 0.05 mm, then the intensity of the formation of the new concrete cracks is less and the stochastic process {XF,RM2(u):u∈(0.05,0.35]mm} may be treated as a Gaussian. The decreasing intensity of the formation of the new concrete cracks whenu > 0.05 mm can be seen very clearly from realization 2 (dotted line) of r.v. XN,br,RM2 shown on Figure 6(b). It should be noted that the considered realizations of r.vs. XF,RM2 and XN,br,RM2 are obtained by assuming that the critical concrete tensile strain r.vs. Xε,crit,i, j∈{1,2,3} are distributed according to truncated normal variable when the coefficient of the variationCvar(Xε,crit,i) = 0.1776. For the other coefficients of variation, the intervals in which the stochastic processes can be treated as a non-Gaussian or Gaussian could be different. However, the performed analysis allows us to claim that r.v. XF,RM2 is not Gaussian within the intervals in which the intensity of crack formation is very big. And on the contrary, r.v. XF,RM2 becomes Gaussian with decreasing the intensity of the formation of the new concrete cracks. It should be noticed that in the common reinforced structures, there are at least three characteristic intervals of the displacements of the cracking of the reinforcement of ties; see Model Code 2010 [41]: the first uncracked interval starts at the zero and ends when the first crack appears. Then the next, second cracking interval starts in which the new cracks appear very intensively. Usually, the force-displacement curve of this interval is uniform like yield plateau of mild steel. The intensity of the formation of the cracks in this displacement interval is of the most intensity. The next, third interval starts when the intensity of the formation of the new cracks becomes lesser. In this interval, the force-displacement curve increases with increasing the displacement of the reinforcement. The obtained force-displacement relationships of the 2D model also have similar intervals, and the relationships are not Gaussian at the big crack formation intensity. Therefore, the conclusion that the force-displacement relationship is not Gaussian at least within the second interval of the cracking can be extended in general also for the common reinforced concrete ties. However, this question demands more deep investigations.The correlation matrices of the correlation coefficientsρi,j = ρ(XF,RM2(ui), XF,RM2(uj)) of the tensile force r.v. XF,RM2 are shown in Figure 13. Figure 13(a) shows ρi,j when ui,uj∈{0.01+0.0025k:k∈0,…,36} and Figure 14(b) when ui,uj∈{0.1+0.01k:k∈{0,…,20}}. ρi,j attains the relatively low values at smallu∈[0.01,0.1] mm in comparison to the ρi,j at the bigger displacementsu∈[0.1, 0.3] mm. For example, ρ(XF,RM2(0.06), XF,RM2(0.1)) = 0.209 while ρ(XF,RM2(0.26), XF,RM2(0.3)) = 0.662. The coefficient of correlation ρ(XF,RM2(ui), XF,RM2(uj)) increases with the decreased distance abs(ui - uj).Figure 13
The correlation matrices of r.v.XF,RM2: when u∈{0.01+0.0025k:k∈{0,…,36}} (a) and u∈{0.1+0.01k:k∈{0,…,20}} (b).
(a) (b)Figure 14
Dependencies of the axial forceF (a) and number of broken springs Nbr (b) on the displacementu.
(a) (b)The weak correlationsρ(XF,RM2(ui), XF,RM2(uj)) shows weak linear dependence between tensile forces FRM2(ui) and XF,RM2(uj). This linear dependence decreases with increasing the distance abs(ui - uj). In the majority of cases, ρi,j > 0. It shows that the increase tensile force FRM2(ui) leads to the increased FRM2(uj), where ui < uj. However, this tendency decreases with the increased distance abs(ui - uj). Despite the fact that the correlation coefficients ρ(XF,RM2(ui), XF,RM2(uj)) decrease with the increased distance abs(ui - uj), we cannot claim that r.v. XF,RM2(uj) becomes less dependent on r.v. XF,RM2(ui) with increasing abs(ui - uj), where uj > ui since Finf (u) ≤F(u) ≤ Fsup (u) and Finf(u) as well as Fsup(u) increase with increasing the displacementu. It should be noticed that the negative coefficients of correlation ρ(XF,RM2(uξ), XF,RM2(uη)) < 0 exist at the displacements (uξ,uη) at which the increasing tensile forceF(uξ) leads to the decreasedF(uη). This decrease is caused by the formation of the new cracks in concrete whenu∈ [uξ,uη].The dependencies of the tensile forcesF∈FRM1,FRD1,FRM2,FDM2,FRM3,FDM3 and the number of the broken springs Nbr∈Nbr,RM1,Nbr,RD1, Nbr,RM2, Nbr,DM2, Nbr,RM3, Nbr,DM3 on the displacement u are depicted in Figure 14. In this figure FRMi and Nbr,RMi, i∈{1,2,3}, denote single realizations of r.vs. XF,RMi and XN,br,RMi. As we can see, Nbr,RMi(u)≥Nbr,DMi(u), where i∈{1,2,3}, u∈[0,0.2] mm, in major part of the interval [0,0.2] mm. However there is no unique tendency on FDMi(u)≥FRMi(u). It is evident that FRM1 ≥ FRM2 ≥ FRM3, FDM1 ≥ FDM2 ≥ FDM3, Nbr,RM1 ≤ Nbr,RM2 ≤ Nbr,RM3, and Nbr,DM1 ≤ Nbr,DM2 ≤ Nbr,DM3 since με,crit,1 > με,crit,2>με,crit,3.Figure14 shows that whenu > 0.162 mm, then Nbr,DM2 < Nbr,RM1 even though με,crit,2>με,crit,1. From Figure 14 we can also see that with decreasing critical tensile strain εcrit the tensile forceF approaches to the Finf. The intense cracks formation stage starts at the bigger displacements of the reinforcement for the specimens with bigger critical concrete tensile strains εcrit, and these displacements are bigger for the deterministic specimens than for random.


### 5.2. Cracking of the Tie

The evaluation of the cracks pattern of the specimens DM2 and DM3, when the critical concrete tensile strainεcrit,2 = 1.298·10-3 and εcrit,3 = 0.322·10-3, at the displacementsu∈(0.04, 0.072, 0.104, 4.0) mm, is shown in Figures 15 and 16. As we can see from these pictures, whenu = 0.04 mm, the debonding crack (the crack between the steel reinforcement and the concrete) is longer for the model DM3 than for DM2. For this case, the debonding crack reaches the second ridge from the left for the model DM2 and the first ridge from the left for the model DM3. Also, we can see that the cracks are higher and more inclined for the model DM2 than for DM3.Figure 15
Cracks patterns of the specimen DM2, whenεcrit,2 = 1.298·10-3, at the displacementsu∈ (0.04, 0.072, 0.104, 4.0) mm.
(a) (b) (c) (d)Figure 16
Cracks patterns of the specimen DM3, whenεcrit,3 = 0.322·10-3, at the displacementsu∈ (0.04, 0.072, 0.104, 4.0) mm.
(a) (b) (c) (d)Cracking patterns of the deterministic, DM1, DM2 and DM3, and random, RM1, RM2 and RM3, models at the displacementu = 0.2 mm are shown in Figure 17.Figure 17
Cracking patterns of the specimens DMi and RMi,i∈ {1,2,3} atu = 0.2 mm: (a), (c), and (e) denote the deterministic models RD1, RD2, and RD3, respectively; (b), (d), and (f) denote the random models RM1, RM2, and RM3, respectively.
(a) (b) (c) (d) (e) (f)From Figure17 we can notice that the spatial variability of the cracks of all single realizations of the random specimens RM1, RM2, and RM3 is more significant in comparison to the deterministic specimens DM1, DM2, and DM3. Also, we can see that the spatial variability of the cracks is more significant for a realization of the specimen RM3 than for RM2 and RM1, Figures 17(f), 17(d), and 17(b) respectively. Therefore, the decreasing critical tensile strains also increase the cracks variability of the concrete. It should be noted that in principle the location and inclination of the modelled concrete cracks, Figures 15, 16, and 17, near the ribs of the reinforcement, are very similar to the cracks obtained experimentally; see Goto and Otsuka [42]. That is the cracks are inclined in the side of the imposed reinforcement force. Also, the angle of the inclination and location of the modelled cracks near the ribs were compared with the numerical results obtained by modelling the identical yet homogeneous and unraced tie using the finite element method. The comparison showed that the location of the cracks near the ribs obtained by the lattice method corresponds to the location of the maximum principal tensile stresses in the ribs region obtained by the finite element method. Moreover, the inclination of the principal planes of the maximum principal tensile stresses is similar to the inclination of the modelled cracks.


## 5.1. Force-Elongation Behaviour of the Steel Reinforcement

The influence of the random critical concrete tensile strainεcrit on the force-displacement and fracture behaviour of the tie is presented in the current section. Also, a comparison of the force-displacement and fracture behaviour of the ties with the stochastic and deterministic critical concrete tensile strains εcrit is presented in this section. 72 realizations of the specimen RM2 and one realization for each specimen RM1 and RM3 were calculated.The dependence of the tensile axial forcesFDM2 and estimated average force F-RM2 of the models DM2 and RM2 on the displacement u∈0,0.1 mm of the free end of the steel reinforcement is depicted in Figure 5. The estimated average force F-RM2u=1/n∑i=1nFRM2,iu,n = 72, where FRM2,i is realizationi of the axial forceF of the specimen RM2. In this figure, the supremum Fsup is calculated by assuming that εcrit=∞ and the infimum Finf (u) = (AsˑEs)ˑu/L0 = (As,efffˑEs,eff)ˑu/L0.Figure 5
Tensile axial force-displacement relationships of the tie: deterministicFDM2, estimated average force F-RM2 of the model RM2 as well as infimum Finf and supremum Fsup.As we can see from Figure5, the estimated average tensile axial force F-RM2 is smooth in comparison to the FDM2. Also, from Figure 5 we can see that FDM2 and F-RM2 approach to Finf as the displacementu increases and FDM2 and F-RM2 approach to Fsup as the displacementu decreases. It is clear that Finf<F<Fsup, where F∈FDM2,F-RM2. The FDM2 and F-RM2 forces almost coincide when the displacementu > 0.016 mm andu < 0.07 mm. However, when the displacementu∈(0.016, 0.07), the force FDM2 may differ from F-RM2 significantly.The dependencies of the axial tensile forceFRM2 and the number of broken springs Nbr,RM2 of the three different realizations of the specimen RM2 on the displacementu are shown in Figure 6. Also, in this figure the following dependencies on the displacementu are shown: the minimum FRM2,min(u) = min⁡FRM2,iu,i∈IRM2 and the maximum FRM2,max(u) = max⁡{FRM2,i(u),i∈IRM2} of the axial tensile force, as well as the minimum Nbr,RM2,min(u) = min⁡Nbr,RM2,iu,i∈IRM2 and the maximum Nbr,RM2,max(u) = max⁡{Nbr,RM2,i(u) of the number of the broken springs, where FRM2,i and Nbr,RM2,i are realizationi of the r.vs. XF,RM2 and XN,br,RM2 respectively, i∈IRM2, IRM2 = {1,2,…,72} is an indexing set of the realizations of the specimen RM2. In Figure 6, the vertical dotted lines denote the displacementsu at which the estimated standard deviation and coefficient of variation of r.v. XN,br,RM2 attain local maximums.Figure 6
Dependencies of the tensile forceFRM2  (a) and number of broken springs Nbr,RM2  (b) on the displacement of the steel reinforcementu of model RM2 obtained by the stochastic analysis when the critical tensile strain mean με,crit,2 = 1.298·10-3 and the standard deviation √D(Xε,crit,2) = 2.305·10-4.Figure6 shows that the variation of the force FRM2 and the number of broken springs Nbr,RM2 is small at the small displacementsu. When the spring, which represents concrete tensile strain εc, is very small, then the probabilityPr(Xε,crit ≤ εc) of the breaking of the spring is also very small. Therefore, the number Nbr,RM2 of broken springs is also very small. The tie is almost monolithic. The density of the probabilityPr(Xε,crit ≤ εc) at the small displacementu is also small. As a consequence, the tensile force FRM2 relationship on the displacementu practically does not depend on the critical concrete tensile strains. Therefore, at the beginning of the deformation, we have that the relationship between the tensile forceF and the displacementu is almost deterministic rather that Finf(u) ≤F(u) ≤ Fsup(u) and Finf(u) → Fsup(u) asu → 0. With increasing the displacementu, the intensity of the breaking of the springs increases and the variability of the tensile forceF increases too. It can be clearly seen from Figures 6(a) and 6(b) whenu∈ [0.03, 0.07] mm. In this interval, the increase of the number of the broken springs is the biggest in comparison to the other intervalsu. When the displacement increases, i.e., whenu > 0.06 mm, the probabilityPr(Xε,crit ≤ εc) of the spring breaking is big; however, the density of the probabilityPr(Xε,crit ≤ εc) is decreased and is small. In addition, the tensile forceF(u) approaches to the Finf(u) asu increases since the damaged bond between the reinforcement and the concrete decreases the stiffness of the tie. These reasons lead to the decreased variability of the tensile forceF with increasing the displacementu. This statement is clearly illustrated in Figures 7, 8, 9, and 10.Figure 7
Coefficient of variationCvarF,RM2 of the force r.v. XF,RM2.Figure 8
Estimated coefficient of variationCvarN,br,RM2 of the broken spring number r.v. XN,br,RM2.Figure 9
The differencesFRM2,max - FRM2,min and the estimated standard deviation sF,RM2 of the force r.v. XF,RM2.Figure 10
The differencesNbr,RM2,max – Nbr,RM2,min and the estimated standard deviation sN,br,RM2 of the number of broken springs r.v. XN,br,RM2.The dependencies of the estimated coefficients of variationCvarF,RM2 = sF,RM2/mF,RM2 and CvarN,br,RM2 = sN,br,RM2/mN,br,RM2 of the force r.v. XF,RM2 and the number of the broken spring r.v. XN,br, respectively, are shown in Figures 7 and 8, where sF,RM2, sN,br,RM2, mF,RM2, and mN,br,RM2 are estimations of the standard deviations and the means of the force r.v. XF,RM2 and the number of the broken spring r.v. XN,br,RM2, respectively, mj,RM2=1/n∑i=1njRM2,i, s2j,RM2=1/n-1∑i=1njRM2,i-mj,RM22, j∈F,Nbr,n = 72. In Figures 6, 7, and 8, the first four vertical dotted lines show the displacementsu∈0.016397, 0.035185, 0.04509, 0.060122 mm at which CvarF,RM2 attains local maximums.As we can see, from Figures7 and 8, the maximums of CvarF,RM2 are very close to the maximums of the estimated coefficients of variations CvarN,br,RM2 of the number of the broken springs. From Figures 7 and 8 we also can see that CvarF,RM2 in general decreases while CvarN,br,RM2 is almost constant whenu∈[0.07, 0.35] mm. The same statement is valid for the estimation of the standard deviation sF,RM2 of r.v. XF,RM2, Figure 9.The dependencies of the differencesFRM2,max – FRM2,min and Nbr,RM2,max – Nbr,RM2,min of the axial tensile force FRM2 and the number of broken springs Nbr,RM2 of the specimen RM2 on the displacementu are shown in Figures 9 and 10. From Figure 10 we can see that, in general, sN,brRM2 increases whenu∈[0, 0.35] mm, while the difference Nbr,RM2,max – Nbr,RM2,min does not have the unambiguous tendency. However, we cannot say that the difference Nbr,RM2,max – Nbr,RM2,min decreases with increasing the displacementu.The histograms and the kernel density estimations of the realizationsFRM2,i(u),i∈IRM2, of the force r.v. XF,RM2(u), atu∈{0.016397,0.035185,0.04509,0.060122} mm at which the estimations CvarF,RM2 and sN,br,RM2 attain local maximums, are shown in Figure 11. As we can see, the histograms and the estimations of the densities are not symmetrical with respect to the mean. The coefficient of skewness varies between -0.0172 and 0.4306. Also, the histograms and the kernel density estimations show that r.v. XF,RM2 is neither a unimodal nor a Gaussian. For example, thepvalues,p, of the Shapiro-Wilk test for the hypothesis on the normality of the models of realizations of the r.v. XF,RM2(u) are the following: p=5.98·10-7 foru = 0.016397 mm, p=2.77·10-9 foru = 0.03518 mm, p=8.6·10-3 foru = 0.04509 mm, andp = 0.131 foru = 0.060122 mm. Only whenu = 0.060122 the hypothesis about the normality of r.v. XF,RM2 can be accepted. For the first three cases, Figures 11(a), 11(b), and 11(c), thep << α = 0.05, where α is significance level. Only for the fourth case whenu = 0.06012, Figure 11(d),p > α = 0.05 and the hypothesis on normality of r.v. XF,RM2 can be accepted.Figure 11
Histograms of the tensile forceFRM2 and kernel density estimations of 72 realizations of the force random variable XF,RM2 whenu∈ {0.016397,0.035185,0.04509,0.060122} mm.
(a) (b) (c) (d)The histograms and kernel density estimations of the realizationsFRM2,i(u),i∈IRM2, of the force r.v. XF,RM2(u), at u∈0.10,0.1503,0.20,0.25 mm, are shown in Figure 12. As we can see from these figures, the histograms and the estimations of the densities also are not symmetrical with respect to the estimated means. The skewness varies between -0.339 and 1.180. In the present case r.v. XF,RM2 can be treated as a unimodal and Gaussian. Thepvalues of the Shapiro-Wilk test for the hypothesis on the normality of the models of the realizations of the r.v. XF,RM2(u) are the following:p = 0.953 foru = 0.10 mm, p=4.795·10-5 foru = 0.1503 mm,p = 0.5378 foru = 0.200 mm,p = 0.153 foru = 0.250 mm. Only whenu = 0.1503 mm, Figure 12(b), the hypothesis of the normality of r.v. XF,RM2 must be rejected at the significance level α = 0.05; for the rest cases the hypotheses about normality can be accepted; see Figures 12(a), 12(c), and 12(d).Figure 12
Histograms of the tensile forceFRM2 and kernel density estimations of 72 realizations of the force random variable XF,RM2 when u∈{0.1,01503,020,0.250} mm.
(a) (b) (c) (d)The analysis of the normality of the realizations of the r.v.XF,RM2(u), where u∈{0.01+0.0025k:k∈{0,…,36}}={0.01,…,0.1} mm, showedpvalue of the Shapiro-Wilk test,p < α = 0.05, whenu ≤ 0.05 mm. Therefore, it is a reason to claim that r.v. XF,RM2(u) is not a Gaussian whenu ≤ 0.05 mm. Whenu > 0.05 mm, then in majority of cases, at 15 points of the remaining 20,p ≥ 0.05. Therefore, whenu > 0.05 mm, there is a reason to treat r.v. XF,RM2(u) as a Gaussian r.v. Also, on the basis of the performed analysis it is reasonable to claim that the stochastic process {XF,RM2(u):u∈[0,0.05]mm} is not Gaussian, while the stochastic process {XF,RM2(u):u∈(0.05,0.35]mm} may be treated as a Gaussian. Combining these conclusions and Figure 6(b) we can also state that the stochastic process {XF,RM2(u):u∈[0,0.05]mm} is not a Gaussian within the interval in which the new concrete cracks form most intensive. Whenu > 0.05 mm, then the intensity of the formation of the new concrete cracks is less and the stochastic process {XF,RM2(u):u∈(0.05,0.35]mm} may be treated as a Gaussian. The decreasing intensity of the formation of the new concrete cracks whenu > 0.05 mm can be seen very clearly from realization 2 (dotted line) of r.v. XN,br,RM2 shown on Figure 6(b). It should be noted that the considered realizations of r.vs. XF,RM2 and XN,br,RM2 are obtained by assuming that the critical concrete tensile strain r.vs. Xε,crit,i, j∈{1,2,3} are distributed according to truncated normal variable when the coefficient of the variationCvar(Xε,crit,i) = 0.1776. For the other coefficients of variation, the intervals in which the stochastic processes can be treated as a non-Gaussian or Gaussian could be different. However, the performed analysis allows us to claim that r.v. XF,RM2 is not Gaussian within the intervals in which the intensity of crack formation is very big. And on the contrary, r.v. XF,RM2 becomes Gaussian with decreasing the intensity of the formation of the new concrete cracks. It should be noticed that in the common reinforced structures, there are at least three characteristic intervals of the displacements of the cracking of the reinforcement of ties; see Model Code 2010 [41]: the first uncracked interval starts at the zero and ends when the first crack appears. Then the next, second cracking interval starts in which the new cracks appear very intensively. Usually, the force-displacement curve of this interval is uniform like yield plateau of mild steel. The intensity of the formation of the cracks in this displacement interval is of the most intensity. The next, third interval starts when the intensity of the formation of the new cracks becomes lesser. In this interval, the force-displacement curve increases with increasing the displacement of the reinforcement. The obtained force-displacement relationships of the 2D model also have similar intervals, and the relationships are not Gaussian at the big crack formation intensity. Therefore, the conclusion that the force-displacement relationship is not Gaussian at least within the second interval of the cracking can be extended in general also for the common reinforced concrete ties. However, this question demands more deep investigations.The correlation matrices of the correlation coefficientsρi,j = ρ(XF,RM2(ui), XF,RM2(uj)) of the tensile force r.v. XF,RM2 are shown in Figure 13. Figure 13(a) shows ρi,j when ui,uj∈{0.01+0.0025k:k∈0,…,36} and Figure 14(b) when ui,uj∈{0.1+0.01k:k∈{0,…,20}}. ρi,j attains the relatively low values at smallu∈[0.01,0.1] mm in comparison to the ρi,j at the bigger displacementsu∈[0.1, 0.3] mm. For example, ρ(XF,RM2(0.06), XF,RM2(0.1)) = 0.209 while ρ(XF,RM2(0.26), XF,RM2(0.3)) = 0.662. The coefficient of correlation ρ(XF,RM2(ui), XF,RM2(uj)) increases with the decreased distance abs(ui - uj).Figure 13
The correlation matrices of r.v.XF,RM2: when u∈{0.01+0.0025k:k∈{0,…,36}} (a) and u∈{0.1+0.01k:k∈{0,…,20}} (b).
(a) (b)Figure 14
Dependencies of the axial forceF (a) and number of broken springs Nbr (b) on the displacementu.
(a) (b)The weak correlationsρ(XF,RM2(ui), XF,RM2(uj)) shows weak linear dependence between tensile forces FRM2(ui) and XF,RM2(uj). This linear dependence decreases with increasing the distance abs(ui - uj). In the majority of cases, ρi,j > 0. It shows that the increase tensile force FRM2(ui) leads to the increased FRM2(uj), where ui < uj. However, this tendency decreases with the increased distance abs(ui - uj). Despite the fact that the correlation coefficients ρ(XF,RM2(ui), XF,RM2(uj)) decrease with the increased distance abs(ui - uj), we cannot claim that r.v. XF,RM2(uj) becomes less dependent on r.v. XF,RM2(ui) with increasing abs(ui - uj), where uj > ui since Finf (u) ≤F(u) ≤ Fsup (u) and Finf(u) as well as Fsup(u) increase with increasing the displacementu. It should be noticed that the negative coefficients of correlation ρ(XF,RM2(uξ), XF,RM2(uη)) < 0 exist at the displacements (uξ,uη) at which the increasing tensile forceF(uξ) leads to the decreasedF(uη). This decrease is caused by the formation of the new cracks in concrete whenu∈ [uξ,uη].The dependencies of the tensile forcesF∈FRM1,FRD1,FRM2,FDM2,FRM3,FDM3 and the number of the broken springs Nbr∈Nbr,RM1,Nbr,RD1, Nbr,RM2, Nbr,DM2, Nbr,RM3, Nbr,DM3 on the displacement u are depicted in Figure 14. In this figure FRMi and Nbr,RMi, i∈{1,2,3}, denote single realizations of r.vs. XF,RMi and XN,br,RMi. As we can see, Nbr,RMi(u)≥Nbr,DMi(u), where i∈{1,2,3}, u∈[0,0.2] mm, in major part of the interval [0,0.2] mm. However there is no unique tendency on FDMi(u)≥FRMi(u). It is evident that FRM1 ≥ FRM2 ≥ FRM3, FDM1 ≥ FDM2 ≥ FDM3, Nbr,RM1 ≤ Nbr,RM2 ≤ Nbr,RM3, and Nbr,DM1 ≤ Nbr,DM2 ≤ Nbr,DM3 since με,crit,1 > με,crit,2>με,crit,3.Figure14 shows that whenu > 0.162 mm, then Nbr,DM2 < Nbr,RM1 even though με,crit,2>με,crit,1. From Figure 14 we can also see that with decreasing critical tensile strain εcrit the tensile forceF approaches to the Finf. The intense cracks formation stage starts at the bigger displacements of the reinforcement for the specimens with bigger critical concrete tensile strains εcrit, and these displacements are bigger for the deterministic specimens than for random.


## 5.2. Cracking of the Tie

The evaluation of the cracks pattern of the specimens DM2 and DM3, when the critical concrete tensile strainεcrit,2 = 1.298·10-3 and εcrit,3 = 0.322·10-3, at the displacementsu∈(0.04, 0.072, 0.104, 4.0) mm, is shown in Figures 15 and 16. As we can see from these pictures, whenu = 0.04 mm, the debonding crack (the crack between the steel reinforcement and the concrete) is longer for the model DM3 than for DM2. For this case, the debonding crack reaches the second ridge from the left for the model DM2 and the first ridge from the left for the model DM3. Also, we can see that the cracks are higher and more inclined for the model DM2 than for DM3.Figure 15
Cracks patterns of the specimen DM2, whenεcrit,2 = 1.298·10-3, at the displacementsu∈ (0.04, 0.072, 0.104, 4.0) mm.
(a) (b) (c) (d)Figure 16
Cracks patterns of the specimen DM3, whenεcrit,3 = 0.322·10-3, at the displacementsu∈ (0.04, 0.072, 0.104, 4.0) mm.
(a) (b) (c) (d)Cracking patterns of the deterministic, DM1, DM2 and DM3, and random, RM1, RM2 and RM3, models at the displacementu = 0.2 mm are shown in Figure 17.Figure 17
Cracking patterns of the specimens DMi and RMi,i∈ {1,2,3} atu = 0.2 mm: (a), (c), and (e) denote the deterministic models RD1, RD2, and RD3, respectively; (b), (d), and (f) denote the random models RM1, RM2, and RM3, respectively.
(a) (b) (c) (d) (e) (f)From Figure17 we can notice that the spatial variability of the cracks of all single realizations of the random specimens RM1, RM2, and RM3 is more significant in comparison to the deterministic specimens DM1, DM2, and DM3. Also, we can see that the spatial variability of the cracks is more significant for a realization of the specimen RM3 than for RM2 and RM1, Figures 17(f), 17(d), and 17(b) respectively. Therefore, the decreasing critical tensile strains also increase the cracks variability of the concrete. It should be noted that in principle the location and inclination of the modelled concrete cracks, Figures 15, 16, and 17, near the ribs of the reinforcement, are very similar to the cracks obtained experimentally; see Goto and Otsuka [42]. That is the cracks are inclined in the side of the imposed reinforcement force. Also, the angle of the inclination and location of the modelled cracks near the ribs were compared with the numerical results obtained by modelling the identical yet homogeneous and unraced tie using the finite element method. The comparison showed that the location of the cracks near the ribs obtained by the lattice method corresponds to the location of the maximum principal tensile stresses in the ribs region obtained by the finite element method. Moreover, the inclination of the principal planes of the maximum principal tensile stresses is similar to the inclination of the modelled cracks.


## 6. Conclusions

The force-displacement behaviour of the reinforced tie with spatially variable concrete tensile strains was modelled and investigated as a stochastic process using the regular triangular 2D lattice. The analysis showed that the 2D lattice approach, in general, can be applied to model the force-displacement behaviour of the reinforced concrete ties. The stochastic modelling shows that the distribution of the tensile force with respect to the displacement of the reinforcement is not a Gaussian random variable at the small displacements as well as within the intervals of the displacement within which the formation of the new concrete cracks is very intensive,u ≤ 0.05 mm for the investigated specimen with the significance level α = 0.05. Therefore, the tensile force dependent on the displacement, in general, cannot be treated as a Gaussian random process. In the next stage, when the intensity of the cracking of the concrete is much less,u > 0.05 mm for the investigated specimen, the tensile force can be treated as a Gaussian random process when the significance level α = 0.05. The range, standard deviation and the coefficient of the variation of the tensile force increase at the beginning of the deformation of the reinforcement. However, at the particular value of the displacement, the range, standard deviation and the coefficient of variation of the force decrease with the increased displacement. In general, the autocorrelation function decreases with increased difference of the displacements at which the correlations are measured. The decrease of the autocorrelation is sharper at the small displacements of the reinforcement and slower at the bigger displacements, after the stage of the intense formation of the cracks. Decreasing critical tensile strains increases the cracks variability of the tie concrete.



---

*Source: 6340656-2019-04-10.xml*