Simulation of Crack Propagation in Reinforced Concrete Elements

: The crack widths in reinforced concrete structures are affected by various inﬂuencing parameters, such as the tensile strength of the concrete, which is a highly variable parameter. While safe predictions of crack widths are possible with the existing models supplied in the code provisions accompanied by suitable safety factors, there are large discrepancies between the experimentally measured crack widths and the calculated values. Even within a uniaxial tensile test on reinforced concrete components, several crack parameters (spacings and widths) are usually measured. To further investigate the effects of scattering input parameters on the distribution of crack parameters, a crack propagation model (CPM) is developed, which is a powerful tool that can be used to observe sequential cracking of uniaxially loaded reinforced concrete components. In the present paper, the background of the scattering of the concrete tensile strength will be explained and the procedure for the simulation in the CPM, including the investigation of different bond–slip relationships, will be introduced. Finally, the CPM will be used to calculate experimental crack parameters and the distribution of the calculated crack parameters will be discussed.


Introduction
In reinforced concrete structures, cracks are unproblematic, as long as the crack widths are limited to a specific level depending on the environmental conditions. For this purpose, crack width verification is carried out in the design of a structural component, which ensures the requirements are met in terms of the serviceability and durability of the structural component by means of a normatively defined calculated value of the crack width. Different calculation models exist in the respective standards that can be used to verify the calculated crack width, which are nowadays usually mechanically based. Essential differences in the models result from the reference values used for the crack width. In some models, the mean value of the crack widths is used as a reference value and the characteristic crack width is determined using a corresponding increase factor. A further group of models directly calculate the characteristic crack widths. For the comparison of such characteristic crack widths with the measured crack widths in the tests, statistical considerations are required. The 95th percentile of the distribution is associated with the characteristic crack width value. Therefore, for both procedures, sufficient knowledge of the stress-dependent crack formation process is required, as well as knowledge of the resulting scattering and distribution density of the occurring crack widths and crack spacings.
The experimental determination of the stress-dependent crack width distribution is associated with a very high experimental effort, because uniaxial tension tests on long reinforced concrete specimens are required, in which a statistically relevant number of cracks emerge. Such tests were carried out in the 1960s (e.g., [1,2]) and provided a good data basis for structural elements with steel reinforcement and the common concrete properties at that time. For new concrete mixtures, high reinforcement ratios, large-diameter steel bars, or modern non-metallic reinforcement, the amount of available data on existing tests for the statistical distribution of crack parameters (i.e., crack spacing and crack widths) is very limited or does not even exist.
An alternative approach involves the use of numerical models of the crack formation process involving a consistent mechanical background. The performance and prediction accuracy of such models must be proven via validation with experimental tests. In [3,4], efficient models of the simulation of random crack formation have already been developed, which were each developed and used for specific tasks. For more comprehensive considerations of crack formation in reinforced tensile components or flexural tension zones [5], a crack propagation model (CPM) was developed, with the ability to determine the distribution functions of crack widths and crack spacings in a general form.
In this paper, the model setup and features of the CPM are explained and a comparison with experimental tests is carried out for deformed bars. Furthermore, possible extensions of the CPM for investigation of the crack parameters of recently used concrete mixtures (e.g., steel-fiber-reinforced concrete, high-strength concrete, carbon concrete) are presented.

General Crack Propagation
If an uncracked reinforced concrete specimen (uncracked stage) is subjected to a steadily increasing centric tensile load, cracks will appear at the point of the lowest concrete tensile strength (initial cracking stage). In the crack, the concrete strain value is zero and the steel strain value is maximal. Due to the bond stresses activated on both sides of the crack, tensile stresses are introduced into the concrete and the steel tensile stresses are reduced until the steel and concrete strains are identical again at the end of the transmission length l t (Figure 1). With a further increase of the load, the next crack occurs in the remaining uncracked component zone at the point with the lowest concrete tensile strength (crack formation stage; Figure 1, stage 2). This process is repeated until at no point are the steel and concrete strains identical. At this point the stabilized cracking stage is reached ( Figure 1, stage 3), however it must be noted that with further load increase beyond this state, further cracks can form. diameter steel bars, or modern non-metallic reinforcement, the amount of available data on existing tests for the statistical distribution of crack parameters (i.e., crack spacing and crack widths) is very limited or does not even exist.
An alternative approach involves the use of numerical models of the crack formation process involving a consistent mechanical background. The performance and prediction accuracy of such models must be proven via validation with experimental tests. In [3,4], efficient models of the simulation of random crack formation have already been developed, which were each developed and used for specific tasks. For more comprehensive considerations of crack formation in reinforced tensile components or flexural tension zones [5], a crack propagation model (CPM) was developed, with the ability to determine the distribution functions of crack widths and crack spacings in a general form.
In this paper, the model setup and features of the CPM are explained and a comparison with experimental tests is carried out for deformed bars. Furthermore, possible extensions of the CPM for investigation of the crack parameters of recently used concrete mixtures (e.g., steel-fiber-reinforced concrete, high-strength concrete, carbon concrete) are presented.

General Crack Propagation
If an uncracked reinforced concrete specimen (uncracked stage) is subjected to a steadily increasing centric tensile load, cracks will appear at the point of the lowest concrete tensile strength (initial cracking stage). In the crack, the concrete strain value is zero and the steel strain value is maximal. Due to the bond stresses activated on both sides of the crack, tensile stresses are introduced into the concrete and the steel tensile stresses are reduced until the steel and concrete strains are identical again at the end of the transmission length lt (Figure 1). With a further increase of the load, the next crack occurs in the remaining uncracked component zone at the point with the lowest concrete tensile strength (crack formation stage; Figure 1, stage 2). This process is repeated until at no point are the steel and concrete strains identical. At this point the stabilized cracking stage is reached (Figure 1, stage 3), however it must be noted that with further load increase beyond this state, further cracks can form.
The crack width in the crack formation stage results from the integration of the strain differences between the reinforcement and surrounding concrete over the transmission length on both sides of the observed crack. In the stabilized cracking stage, the crack width is determined by integrating the strain difference over half of the crack spacing on both sides of the considered crack.  The crack width in the crack formation stage results from the integration of the strain differences between the reinforcement and surrounding concrete over the transmission length on both sides of the observed crack. In the stabilized cracking stage, the crack width is determined by integrating the strain difference over half of the crack spacing on both sides of the considered crack.
Hence, the concrete tensile strength values (specifically the lower concrete tensile strength values) over the length are significant in terms of initiation of cracks, while the bond-slip behaviour is the major parameter affecting the transmission of stresses between the steel bar and the concrete. Therefore, to simulate the whole crack propagation process with increasing load level (i.e., increasing steel stress σ s ), the crack propagation model (CPM) is configured with the following features:

