Assessing the Instantaneous Stiffness of Cracked Reinforced Concrete Beams Based on a Gradual Change in Strain Distributions

Concrete cracking causes a gradual change in strain distributions along the cross section height of reinforced concrete beams, which will finally affect their instantaneous stiffness. A method for assessing the stiffness is proposed based on the gradual change, which is considered through modeling different strain distributions for key sections in cracked regions. Internal force equilibria are adopted to find a solution to top strains and neutral axes in the models, and then the inertias of the key sections are calculated to assess the beam stiffness. (e proposed method has been validated using experimental results obtained from tests on five reinforced concrete beams. (e predicted stiffness and displacements are shown to provide a good agreement with experimental data. (e instantaneous stiffness is proven to greatly depend on the crack number and depth. (is dependence can be exactly reflected by the proposed method through simulating the gradual change in concrete strain distributions.


Introduction
To ensure the serviceability of reinforced concrete beams, stiffness control is an important design objective. However, cracks tend to appear in these beams due to the low tensile strength of concrete and greatly affect the beam stiffness through changing the shape or size of cross sections. So in the last four decades, many researchers focused on evaluating the stiffness of cracked concrete beams and predicting their deflection. e effective inertia of cross section is widely adopted by engineers in design for the deflection control of cracked reinforced concrete beams. Branson [1] developed a model based on a weighted average of two types of inertias, respectively, representing the uncracked and fully cracked inertias of reinforced concrete cross sections. But comparisons with experimental data showed that Branson's model underestimated the instantaneous deflection of lightly reinforced beams and slabs after cracking [2,3]. Bischoff [4] modified Branson's model by introducing empirical factors that effectively change the ratio of the two inertias. Bischoff's model is recommended by fib Model Code 2010 [5], which is now the most comprehensive code on concrete structures including their complete life cycle and also a basis for their future codes. ese models depend on empirical data to determine the weight of the cracked inertia and do not consider actual cracking patterns, so they just can be used at design stages.
To consider the effect of crack location on the beam stiffness, François et al. [6] proposed a macro-finite-element method characterized by the homogenized average inertias of cross sections, which were calculated based on the concept of a transfer length necessary for the transmission of tensile loads from steel to concrete, thanks to a steel-concrete bond. Castel et al. [7] conditioned a cracked span by using cracking moments and maximum moments and developed the macro-finite-element method through assuming the distribution of steel and concrete strains in the cracked span. By using a damage variable D ccc , Castel et al. [8] took into account the effects of both primary cracks and damage on the bond in the inertia calculation. Murray et al. [9] considered historical loading to include the time-dependent effects of creep and shrinkage in Castel et al.'s method. Xu et al. [10] modified D ccc by using crack widths and creep coefficients to assess the time-dependent stiffness. e models of François et al. [6] and Castel et al. [7] were characterized by homogenized average inertias, and the inertias of cross sections near cracks were assumed to stay constant. So the strain distributions of these sections were also considered to be unchangeable. In fact, due to the crack presence, these sections do not remain plane, and their strain distributions show obvious nonlinearity [11,12]. Over a sufficient distance from a crack, Saint-Venant's principle indicates that the concrete cross sections should exhibit a gradual transition from highly nonlinear deformation to linear deformation [13]. erefore, the strain distributions of the sections change gradually near the cracks, and their neutral axes and inertias also vary in cracked regions. To exactly analyze the instantaneous stiffness of cracked reinforced concrete beams, the change in the strain distributions is taken into account in this paper. Key sections are selected in the region affected by cracks, and different nonlinear distributions are modeled for the key sections to consider the change. With the aid of internal force equilibria, top strains and neutral axes in the models are solved. Based on the solution, the section inertias can be exactly calculated to assess the instantaneous stiffness.

e Effect Region.
A reinforced concrete beam with distributed cracks is considered and assumed to be subjected to four-point bending, as shown in Figure 1. Due to the crack presence, concrete strain distribution is not linear along the height of cross sections, so a region possessing nonlinear strain distribution can be viewed as the region affected by the crack. e region affected by Crack I in Figure 1 is taken as an example to illustrate the region analysis. It is located around Crack I, but its length is still unknown.
To determine the length, the beam is assumed to be cut by Section C where Crack I lies and Section D where the external load is located. en, the free body of the beam CD is produced, as shown in Figure 2, and equivalent loads acting on the cuts are both equal to the moment M 0 , which is the product of the external load P 0 and its distance l p to the right support. e body CD is only subjected to the moment M 0 at its ends. According to Saint-Venant's principle [14], for each end, the concrete strains of sections near it are affected by the local distribution of the loads and do not agree with the plane-section hypothesis. e length of a local region where these effect sections are located is not bigger than the load length [15]. At Section C, the equivalent load is a resultant force of the stresses of concrete and steel bars, so the load length is not bigger than the beam height h. As a result, the length of the local region is also smaller than h. To simplify the analysis, the length is set as the section height h.
As the crack is located at one end of the body, the local distribution of the equivalent load at this end is directly related to the crack. As a result, the local region near this end affected by the load can be viewed as a region affected by the crack. Because the length of the local region is h, the length l ef of the effect region in the free body is also equal to h. e free body is on the right side of Crack I, so the right half of the entire region affected by the crack is determined. e entire region can be considered to be bilaterally symmetrical about the crack, so its length is 2l ef . For other cracks, the effect regions can be obtained by using the same method.

Models of Nonlinear Strain Distribution for Key Sections.
Concrete strain distributions change with their axial distances to the studied crack. For the section at the end of the effect region, the entire section bears the bending moment, and its strain distribution is linear along the section height. However, for the cracked section, only the concrete above the crack tip bears the moment, and its strain distribution exhibits high nonlinearity due to strain concentration near the tip, as shown in Figure 3.
To describe this change process, three types of key sections, denoted by Section 1, Section 2, and Section 3, are selected in the effect region, and different nonlinear strain distributions are modeled for these sections. Section 1 is the cracked section. Section 2 is closer to the crack compared to Section 3. e detailed characteristics of these sections will be studied as follows.

Strain Distribution of Section 1.
e crack tip appears in this section, and concrete strains near the tip are concentrated, so the strain distribution of this section is highly nonlinear. It is assumed that concrete strains above the neutral axis are linear to the height, while strains between the neutral axis and crack tip are proportional to the square root of their vertical distances to the neutral axis, as shown in Figure 3. So the concrete strain distribution can be written as where ε ct and ε cp denote the concrete strains at the section top and crack tip, respectively, and y n is the y-axis coordinate value of the neutral axis. In this section, concrete strains near the beam bottom are equal to zero. But the steel bar is subjected to a large tensile force and is debonding from the concrete. e steel strain is calculated according to the tension force transfer between the bar and concrete in an effective area around the bar. Wu and Gilbert [16] found that only concrete within a distance of 1.5d b from the bar surface influences the bondslip relationship, where d b denotes the bar diameter, so the depth d ea of the effective area is assumed to be equal to 3d b . e sum of the tensile forces of the bar and concrete in the effective area is considered to stay constant for the key sections and can be expressed as follows: 2 Advances in Materials Science and Engineering where ε cbs is the concrete strain in the effective area of Section 2, ε t1 and ε t2 denote the steel strains of Section 1 and Section 2, respectively, A ce is the effective area, and E t and E c are the elastic moduli of steel and concrete. So ε t1 can be calculated according to the following equation: e calculation of ε t2 will be introduced below.

Strain Distribution of Section 2.
Section 2 is close to the crack, and there is a slipping between the steel bar and concrete, so the steel strain is not equal to the concrete strain ε cl at the level of the bar. However, Manfredi and Pecce [17] considered that a section close to the crack possessed a characteristic that the concrete in compression and steel in tension could be strained according to the plane-section hypothesis, as shown in Figure 3. In this paper, Section 2 is considered to have such a characteristic, and its distance to the crack is assumed to be 0.2l ef . erefore, the steel strain ε t2 can be calculated based on the hypothesis as follows: e concrete strain distribution is still nonlinear along the section height. An inflection point appears at the position whose height is equal to the crack depth d c because the crack appearance causes distinct strain distributions above and below the crack tip. e bottom strain is small as this section is close to the crack, and the strains below the inflection point follow a nonlinear distribution, as shown in Figure 3. e concrete strain distribution can be expressed as where ε cb is the bottom strain.

Strain Distribution of Section 3.
In the three key sections, Section 3 is the farthest from the crack. e distance of Section 3 to the crack is considered to be a transfer length l t , which is an embedment length from the crack to the first point at which the strains of steel and concrete are equal to each other. According to fib Model Code 2010 [5], l t can be calculated by

Advances in Materials Science and Engineering
where c denotes the thickness of the concrete cover, ρ ef denotes the effective reinforcement ratio, f ctm is the mean tensile strength of concrete, and τ bms is the bond strength between steel and concrete; for stabilized cracking stages, τ bms � 1.8f ctm . erefore, the steel strain ε t3 of this section is equal to the concrete strain ε cl at the steel level. As Section 3 is the farthest from the crack, the crack effects on it are weakest. A bilinear distribution is used to simulate the concrete strain distribution, as shown in Figure 3. An inflection point still appears at the position whose height is equal to the crack depth d c . So the strain distribution can be expressed as e sections between two cracks can be considered to be only affected by the crack closer to them, and their strain distribution also can be analyzed by using the above method. However, key sections to be selected may change with the crack spacing l s . When l s is larger than 2l ef , the section at the end of the effect region can be taken as another key section. e strain distribution of this section is not affected by the cracks and can be obtained by using the plane-section hypothesis. When l s is between 2l t and 2l ef , the three types of key sections, Section 1, Section 2, and Section 3, should be considered. When l s is smaller than 2l t , the two types, Section 1 and Section 2, are considered.

Solution of Unknowns in the Models.
Although different models are built to describe the change in concrete strain distributions, there are still some unknowns in these models, such as the bottom strain ε cb , the top strain ε ct , and the neutral axis location y n . To solve so many unknowns, the bottom strain ε cb is needed to be modeled according to its change in the effect region and the internal force equilibria of the cross section are used.
For sections beyond the effect region, concrete strains accord with the plane-section hypothesis, and the bottom strain can be solved by using the classical beam theory as where I 0 denotes the inertia of transformed sections unaffected by cracks, including the contribution of steel. However, for the cracked section, its bottom strain is equal to 0, as only the upper part of the section bears the internal force, and there is no load acting on the crack interface. erefore, the bottom strains change from ε cb0 to 0 in the effect region. is change can be described by a quadratic polynomial [7]. e values of the bottom strains at the crack and the end of the effect region are equal to 0 and ε cb0 , respectively. Besides, the bottom strains at sections beyond the effect region are still ε cb0 , which indicates that the derivative of bottom strains dε cb /dx at the end is equal to 0. By substituting these boundary conditions into the quadratic polynomial, the bottom strains are expressed as where l c is the distance of the crack to the left support. e bottom strain is assumed to be symmetrical to the crack. As a result, two unknowns are left in the models of strain distribution: the top strain ε ct and the location y n of the neutral axis. en, the internal force equilibria of the cross section are used.
where A c is the concrete area. Based on equations (10a) and (10b), the unknowns will be calculated successively for the key sections. For Section 3, the bottom strain can be obtained from equation (9), and the steel strain ε t can be expressed in terms of concrete strains based on equation (7).
By substituting equations (7) and (11) into the force equilibrium equations (10a) and (10b), the following equations are obtained: where A cu and A cb denote concrete areas above and below the inflection point, respectively. Equations (12a) and (12b) are a set of simultaneous equations involving two variables. By solving equations (12a) and (12b), the two unknowns are determined, and the top strain and neutral axis are obtained. For Section 2, the steel strain can be expressed using the concrete strain, as shown in equation (4). By substituting equation (4) and equation (5)  e neutral axis of Section 1 can be considered to have the same height as that of Section 2 due to strain concentrations near the crack tip. From Figure 3, it is found that the concrete tensile strains are much bigger at Section 1 than at Section 2 because of the concentration effect, so the tensile force acting on the concrete at Section 1 is larger. Besides, the steel strain is also bigger at Section 1 according to equation (3). erefore, the total tensile force is larger at Section 1.
As the beam is subjected to the bending, the resultant force of compressed strains should be large to balance the  tensile force, and the compressed area is not smaller at Section 1. So the compressed area should not decrease near the crack, and the neutral axes of sections close to the crack are almost at the same height. is phenomenon was validated by numerical studies [11,12].
However, there is a new unknown in the distribution model of Section 3: the concrete strain ε cp at the crack tip. So two unknowns are still needed to be solved in the models: the concrete strains ε ct and ε cp . Equation (1) and equation (3) are substituted into equations (10a) and (10b), and an equation set similar to equations (12a) and (12b) is obtained. To solve this set, ε ct and ε cp are obtained.

Formulation of Stiffness Matrix.
In the effect region, concrete strains change linearly along the section height from the beam top to the neutral axis according to the strain distribution models in Figure 3. So the inertias I eq of studied sections can be expressed using the top strain and neutral axis as Using equation (13), the inertias are exactly obtained for the key sections. e variation curve of section inertias in the effect regions can be obtained by applying the piecewise polynomial fitting of the inertias obtained by equation (13), and the changing trend of inertias can be represented by the fitted curve. e finite-element method is used to analyze the stiffness characteristics of cracked concrete beams in this paper. Each element is viewed as a beam element with four degrees (node displacements at each end, transversal displacement, and rotation), and its local stiffness matrix can be written as where I f is the fitted curve of inertias along the element length, l me denotes the element length, ξ is equal to (x− x ml )/ l me , where x ml is the x-axis coordinate of the left node of the element, and N is the Hermite interpolation function [18]. e local stiffness matrix of each element is assembled as usual by the addition of stiffness coefficients corresponding to the same node, and then the global stiffness matrix of the cracked beam is obtained. is global matrix can be used to calculate the displacements of the beam subjected to a given load.

Experimental Program.
To test the proposed method, a comparison is made with experimental data from Castel et al. [8] on the instantaneous stiffness of five reinforced concrete beams labeled B1, B2, B3, B5, and B6. All the beams were simply supported beams with a span of 3.3 m, and their cross sections were 400 × 300 mm rectangular sections. e reinforcement layout of the beams is shown in Figure 4. e elastic modulus E c of concrete measured in standard conditions was 33 GPa, and the flexural tensile strength f ct f was 3.5 MPa. e yield stress f y of steel bars was 500 MPa.   All the beams were tested in four-point bending for all test cases, as shown in Figure 4. e test of the instantaneous stiffness was divided into two stages. e 1st stage test was carried out after a precracking loading test, which was done at 28 days after casting to establish a stabilized cracking pattern. In the first stage test, each beam was subjected to 10 loading and unloading cycles, and the cracking pattern and displacements were measured in each cycle. e cracking pattern is shown by the solid lines in Figure 5.
en, the beams were subjected to a sustained load for a period of 6 months. After 6 months, the cracking pattern, including the old cracks with a larger depth and new cracks that formed during the period of sustained loading, was measured, as shown in Figure 5 by the dotted line. e 2nd stage test was implemented, and the beams were then subjected to the same cycles of loading and unloading to again measure the instantaneous stiffness.

Validation and Discussion.
In the stiffness assessment, crack spacing is an important parameter. However, cracks in the tests were not vertical and had curved shapes. At different loading stages, cracks appeared with different spacing and depths. To objectively estimate the crack spacing, a method suggested by Gribniak et al. [19] is adopted in this paper. In this method, selected points in the cracks are mapped to the x-axis; then, a clustering technique is used to identify the points on the x-axis that are related to the same crack, and the average distance of the adjacent point cluster is taken as the crack spacing.
Besides, cracks do not show exactly the same pattern when observed from either side of the beams, so the crack spacing should be estimated from both the south and north sides. e inertias along the beam span also should be assessed from the two sides. e beam B1 is taken as an example to illustrate the assessment process. e measured depths and estimated locations of cracks in the beam B1 are shown in Tables 1 and  2. According to these depths and locations, the key sections in the region affected by cracks are determined, as shown in Table 3, and their nonlinear strain distributions are modeled using the proposed method. By solving these models, the concrete strains and neutral axis are obtained, and the inertias of the key sections are calculated by using equation (13). en by applying the piecewise polynomial fitting of these inertias, the variation curves of inertias I f s and I f n are estimated for the north and south sides. Finally, the average of the two fitted inertia curves is taken as the inertia of the beam B1, as shown in Figure 6.
From Figure 6, it is observed that the strains do not remain constant near the cracks. e closer the studied section is to the cracks, the smaller its bottom strain becomes but the larger its top strain gets. e top strains reach their local maximum at every crack. Besides, the maximums are e neutral axis also changes continuously along the beam span and becomes higher as it gets closer to the cracks.
As the strains and neutral axis change gradually in the effect regions, the inertias I f s and I f n and their average vary continuously according to equation (13). e inertias are compared with the homogenized inertias calculated by Xu et al. [10] in Figure 6(e), which were obtained by assuming that inertias and strain distributions were constant near cracks. In fact, more and more studies show that the strain distributions near cracks change gradually [11][12][13], so the changing inertias are more suitable for cracked beams.
From Figure 6(e), it is also seen that the beam inertias are smaller at the 2nd stage than at the 1st stage. e larger depth of the old cracks at the 2nd stage is considered to be the main reason. Although there are new cracks that formed during the period of sustained loading, their number and depth are both small. So the effect of the new cracks on the inertias is small. However, the larger depth of the old cracks at the 2nd stage causes the bigger top strains and higher neutral axis, as shown in Figures 6(b) and 6(c). As a result, the inertias at the 2nd stage are smaller. Similar phenomena also can be observed in the inertias of the other beams shown in Figure 7, which are also calculated through the proposed method.
Based on these inertias, the stiffness matrices of these cracked beams are estimated by using equation (14) and adopted to calculate the beam displacements. ese calculated results are compared with the data measured in the experiments, as shown in Figure 8. It is found that the calculated results provide a good agreement with the measured data for both test stages. If an equivalent stiffness is considered to be equal to the load divided by the displacement, it is evaluated by using the experimental data and the calculated displacements, respectively. e evaluated results are compared with the stiffness obtained by Xu et al. [10], as shown in Table 4.
From Table 4, it is seen that, at the 2nd stage, the difference ratios of the experimental and calculated stiffness are not more than 4%, while the ratios at the 1st stage increase but still are not more than 7%, which indicates the exactness of the proposed method for assessing the stiffness of reinforced concrete beams.
For the beams B5 and B6, the stiffness obtained by the proposed method is closer to the experimental stiffness, compared to the results predicted by Xu et al. [10]. It is because Xu et al. [10] assessed the stiffness using the coefficients of shrinkage and creep, but these coefficients had an error for the beams B5 and B6 [10]. On the contrary, the instantaneous stiffness is taken as a deformation resistance at an instant, and the shrinkage and creep are long-term deformations, so the shrinkage and creep could not take their effects at an instant, and their effects are not considered in this paper. erefore, the results obtained by the proposed method are not disturbed by the error.   8 Advances in Materials Science and Engineering For the same beam, the stiffness is smaller at the 2nd stage than at the 1st stage, which coincides with the smaller inertia at the 2nd stage shown in Figure 7. As mentioned earlier, the larger crack depth at the 2nd stage is the main reason for the smaller inertia, so it also causes the smaller stiffness.
e stiffness of the beam B3 is smaller compared to that of the beam B1. It is due to the fact that the beam B3 has more cracks and its entire region affected by the cracks is longer. As the inertias of sections in the effect region decrease, the stiffness of the beam B3 becomes smaller. Similar phenomena also can be found by comparing the stiffness of the beams B2 and B6. erefore, the instantaneous stiffness of cracked reinforced concrete beams greatly depends on the crack depth and number. e proposed method models different nonlinear strain distributions for key sections to consider the gradual change in strain distributions caused by cracks, and its calculated results can exactly reflect the dependence of the beam stiffness on the cracking pattern, including the crack depth and number.

Conclusions
e following conclusions can be drawn from this paper: (1) A method for assessing the instantaneous stiffness of reinforced concrete beams is proposed based on the gradual change in concrete strain distributions near cracks, which is considered through modeling   different strain distributions for key sections in the regions affected by the cracks. With the aid of internal force equilibria, the top strain and neutral axis are solved. Based on the solution, the section inertias can be exactly calculated to assess the beam stiffness. (2) e effectiveness of the proposed method has been validated experimentally using five reinforced concrete beams with cracks. At the 1st and 2nd stage tests, the difference ratios of the experimental and estimated stiffness are not more than 7% and 4%, respectively. (3) e instantaneous stiffness is taken as a deformation resistance at an instant, and its value is proven to greatly depend on the crack depth and number. e proposed method can exactly reflect this dependence through simulating the change in concrete strain distributions near cracks. (4) e proposed method only utilizes the information of the current cracking pattern, including the crack number and depth, and does not require the knowledge on the previous load history, shrinkages, and creeps. So it is applicable to assess the instantaneous stiffness of existing cracked reinforced concrete beams under stabilized cracking stages.
Data Availability e data analyzed during the current study include the responses of cracked reinforced concrete beams. ey were derived from the following public domain resource: https:// doi.org/10.5281/zenodo.3244472.

Conflicts of Interest
e authors declare that they have no conflicts of interest.