•
For random cracking, concrete tensile strength scattering along the longitudinal axis of the specimen is required. Therefore, the CPM can map different statistical distributions; • It must be possible to calculate steel and concrete strains separately over the entire length of the specimen. Thus, both materials shall be represented separately; • A bond stress-slip relationship is required for the interaction of both materials. This should be definable in a general mathematical form to be able to take into account different loading levels and material configurations.
Using these features, the crack formation can be considered in a very variable and specific way (e.g., through adaptation of input parameters and calculation steps and isolated consideration of influence parameters with almost any number of calculation runs). Using a random definition for the concrete tensile strength over the specimen length, the random locations of the crack initiation and random crack spacings can be determined. Furthermore, the statistical distribution of the output crack widths and crack spacings can be investigated with several runs and numerous random parameters within a specific range.
However, to obtain a mechanically reasoned model with very good prediction accuracy, specific considerations concerning the tensile strength distribution are needed to reach element-size-independent and reasonable results. Such aspects of the CPM model will be explained in the following.

Model Setup
To implement the features listed in the previous section, the rheological model described in [6,7] was extended. The basic model consists of two uniaxial bars (concrete and reinforcing steel), which are connected with springs at a defined distance ( Figure 2). Here, the springs represent the bond between the steel bar and the concrete bar. The extended crack propagation model (CPM) allows "free crack formation" of the concrete tensile strength scattering along the bar or component axis (cf. models in [3,4]). The scattering of the concrete tensile strength is randomly generated according to a given distribution (e.g., normal distribution, logarithmic normal distribution).
Hence, the concrete tensile strength values (specifically the lower concrete tensile strength values) over the length are significant in terms of initiation of cracks, while the bond-slip behaviour is the major parameter affecting the transmission of stresses between the steel bar and the concrete. Therefore, to simulate the whole crack propagation process with increasing load level (i.e., increasing steel stress s ), the crack propagation model (CPM) is configured with the following features:


For random cracking, concrete tensile strength scattering along the longitudinal axis of the specimen is required. Therefore, the CPM can map different statistical distributions;  It must be possible to calculate steel and concrete strains separately over the entire length of the specimen. Thus, both materials shall be represented separately;  A bond stress-slip relationship is required for the interaction of both materials. This should be definable in a general mathematical form to be able to take into account different loading levels and material configurations.
Using these features, the crack formation can be considered in a very variable and specific way (e.g., through adaptation of input parameters and calculation steps and isolated consideration of influence parameters with almost any number of calculation runs). Using a random definition for the concrete tensile strength over the specimen length, the random locations of the crack initiation and random crack spacings can be determined. Furthermore, the statistical distribution of the output crack widths and crack spacings can be investigated with several runs and numerous random parameters within a specific range.
However, to obtain a mechanically reasoned model with very good prediction accuracy, specific considerations concerning the tensile strength distribution are needed to reach element-size-independent and reasonable results. Such aspects of the CPM model will be explained in the following.

Model Setup
To implement the features listed in the previous section, the rheological model described in [6,7] was extended. The basic model consists of two uniaxial bars (concrete and reinforcing steel), which are connected with springs at a defined distance ( Figure 2). Here, the springs represent the bond between the steel bar and the concrete bar. The extended crack propagation model (CPM) allows "free crack formation" of the concrete tensile strength scattering along the bar or component axis (cf. models in [3,4]). The scattering of the concrete tensile strength is randomly generated according to a given distribution (e.g., normal distribution, logarithmic normal distribution).  In contrast to the rheological model in [6,7], the specified model length l of the CPM corresponds to the entire length of the simulated component. The crack formation is recorded via stiffness reduction of the corresponding concrete element (EA) c,I = 0 ( Figure 3). In contrast to the rheological model in [6,7], the specified model length l of the CPM corresponds to the entire length of the simulated component. The crack formation is recorded via stiffness reduction of the corresponding concrete element (EA)c,I = 0 ( Figure 3). Since a crack can only emerge between two nodes, i.e., in the middle of a concrete bar element (Figure 3), the crack width is calculated by using the slip difference of the left and right nodes, taking into account the steel bar strains. Considering that the strain values in the remaining "concrete stumps" (between concrete nodes KCi and KCi+1) equal zero and the steel strain values in the steel bar element i (between steel nodes KSi and KSi+1) correspond to εs,i = (uKS,i+1 − uKS,i) ⁄ li, the calculated crack width can be simplified as the difference between the translation u of the concrete nodes KCi+1 and KCi: where si is the slip at node i; uKS,i is the translation of the steel bar at node i; uKC,i is the translation of the concrete bar at node i; εs,i is the strain of the steel bar i; and li is the length of element i.
Using such simplified linear displacement methods, the element size should be selected carefully to obtain reasonable results. In the CPM, the error of the determined crack widths is usually less than 3%, with 6 to 8 elements between two cracks. Therefore, the element length le is generally selected with regard to the expected crack spacing sr, which is assumed over a reasonable range for the investigated specimen (e.g., considering the reinforcement ratio, rebar diameter, and bond-slip relationship).

Distribution of Tensile Strength
The experimental determination of tensile strength has been the subject of recent research (e.g., [8,9]), and several indirect and direct testing methods exist. Uniaxial concrete tensile strength is determined most effectively using a direct (uniaxial) tensile test ( Figure 4). To determine the scattering of the concrete tensile strength of a single batch, the statistical distribution of several direct tensile tests is required. This scattering is caused by a wide variety of reasons, such as flaws due to insufficient compaction, the presence of inhomogeneity and anisotropy caused by aggregates, as well as internal stress states caused by shrinkage and hydration of cement.
To describe the distribution of tensile strength, a probability function is required. In [10], a logarithmic normal distribution is assumed. In contrast, [11] and some codes in [12,13] assume a normal distribution, while a summary is provided in [14]. Both distributions are shown in Figure 5   Since a crack can only emerge between two nodes, i.e., in the middle of a concrete bar element (Figure 3), the crack width is calculated by using the slip difference of the left and right nodes, taking into account the steel bar strains. Considering that the strain values in the remaining "concrete stumps" (between concrete nodes KC i and KC i+1 ) equal zero and the steel strain values in the steel bar element i (between steel nodes KS i and KS i+1 ) correspond to ε s,i = (u KS,i+1 − u KS,i ) ⁄ l i , the calculated crack width can be simplified as the difference between the translation u of the concrete nodes KC i+1 and KC i : where s i is the slip at node i; u KS,i is the translation of the steel bar at node i; u KC,i is the translation of the concrete bar at node i; ε s,i is the strain of the steel bar i; and l i is the length of element i. Using such simplified linear displacement methods, the element size should be selected carefully to obtain reasonable results. In the CPM, the error of the determined crack widths is usually less than 3%, with 6 to 8 elements between two cracks. Therefore, the element length l e is generally selected with regard to the expected crack spacing s r , which is assumed over a reasonable range for the investigated specimen (e.g., considering the reinforcement ratio, rebar diameter, and bond-slip relationship).

Distribution of Tensile Strength
The experimental determination of tensile strength has been the subject of recent research (e.g., [8,9]), and several indirect and direct testing methods exist. Uniaxial concrete tensile strength is determined most effectively using a direct (uniaxial) tensile test ( Figure 4). To determine the scattering of the concrete tensile strength of a single batch, the statistical distribution of several direct tensile tests is required. This scattering is caused by a wide variety of reasons, such as flaws due to insufficient compaction, the presence of inhomogeneity and anisotropy caused by aggregates, as well as internal stress states caused by shrinkage and hydration of cement.
To describe the distribution of tensile strength, a probability function is required. In [10], a logarithmic normal distribution is assumed. In contrast, [11] and some codes in [12,13] assume a normal distribution, while a summary is provided in [14]. Both distributions are shown in Figure 5    It can be seen that both distributions are very similar for a coefficient of variation vk,1 = 10%, with an overlap of 96.2%. For vk,2 = 20%, the overlap is 92.4%, but the visible deviation in Figure 5 is significantly larger than for vk,1 = 10%. Consequently, the logarithmic normal distribution for vk < 20% can be approximately represented by a normal distribution.
For ranges of variation for different expected values, a distinction must be made between a constant standard deviation and a constant coefficient of variation. For the concrete tensile strength-in contrast to the concrete compressive strength (standard deviation constant [15,16])-the coefficient of variation is considered a constant quantity. With regard to the scattering of the tensile strength at failure, evaluations in [11,17,18] show a coefficient of variation vk < 10%. In Model Code 2010 [12] and Eurocode 2 [13], the coefficient of variation of the tensile strength of the component is vk = 18% (see [4]). Similar results are demonstrated in [19], by considering the influences of different concrete mixtures [17]. Thus, for one concrete batch, a coefficient of variation vk < 10% seems to be realistic. For several concrete batches in a component or in the absence of tensile tests, vk = 18% may be reasonable. For the probability function of the tensile strength, however, a normal distribution can be assumed approximately in both cases.
The aforementioned distribution function applies to the tensile strength values obtained from several uniaxial tensile tests and cannot be directly implemented to define a tensile strength distribution in the CPM to simulate a tensile bar. This is due to the fact  It can be seen that both distributions are very similar for a coefficient of variation vk,1 = 10%, with an overlap of 96.2%. For vk,2 = 20%, the overlap is 92.4%, but the visible deviation in Figure 5 is significantly larger than for vk,1 = 10%. Consequently, the logarithmic normal distribution for vk < 20% can be approximately represented by a normal distribution.
For ranges of variation for different expected values, a distinction must be made between a constant standard deviation and a constant coefficient of variation. For the concrete tensile strength-in contrast to the concrete compressive strength (standard deviation constant [15,16])-the coefficient of variation is considered a constant quantity. With regard to the scattering of the tensile strength at failure, evaluations in [11,17,18] show a coefficient of variation vk < 10%. In Model Code 2010 [12] and Eurocode 2 [13], the coefficient of variation of the tensile strength of the component is vk = 18% (see [4]). Similar results are demonstrated in [19], by considering the influences of different concrete mixtures [17]. Thus, for one concrete batch, a coefficient of variation vk < 10% seems to be realistic. For several concrete batches in a component or in the absence of tensile tests, vk = 18% may be reasonable. For the probability function of the tensile strength, however, a normal distribution can be assumed approximately in both cases.
The aforementioned distribution function applies to the tensile strength values obtained from several uniaxial tensile tests and cannot be directly implemented to define a tensile strength distribution in the CPM to simulate a tensile bar. This is due to the fact It can be seen that both distributions are very similar for a coefficient of variation v k,1 = 10%, with an overlap of 96.2%. For v k,2 = 20%, the overlap is 92.4%, but the visible deviation in Figure 5 is significantly larger than for v k,1 = 10%. Consequently, the logarithmic normal distribution for v k < 20% can be approximately represented by a normal distribution.
For ranges of variation for different expected values, a distinction must be made between a constant standard deviation and a constant coefficient of variation. For the concrete tensile strength-in contrast to the concrete compressive strength (standard deviation constant [15,16])-the coefficient of variation is considered a constant quantity. With regard to the scattering of the tensile strength at failure, evaluations in [11,17,18] show a coefficient of variation v k < 10%. In Model Code 2010 [12] and Eurocode 2 [13], the coefficient of variation of the tensile strength of the component is v k = 18% (see [4]). Similar results are demonstrated in [19], by considering the influences of different concrete mixtures [17]. Thus, for one concrete batch, a coefficient of variation v k < 10% seems to be realistic. For several concrete batches in a component or in the absence of tensile tests, v k = 18% may be reasonable. For the probability function of the tensile strength, however, a normal distribution can be assumed approximately in both cases.
The aforementioned distribution function applies to the tensile strength values obtained from several uniaxial tensile tests and cannot be directly implemented to define a tensile strength distribution in the CPM to simulate a tensile bar. This is due to the fact that in uniaxial tensile tests (of unnotched specimens), the tensile strength applies to the weakest point of the specimen. In the following, this value will be designated as the fracture tensile strength f ct,f . Previous experimental investigations (e.g., [17,20]) show a difference between the distribution of the tensile strength within a specimen and the distribution of the fracture tensile strength f ct,f . Regarding this point, Empelmann [3] defined a term for the specimen tensile strength f ct , which is explained below. Figure 6a shows the fracture process in a uniaxial tensile specimen. The existing micro-cracks accumulate at the weakest point of the specimen shortly before the maximum load is reached [21]. The minimum tensile strength shown in Figure 6b corresponds to the fracture tensile strength f ct,min = f ct,f . Figure 6b schematically shows the tensile strength values in representative elements (which will be discussed later). Based on this point, it is logical that the mean value of all f ct,i values is higher than the mean value of several f ct,f values obtained using several uniaxial tensile tests. that in uniaxial tensile tests (of unnotched specimens), the tensile strength applies to the weakest point of the specimen. In the following, this value will be designated as the fracture tensile strength fct,f. Previous experimental investigations (e.g., [17,20]) show a difference between the distribution of the tensile strength within a specimen and the distribution of the fracture tensile strength fct,f. Regarding this point, Empelmann [3] defined a term for the specimen tensile strength fct, which is explained below. Figure 6a shows the fracture process in a uniaxial tensile specimen. The existing micro-cracks accumulate at the weakest point of the specimen shortly before the maximum load is reached [21]. The minimum tensile strength shown in Figure 6b corresponds to the fracture tensile strength fct,min = fct,f. Figure 6b schematically shows the tensile strength values in representative elements (which will be discussed later). Based on this point, it is logical that the mean value of all fct,i values is higher than the mean value of several fct,f values obtained using several uniaxial tensile tests. Consequently, an adaption is required to derive the distribution of the specimen tensile strength fct (all values in a component) based on the given distribution of the fracture tensile strength values fct,f (weakest values in a component in uniaxial tensile specimens without a notch).
According to [3], such adaption can be built by assuming a normal distribution for fct,f and fct and setting the mean value of the specimen tensile strength fct,m equal to the 95th percentile value of the fracture tensile strength fct,f,95%.
For the mean value of the fracture concrete tensile strength fct,f,m = 2.0 N/mm 2 , the following percentile values are obtained by setting the assumed coefficient of variation (of the fracture tensile strength) as equal to vk,f = 10%: The mean value of the specimen tensile strength corresponds to: ct,m = ct,f,95% = 2.33 N mm 2 ⁄ Assuming equal coefficient variations vk = vk,f = 10%: Consequently, an adaption is required to derive the distribution of the specimen tensile strength f ct (all values in a component) based on the given distribution of the fracture tensile strength values f ct,f (weakest values in a component in uniaxial tensile specimens without a notch).
According to [3], such adaption can be built by assuming a normal distribution for f ct,f and f ct and setting the mean value of the specimen tensile strength f ct,m equal to the 95th percentile value of the fracture tensile strength f ct,f,95% .
For the mean value of the fracture concrete tensile strength f ct,f,m = 2.0 N/mm 2 , the following percentile values are obtained by setting the assumed coefficient of variation (of the fracture tensile strength) as equal to v k,f = 10%: The mean value of the specimen tensile strength corresponds to: Assuming equal coefficient variations v k = v k,f = 10%: where f ct,f,m is the mean value of the fracture tensile strength; v k,f is the coefficient of variation of the fracture tensile strength; f ct,m is the mean value of the specimen tensile strength; v k is the coefficient of variation of the specimen tensile strength; 5% is the 5th percentile value; and 95% is the 95th percentile value. This calculation method is depicted in Figure 7.
where fct,f,m is the mean value of the fracture tensile strength; vk,f is the coefficient of variation of the fracture tensile strength; fct,m is the mean value of the specimen tensile strength; vk is the coefficient of variation of the specimen tensile strength; 5% is the 5th percentile value; and 95% is the 95th percentile value. This calculation method is depicted in Figure 7.

Element Size and Autocorrelation
The concrete bar of the CPM consists of a finite elements series, with each element being assigned a random tensile strength fct,i according to the applied distribution. As shown in Section 3.2 ( Figure 6), each tensile strength value is assigned according to a mechanically reasoned representative element size ∆x, which is often defined based on the maximum aggregate size dg. Here, ∆x is assumed to be in the range of 2•dg to 3•dg [3,22,23]. A maximum aggregate size of 16 mm results in a representative element size ∆x of about 32-48 mm. The element size affects directly the prediction accuracy of the model and the accuracy of the crack spacing values, as a crack can only emerge in the mid-length area of each element. Therefore, to enable the simulation of fine cracking patterns, the element length le should also be very small. To consider these aspects and simultaneously achieve a nearly mesh-independent model, a smaller element length le is selected and the concrete tensile strengths of the adjacent elements within a correlation length Lx are correlated to avoid mechanically inconsistent oscillation of the tensile strengths in the smaller elements. This approach is known as spatial autocorrelation [24].
The expression of the linear dependence or correlation between two points x1 and x2 is described by the correlation function ρ(x1,x2). A common form is the quadratic exponential (Gaussian) correlation function [4,25,26]. Figure 8 illustrates the correlation function according to Equation (7). For ρ(x1,x2) = 0, there is no linear correlation between x1 and x2, whereas ρ(x1,x2) = 1 shows a full dependency of both points.

Element Size and Autocorrelation
The concrete bar of the CPM consists of a finite elements series, with each element being assigned a random tensile strength f ct,i according to the applied distribution. As shown in Section 3.2 ( Figure 6), each tensile strength value is assigned according to a mechanically reasoned representative element size ∆x, which is often defined based on the maximum aggregate size d g . Here, ∆x is assumed to be in the range of 2·d g to 3·d g [3,22,23]. A maximum aggregate size of 16 mm results in a representative element size ∆x of about 32-48 mm. The element size affects directly the prediction accuracy of the model and the accuracy of the crack spacing values, as a crack can only emerge in the mid-length area of each element. Therefore, to enable the simulation of fine cracking patterns, the element length l e should also be very small. To consider these aspects and simultaneously achieve a nearly mesh-independent model, a smaller element length l e is selected and the concrete tensile strengths of the adjacent elements within a correlation length L x are correlated to avoid mechanically inconsistent oscillation of the tensile strengths in the smaller elements. This approach is known as spatial autocorrelation [24].
The expression of the linear dependence or correlation between two points x 1 and x 2 is described by the correlation function ρ(x 1 ,x 2 ). A common form is the quadratic exponential (Gaussian) correlation function [4,25,26]. Figure 8 illustrates the correlation function according to Equation (7). For ρ(x 1 ,x 2 ) = 0, there is no linear correlation between x 1 and x 2 , whereas ρ(x 1 ,x 2 ) = 1 shows a full dependency of both points.  By multiplying the correlation function with the standard deviations σ(x1) and σ(x2) (here σ(x1) = σ(x2) = σ), the value of covariance is obtained and the positive symmetric covariance matrix C can be composed: The generation of normally distributed correlated random numbers can be performed by means of the Cholesky decomposition. When using the Cholesky decomposition, normally distributed uncorrelated random values (a1, a2, … am) with a mean value of zero and a standard deviation of one must first be generated. These represent the size of the scattering area.
The Cholesky decomposition of the correlation matrix leads to: Subsequently, the normally distributed random values (a1, a2, … am) are multiplied with the triangular matrix L: This results in a normally distributed correlated distribution, which must be added to the mean value of the desired distribution (here this is fct,m) [25,27]: As laid out before, Lx should be defined in order to give a meaningful correlation by considering the representative element size and so as to generate a good dependency with regard to the correlation function at the same time. According to the correlation function in Figure 8, there is almost no correlation in the case of le/Lx = 2. This means that almost uncorrelated random values are generated for Lx = le/2. Based on Section 3.2, this should be valid for a representative element size ∆x of 2•dg to 3•dg. Therefore, the correlation length can subsequently be assumed according to: By multiplying the correlation function with the standard deviations σ(x 1 ) and σ(x 2 ) (here σ(x 1 ) = σ(x 2 ) = σ), the value of covariance is obtained and the positive symmetric covariance matrix C can be composed: The generation of normally distributed correlated random numbers can be performed by means of the Cholesky decomposition. When using the Cholesky decomposition, normally distributed uncorrelated random values (a 1 , a 2 , . . . a m ) with a mean value of zero and a standard deviation of one must first be generated. These represent the size of the scattering area.
The Cholesky decomposition of the correlation matrix leads to: Subsequently, the normally distributed random values (a 1 , a 2 , . . . a m ) are multiplied with the triangular matrix L: This results in a normally distributed correlated distribution, which must be added to the mean value of the desired distribution (here this is f ct,m ) [25,27]: As laid out before, L x should be defined in order to give a meaningful correlation by considering the representative element size and so as to generate a good dependency with regard to the correlation function at the same time. According to the correlation function in Figure 8, there is almost no correlation in the case of l e /L x = 2. This means that almost uncorrelated random values are generated for L x = l e /2. Based on Section 3.2, this should be valid for a representative element size ∆x of 2·d g to 3·d g . Therefore, the correlation length can subsequently be assumed according to: which means a correlation length L x of between 16 and 24 mm for an aggregate size d g = 16 mm. Figure 9 illustrates the influence of the correlation length for uncorrelated and correlated distributions of the tensile strength along the model axis for element sizes of 5 and 20 mm. All distributions are generated randomly according to the values specified in Figure 9. For uncorrelated values, the distribution of l e = 5 mm has too many local minima and maxima compared to l e = 20 mm. The distributions along the axis are not comparable (although the probability density functions are almost the same), leading to a mesh dependency. In the case of correlated values, there is almost no mesh dependency, as the numbers of local minima and maxima are similar for l e = 20 mm and l e = 5 mm.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 19 1.0 ⋅ g ≤ x ≤ 1.5 ⋅ g (12) which means a correlation length Lx of between 16 and 24 mm for an aggregate size dg = 16 mm. Figure 9 illustrates the influence of the correlation length for uncorrelated and correlated distributions of the tensile strength along the model axis for element sizes of 5 and 20 mm. All distributions are generated randomly according to the values specified in Figure 9. For uncorrelated values, the distribution of le = 5 mm has too many local minima and maxima compared to le = 20 mm. The distributions along the axis are not comparable (although the probability density functions are almost the same), leading to a mesh dependency. In the case of correlated values, there is almost no mesh dependency, as the numbers of local minima and maxima are similar for le = 20 mm and le = 5 mm.

Bond-Slip Relationship
The bond-slip behaviour between the concrete bar and reinforcement bar is simulated in the CPM using springs with a stiffness equal to: where τb(s) is the local bond stress; s is the local slip (relative displacement between concrete and rebar); us is the perimeter of the reinforcement bar; and le is the element length. A bond-slip relationship according to Equation (14) is chosen in the first instance, with which a large variety of bond laws can be integrated: This correlation corresponds to the ascending branch of the bond law according to Model Code 2010 [12]. For calculations in the serviceability limit state, using the ascending branch of the bond-slip relationship is sufficient according to [28]. This considers the expected maximum crack width of 0.4 mm, which corresponds to a slip of 0.2 mm at both sides of the cracks, which is clearly within the ascending branch of the bond-slip relationship.
For the general form of the bond-slip relationship according to Equation (14), some reference values based on the technical literature are listed in Table 1 and are illustrated in Figure 10a for concrete strength class C30/37.

Bond-Slip Relationship
The bond-slip behaviour between the concrete bar and reinforcement bar is simulated in the CPM using springs with a stiffness equal to: where τ b (s) is the local bond stress; s is the local slip (relative displacement between concrete and rebar); u s is the perimeter of the reinforcement bar; and l e is the element length. A bond-slip relationship according to Equation (14) is chosen in the first instance, with which a large variety of bond laws can be integrated: This correlation corresponds to the ascending branch of the bond law according to Model Code 2010 [12]. For calculations in the serviceability limit state, using the ascending branch of the bond-slip relationship is sufficient according to [28]. This considers the expected maximum crack width of 0.4 mm, which corresponds to a slip of 0.2 mm at both sides of the cracks, which is clearly within the ascending branch of the bondslip relationship.
For the general form of the bond-slip relationship according to Equation (14), some reference values based on the technical literature are listed in Table 1 and are illustrated in Figure 10a for concrete strength class C30/37.  Even if the real bond behaviour differs significantly from the bond-slip relationship in Equation (14), user-defined bond-slip relationships can be defined in the CPM using slip values and the corresponding bond stress values. Between the points, a linear interpolation is used to determine the missing values. The only prerequisite here is the starting point at the coordinate origin (τb(s = 0) = 0). An exemplarily user-defined bondslip relationship is shown in Figure 10b for a splitting failure.  (14) and Table 1; (b) user-defined.
A reduction of the local bond stress in the vicinity of a separation crack, e.g., due to internal cracking, can be considered by using a linear reduction factor according to [12]: b,red (s, ) = λ( ) ⋅ b ( ) (15) where: λ(x) is the factor for a reduction of the local bond stress at distance x from the separation crack; x is the distance from the crack; and τb(s) is the local bond stress (e.g., according to Equation (14)).
In Model Code 2010 [12], λ(x) is defined by: As an alternative to the linear reduction factor according to Equations (15) and (16), it is also possible to take into account the zone of the disturbed bond using a regionally separated bond stress-slip relationship according to [33][34][35]. It must be noted that the numerical calculation effort and the computation time increase significantly with the used model. Even if the real bond behaviour differs significantly from the bond-slip relationship in Equation (14), user-defined bond-slip relationships can be defined in the CPM using slip values and the corresponding bond stress values. Between the points, a linear interpolation is used to determine the missing values. The only prerequisite here is the starting point at the coordinate origin (τ b (s = 0) = 0). An exemplarily user-defined bond-slip relationship is shown in Figure 10b for a splitting failure.
A reduction of the local bond stress in the vicinity of a separation crack, e.g., due to internal cracking, can be considered by using a linear reduction factor according to [12]: where: λ(x) is the factor for a reduction of the local bond stress at distance x from the separation crack; x is the distance from the crack; and τ b (s) is the local bond stress (e.g., according to Equation (14)).
In Model Code 2010 [12], λ(x) is defined by: As an alternative to the linear reduction factor according to Equations (15) and (16), it is also possible to take into account the zone of the disturbed bond using a regionally separated bond stress-slip relationship according to [33][34][35]. It must be noted that the numerical calculation effort and the computation time increase significantly with the used model.

Model Features and Advantages
The model is characterized by its fast and stable calculation of cracking. In addition, the CPM is easily extendable, e.g., to simulate cracking of fiber-reinforced concrete (FRC) by considering the fiber effects as a crack-width-dependent additional element load. Furthermore, simulation of high-performance concrete (HPC)and ultra-high performance concrete (UHPC) is possible by setting the properties of the concrete bar and adjusting the bond-slip relationship. Additionally, an extension to a three-bar model (see Figure 11) is possible, which enhances the simulation of elements with different reinforcing elements (e.g., elements with prestressed conventional steel rebars). Non-metallic reinforcement (such as carbon concrete, basalt, or glass fiber reinforcements) can be accommodated by adjusting the material parameters of the "reinforcement bar" and the bond-slip relationship. By adjusting the support conditions, constrained stresses can also be investigated. A transfer to components subjected by flexure is possible if the tension and compression zone are treated separately [5].

Model Features and Advantages
The model is characterized by its fast and stable calculation of cracking. In addition, the CPM is easily extendable, e.g., to simulate cracking of fiber-reinforced concrete (FRC) by considering the fiber effects as a crack-width-dependent additional element load. Furthermore, simulation of high-performance concrete (HPC)and ultra-high performance concrete (UHPC) is possible by setting the properties of the concrete bar and adjusting the bond-slip relationship. Additionally, an extension to a three-bar model (see Figure 11) is possible, which enhances the simulation of elements with different reinforcing elements (e.g., elements with prestressed conventional steel rebars). Non-metallic reinforcement (such as carbon concrete, basalt, or glass fiber reinforcements) can be accommodated by adjusting the material parameters of the "reinforcement bar" and the bond-slip relationship. By adjusting the support conditions, constrained stresses can also be investigated. A transfer to components subjected by flexure is possible if the tension and compression zone are treated separately [5]. In addition to the listed expansion possibilities for different materials, the CPM can be expanded with regard to the scattering of the material parameters. In particular, the scattering of the bond stress should be mentioned here. Other scattering parameters (e.g., the modulus of elasticity of the rebar and concrete) are expected to be less significant or can be represented by already existing scattering parameters (e.g., scattering dimensions of the concrete bar by scattering of the tensile strength, scattering diameter of the rebar by scattering of bond stress).

Simulation and Calculation Procedure
With the preceding explanations, the CPM can be used to simulate the cracking of a reinforced concrete tensile specimen. In the following, the programme sequence will be briefly explained. A representation can be taken from Figure 12.
In the first calculation step, the entire model or component is uncracked. A load increase follows until the concrete tensile strength is exceeded at the "weakest" point. Then, the stiffness of the cracked concrete element is reduced to zero ((EA)c,i = 0) (also see Figure 3). The model is recalculated with the adjusted stiffness value (cracked concrete element). The load is then increased further until the concrete tensile strength is exceeded at another point. The process is repeated until the desired load level Ns,end is reached.
The iterative calculation procedure in the CPM will be illustrated below using an example. The calculations are based on the parameters listed in Table 2. In addition to the listed expansion possibilities for different materials, the CPM can be expanded with regard to the scattering of the material parameters. In particular, the scattering of the bond stress should be mentioned here. Other scattering parameters (e.g., the modulus of elasticity of the rebar and concrete) are expected to be less significant or can be represented by already existing scattering parameters (e.g., scattering dimensions of the concrete bar by scattering of the tensile strength, scattering diameter of the rebar by scattering of bond stress).

Simulation and Calculation Procedure
With the preceding explanations, the CPM can be used to simulate the cracking of a reinforced concrete tensile specimen. In the following, the programme sequence will be briefly explained. A representation can be taken from Figure 12.
In the first calculation step, the entire model or component is uncracked. A load increase follows until the concrete tensile strength is exceeded at the "weakest" point. Then, the stiffness of the cracked concrete element is reduced to zero ((EA) c,i = 0) (also see Figure 3). The model is recalculated with the adjusted stiffness value (cracked concrete element). The load is then increased further until the concrete tensile strength is exceeded at another point. The process is repeated until the desired load level N s,end is reached.
The iterative calculation procedure in the CPM will be illustrated below using an example. The calculations are based on the parameters listed in Table 2.     Figure 13 shows the strains determined with the CPM for different load levels (expressed by the maximum steel stress σ s,max ). At the cracks, the concrete stress is zero. The first crack occurs with a steel stress of σ s,max = 100 N/mm 2 at the point of lowest concrete tensile strength, depending on the distribution of the specimen tensile strength (for comparison, without scattering of the concrete tensile strength v k = 0, cracking would not occur until σ s,max = 125 N/mm 2 ). Outside the transmission length, the steel and concrete stresses between the cracks are identical until the completed crack pattern with σ s,max = 114 N/mm 2 is reached. Up to this loading level five cracks are formed, mainly at the low points of the strength distribution. Due to the non-linear bond law, a further increase in the load level leads to an increase in the bond forces, and thus also to a greater introduction of force into the concrete. This leads to further cracking occurring in the stabilized cracking stage (progressive stabilized cracking). These cracks predominantly occur centrally between the existing cracks, so that under certain circumstances, regions with higher strength values are also affected by crack formation, i.e., "late" cracks can emerge under relatively high stresses.

Comparison with Experiments
The performance of the CPM will be demonstrated in the following using a series of tests from [36]. The test series included three identical large-scale tests. The focus of the The individual crack spacing and crack widths for each stress level can then be read out from the CPM calculation results. If several calculation runs are carried out with the newly generated specimen tensile strength distributions along the model axis in each case, extensive distributions of crack spacing and crack widths can be generated as functions of the stress level and a large data basis can be obtained. In contrast to the stages defined in Figure 1, the switch between crack formation stages and the stabilized cracking stage is quite smooth. Hence, stress regions causing fast crack formation can be identified.

Comparison with Experiments
The performance of the CPM will be demonstrated in the following using a series of tests from [36]. The test series included three identical large-scale tests. The focus of the tests was on the generation of many cracks and investigation of the crack spacing distribution. For this purpose, a very high reinforcement ratio of ρ = 8% with a bar diameter of 32 mm was chosen. A tensile force of 900 kN was required to reach a stress level of σ s,max = 280 N/mm 2 . Figure 14 shows the geometrical parameters of the test specimens used in [36].

Comparison with Experiments
The performance of the CPM will be demonstrated in the following using a series of tests from [36]. The test series included three identical large-scale tests. The focus of the tests was on the generation of many cracks and investigation of the crack spacing distribution. For this purpose, a very high reinforcement ratio of ρ = 8% with a bar diameter of 32 mm was chosen. A tensile force of 900 kN was required to reach a stress level of σs,max = 280 N/mm 2 . Figure 14 shows the geometrical parameters of the test specimens used in [36]. The material parameters of the concrete were determined in [36] using compression tests and splitting tensile tests (no data available on scattering). The maximum aggregate size was 16 mm. The parameters used in the CPM are listed in Table 3. A correction of the tensile strength according to Section 3.2 was not performed, as splitting tensile tests were conducted (here f ct,f,m = f ct,m ). The calculations were first carried out for a coefficient of variation of v k = 10% and a correlation length of 20 mm (~1.25 d g , cf. Equation (12)). In a second calculation run, a correlation length of 10 mm was applied for comparison and to consider the small crack spacings. To account for the uncommonly high reinforcement ratio and large reinforcement diameter in combination with the relatively small crack spacing, a coefficient of variation of v k = 15% was also applied. Since in [36] the crack spacing was recorded on 3 sides at the depths of each reinforcement layer (3 specimens × 3 surfaces × 2 bars = 18 sets), 18 generated distributions were recorded for each calculation run. To take the disturbed stress region at loading points into account, both the values presented in [36] and the CPM did not take into account the last crack spacing on both sides of the test specimens.
The distribution of experimental crack spacing and the calculated crack spacing with the CPM is shown in Figure 15a for v k = 10% and in Figure 15b for v k = 15%. The statistical parameters are summarized in Table 4. A calculation with twice as many elements (l e = 5 mm) showed similar results, with the exception of small deviations in the distribution of the crack spacings due to the randomly generated tensile strength distributions.  The comparison of the experimental and calculated crack spacing shows the following results:


In general, the test results can be reproduced very well with the CPM (see Figure 15 and Table 4);  Varying the coefficient of variation and correlation length to determine the distribution of the tensile strength along the model axis has only a limited effect on the results;  The range of very small crack spacings ( r ≤ 100 mm) is underestimated for vk = 10% and vk = 15%, as well as for Lx = 20 mm and Lx = 10 mm (see Figure 15);  The experimental value sr,0.95/sr,m = 1.5 is estimated with adequate accuracy by the CPM (sr,0.95/sr,m = 1.4) (see Table 4);  The comparison of the median and mean values of the experimental results show a positive skew (skewed to the right) (see Table 4) indicating a log-normal distribution of the crack spacing sr.  The comparison of the experimental and calculated crack spacing shows the following results:

•
In general, the test results can be reproduced very well with the CPM (see Figure 15 and Table 4); • Varying the coefficient of variation and correlation length to determine the distribution of the tensile strength along the model axis has only a limited effect on the results; • The range of very small crack spacings (s r ≤ 100 mm) is underestimated for v k = 10% and v k = 15%, as well as for L x = 20 mm and L x = 10 mm (see Figure 15);

•
The experimental value s r,0.95 /s r,m = 1.5 is estimated with adequate accuracy by the CPM (s r,0.95 /s r,m = 1.4) (see Table 4); • The comparison of the median and mean values of the experimental results show a positive skew (skewed to the right) (see Table 4) indicating a log-normal distribution of the crack spacing s r .
The high number of very small crack spacings (s r ≤ 100 mm) can be attributed to the fact that complete separation cracks do not always form in testing and that overlapping cracks can occur. Both phenomena can be observed in [36]. Furthermore, the effects of shrinkage strains have not been taken into account in the calculation with the CPM, as the respective input values were not documented in [36]. However, the model can consider shrinkage by adding shrinkage strains as an additional load component to reduce the effective tensile concrete strength [6,7]. As a result, the crack spacing would be reduced.

Discussion and Conclusions
In the present paper, a simulation model of the random crack propagation in reinforced concrete elements is presented. The model consists of two uniaxial bars (concrete and reinforcing steel) divided into finite elements, which are connected with springs. For the random cracking, random concrete tensile strength values are assigned to the element. Several features and advantages of this model are introduced in the paper, the main objective of which is to discuss the major role of the assumed tensile strength and its distribution, as well as the element size considerations in the CPM.
The most important parameter to simulate random cracking is the uniaxial tensile strength along the model axis. It was pointed out that a normal distributed concrete tensile strength can be used for a common range of coefficients of variation (v k < 20%). This was confirmed by investigations in [4], according to which the assumed distribution type of the tensile strength has only a minor effect on the calculated crack widths. In CPM, different coefficients of variation can be assigned to the normal distribution. The model also distinguishes between the mean value of the fracture tensile strength f ct,f and the mean value of specimen tensile strength f ct,m . If several tests using a test method where the fracture plane develops randomly at the weakest point of the specimen (e.g., direct tensile test without a notch) are available, then the test results represent only the distribution of the weakest tensile strength values f ct,f (here this is the fracture tensile strength). To capture the tensile strength distribution in an entire bar along the longitudinal axis with weak and strong parts (here this is the specimen tensile strength), increasing the mean value of specimen tensile strength f ct,m seems logical. For this purpose, an approach was proposed in Section 3.2, setting f ct,m as equal to the 95th percentile of fracture tensile strength values f ct,f based on [3]. This approach requires further experimental evidence; however, there are no common alternative approaches available in the literature. Since the approach was applied for direct tensile tests (without a notch), an application for splitting tensile tests (as in experiments in [36]) is not recommended without further research.
For numerical models, the results are supposed to be mesh-independent. Since the number of tensile strength values along the model axis in the CPM is dependent on the number of elements (and the respective element lengths), tensile strength values are not supposed to scatter randomly between adjacent elements. For this purpose, a dependency between adjacent elements was defined in the model (autocorrelation). The correlation length-present in the correlation function-specifies the scope of this dependency. The correlation length is a mathematical value, which when combined with mechanical values, enables a mesh-independent and realistic distribution of the tensile strength along the longitudinal axis of a specimen. Since realistic data on the correlation length for tensile strength are scarce [4], for each problem different correlation length values are available in the literature [37]. Therefore, the representative length was introduced in this paper with a mechanically based value considering the aggregate size in order to stay within the macroscale.
The comparison with the experimental test in [36] shows a good reproduction of the results with the CPM, although small crack spacings are slightly underestimated by the CPM. The main reason for this underestimation is probably the influence of shrinkage. The investigations in [38] confirm a high influence of shrinkage strains, especially for high reinforcement ratios. However, no information on shrinkage strain is available in [36]. Nevertheless, the tests in [36] were selected for the comparison, as detailed information on crack spacing distributions is very rare in recent studies (the same applies for crack widths, which were not recorded in [36]).
The CPM also enables investigation of the entire sequential cracking process, which is expected to show a higher dependency of the results on the tensile strength distribution under low load levels. In further investigations, the CPM will be used to generate load deformation curves of reinforced concrete ties (as in [39]). The results can be used to define a damage propagation index based on an increasing number of cracks or decreasing crack spacing (compare Figure 16) (cf. [40]). Based on the index, the load-level-dependent progress of the damage can be derived, including the statistical distribution of the damage parameters.