Abstract
The crack presence causes nonlinear stress distributions along the sections of a beam, which change the neutral axis of the sections and further affect the beam stiffness. Thus, this paper presents a method for the stiffness estimation of cracked beams based on the stress distributions. First, regions whose stresses are affected by the crack are analyzed, and according to the distance to the crack, different nonlinear stress distributions are modeled for the effect regions. The inertia moments of section are evaluated by substituting these stress distributions into the internal force equilibrium of section. Then the finite-element technique is adopted to estimate the stiffness of the cracked beam. The estimated stiffness is used to predict the displacements of simply supported beams with a crack, and the results show that both static and vibrational displacements are accurately predicted, which indicates that the estimated stiffness is precise enough. Besides, as the section shape of beam is not limited in the process of modeling the stress distributions, the method could be applicable not only to the stiffness estimation of cracked beams with a rectangular section, but also to that of the beams with a T-shaped section if the crack depth ratio is not larger than 0.7.
1. Introduction
Cracks tend to appear in beam-like structures owing to overloading, improper construction, temperature effects, and other causes. The presence of a crack induces changes in the structural stiffness and may result in adverse effects on the static and dynamic behaviors of the beam. So in the last four decades, numerous researchers have paid considerable attention to cracked beams [1], and different methods are reported for crack modeling in the beams to describe the changes in the beam stiffness. A majority of the considered methods can be attributed to one of the following categories: rotational spring model and continuous flexibility model.
In numerous works including those using crack models with rotational springs, the cracked beam was treated as two undamaged beams connected by a rotational spring at the cracked section. The stiffness of the spring was related to the crack depth and section geometry, and its value was estimated by using fracture mechanics methods to quantify the relation between the applied load and stress concentration around the tip of the crack [2–7]. This spring caused the rotation discontinuity owing to the concentrated flexibility and produced an additional relative rotation of the cracked section.
Based on this rotational spring model, the behaviors of cracked beams were studied. Mahmoud and Abou Zaid [8] adopted the spring model for a crack and developed an iterative modal approach to analyze the vibration of cracked beams subjected to moving loads. Jaksic et al. [9] deduced the eigenvalues of a cracked beam based on the assumption that the spring stiffness was a polynomial function of crack depth ratio and formulated the motion equations of the beam–vehicle interaction. By modeling the crack as a spring, Ariaei et al. [10] used the finite-element approach to analyze the responses of cracked beams under moving masses.
Various continuous flexibility models were also proposed for modeling cracks to develop structural vibration equations. Christides and Barr [11] firstly assumed that the stress around a crack decayed exponentially with the distance from the crack and exhibited a stiffness reduction in the region near the crack tip with an exponential variation. However, Sinha et al. [12] proposed a stiffness reduction with a local effect governed by a triangular variation. Bilello [13] analyzed the effect of a crack in terms of the ineffective area delimited by a linear reduction of its height, starting from the cracked section. Chondros et al. [14] modeled a crack as continuous flexibility by using the displacement field in the vicinity of the crack tip, found with fracture mechanics methods. Yang et al. [15] computed the equivalent bending stiffness by using the strain energy variation around the crack and considered the cracked beam as a continuous system with varying moments of inertia. Abdel Wahab et al. [16] described the stiffness reduction by using a cosine function with three parameters that were determined by the vibration characteristics extracted from the experimental data.
These continuous flexibility models were adopted to analyze the behaviors of cracked beams. Law et al. [17] used the continuous model proposed by Abdel Wahab et al. [16] to simulate cracks in concrete simply supported beam and adopted the modal analysis approach to build an interaction model between the cracked beam and a moving vehicle. Fu [18] adopted the same approach to analyze the dynamic responses of a continuous cracked beam under moving vehicles. Chondros et al. [19] employed the continuous flexibility model [14] to predict the changes in vibration frequencies of a simply supported beam with a breathing crack. Caddemi et al. [20] modeled the crack influence by means of generalized functions and studied the nonlinear dynamic response of a beam with switching cracks.
From the crack models used in the literature, it is found that all the models considered the local flexibility due to a crack appearance and most of them assessed the flexibility magnitudes either by experiment or by fracture mechanics methods [2, 16]. If experiments are desired, then the responses of cracked beams cannot be predicted only by theoretical calculations [16]. When fracture mechanics methods are adopted, the cracked beam should have a rectangular section [2, 16]. This is because in fracture mechanics methods, the beam is considered as a two-dimensional plane with the stress concentration being focused near the crack tip.
In fact, the crack presence does not affect only the stress of the vicinity of the crack tip but also the stress of regions adjacent to the vicinity. As the neutral axis of section is directly related to the stress distributions along the section height, there is a change in the neutral axis, which also can be observed from two-dimensional numerical simulations of cracked beams [21, 22]. So the inertia moments of the sections and the stiffness of these regions also vary with the stress distributions.
Meanwhile, the length of these regions may be larger than that of the vicinity, so the stiffness of these regions will influence the global stiffness of the beam, and it is necessary to estimate the global stiffness from the perspective of the stress distributions of the regions. Besides, in the estimation of stress distributions, the section shape of the beam is not limited, and the analysis may be applicable for cracked beams with irregular sections such as T-shaped sections.
In this paper, the nonlinear stress distributions of the beam are modeled. Then, with the aid of the internal force equilibria, the inertia moments of the sections are estimated, whereas the local stiffness matrices are determined by using the finite-element technique. Finally, these matrices are assembled to form the global stiffness of the beam.
2. Stiffness Formulation for a Cracked Beam
A simply supported beam with an open transverse crack is shown in Figure 1(a). The axial direction of the beam is taken as the axis and the vertical direction as the axis. It is assumed that the beam is homogeneous and cross-section is uniform along the axis. The two supports are denoted as Support A and Support B, and a crack with a depth of is located at the bottom of Section C.

(a)

(b)
2.1. Models of Stress Distributions Near the Crack
Suppose that a couple of static moments act on the two ends of the beam, such that both of their values are equal to and the action length of each moment is equal to the height of the section at the end. Next, the beam is subjected to pure moments according to the mechanical characteristics of simply supported beams.
If the beam is cut by a cracked section perpendicular to the axial direction of the beam, the free body of the beam AC is produced as shown in Figure 1(b). The equivalent moment of force acting on the cut is , but its action length is equal to the height () of the cracked section.
The stress of the free body should be firstly analyzed to acquire the stress of the beam. The body is only subjected to the pure moments at its ends, and the solution problem of its stress can be viewed as Saint–Venant’s problem of pure bending. The stress in the part near the loaded ends is affected by the local load [23]. According to the existing findings, the effect region is not longer than the loading action [24]. For a pure bending beam, the action length is equal to the height of the section. Therefore, the length of the effect region is set as the section height .
The stress in the effect region will be investigated in detail. For a bending beam, the stress in the axial direction is strongly related to the load, and it changes both along the axis and cross-section of beam. Therefore, the axial stress is selected as the main study object, and to simplify the analysis, it is split into stresses along the axis of the beam and the cross-section. These two types of stress will be analyzed.
First, the stress along the axis direction is studied, and the stress of the beam bottom is taken as the studied object. For Section C in Figure 1(b), this stress is equal to 0, because the equivalent moment acts on the upper part of the section and there is no load near the beam bottom, as shown in Figure 1(b). However, for Section D in Figure 1(b) which is just located at the end of the effect region, the stress is no longer affected by the loaded end and its value can be obtained using the classical beam theory. where denotes the inertia moment of the uncracked section and is the vertical distance from the bottom to the neutral axis. Therefore, the bottom stress varies gradually from Section C to Section D.
In the proof of Saint–Venant’s principle, Toupin found that the stress near the loaded end changed exponentially with the distance from the end [25]. Referring to this finding and taking the effects of the crack depth into account, this study assumes the stress of the beam bottom in the effect region is expressed aswhere denotes the axial distance from the studied section to Section C, denotes the change rate, and is set as 6; () denotes the crack effect on the change rate, where the depth ratio of the crack is equal to , and the constant is set as 0.5.
Next, the stress along the cross-section of the beam is analyzed. This stress is assumed to be constant along the width of the section, and it changes only along the section height. As Section D is sufficiently far from the loaded end, its stress is not affected by the local load and, according to the classical beam theory, can be expressed as [26]. where is the vertical coordinate of the neutral axis. So the stress of Section D is linearly related to the height.
However, for sections closer to Section C, the stress along the height exhibits a nonlinear characteristic, which becomes more obvious with the decrease in the distance to Section C owing to the local load or, more precisely, the crack effects. To simulate the variable nonlinear characteristics, three models of the stress distribution are built.
2.1.1. Bilinear Distribution Model
The crack opens in the lower part of Section C, but the upper part is intact. This behavior will affect the stress distribution in the studied section and cause different distributions to appear in the lower and upper parts of the section. To simulate this difference, a bilinear distribution, the simplest nonlinear distribution, is used. The stresses in both the parts are linear, but the slopes are distinct. There is an inflection point between these parts, whose -axis coordinate is set to be equal to , as shown in Figure 2(a). The stress along the height is expressed aswhere denotes the axial stress of the beam top.

(a)

(b)

(c)
The nonlinearity in the bilinear distribution model is slight, and the stress is moderately affected by the crack. Therefore, this model is applicable to the sections that are in the effect region but closer to Section D rather than Section C.
2.1.2. Curve Distribution Model
With decreasing the distance to the crack, the stress of the beam bottom reduces according to (2), but the stress of the inflection point increases. If the former is higher than the latter, the bilinear distribution model can be applicable; otherwise, the nonlinearity is considered, and for the part of the section below the inflection point, the change in the stress along the height follows a quadratic curve, as shown in Figure 2(b). The stress along the height is expressed as
Crack presence obviously affects the stress distribution, particularly the distribution on the part below the inflection point. Therefore, this model is applicable to the sections near the middle of the effect region.
2.1.3. Highly Nonlinear Distribution Model
When the distance to Section C is less than 0.2 , the section is regarded to be in the vicinity of the crack tip and the stress distribution becomes highly nonlinear because the crack tip exhibits the stress concentration phenomenon.
As the stress near the tip is concentrated, its value is much higher than that of the top stress, so that the stress of the upper part of the section becomes nonlinear to the section height. It is assumed that the stress distribution on this part is bilinear and the inflection point is located at the neutral axis of the section.
For Section C, the stress of the part below the crack tip is equal to 0, because there is no force on the crack surface. The closer the studied section to Section C, the larger the part whose stress is approximately equal to 0. The height of this part is assumed to linearly vary along the axial direction of the beam. At Section C, the height is equal to be the crack depth . At the section whose distance to Section C is 0.2 , the height is equal to 0.
Concurrently, the concentrated stresses are tensile, and the resultant force of these stresses is large. As the beam is subjected to pure moments, to balance this large tensile force, the resultant force of the compressed stress should be large and the compared area of the section should not be small. In this model, the compressed area is assumed to be constant and set as the compressed area of the section whose distance to Section C is 0.2 . Consequently, the neutral axis of the section is also kept fixed in this model.
The stress in this model based on the above assumptions is shown in Figure 2(c) and can be expressed aswhere denotes the stress of the inflection point whose -axis coordinate is .
Although the above stress distributions are assumed, their validity will be discussed in Section 3. In addition, these assumptions do not limit the shape of the beam cross-section, and the stress distributions are valid irrespective of the shape.
2.2. Calculation of Inertia Moments
The inertia moments of the effect region cannot be obtained directly based on these distribution models, as the stress of the beam top and neutral axis of the section are still unknown. To obtain the top stress and neutral axis, the internal force equilibrium of the section should be adopted. As the studied body is subjected to pure bending, the following force equilibria are satisfied for each section:where denotes the area of the studied section. Using these force equilibriums, the top stress and neutral axis will be calculated successively for the sections in the macro elements.
First, the sections whose distances to Section C are not less than 0.2 are analyzed. The stresses of these sections follow a bilinear or curve distribution. The bilinear distribution model is initially used to simulate the stresses and substituted into (7a) and (7b) so that the following equations are obtained:where and denote the areas of the parts above and below the inflection point, respectively, in the bilinear distribution model.
The bottom stress in (8a) and (8b) can be obtained considering (2); however, there are still two unknowns: the top stress and the coordinate of the neutral axis . Therefore, (8a) and (8b) form a set of nonlinear equations with two unknowns, and Newton’s method can be used to solve this set. After obtaining the solution, these unknowns can be determined and the stress at the inflection point can be obtained as
According to the assumptions of the stress distribution models, if is lower than the bottom stress , the bilinear distribution model is applicable for the studied section and the two unknowns, the top stress and the coordinate of the neutral axis , are obtained. Otherwise, the curve distribution model is adopted and and should be recalculated as follows.
By substituting (5), the stress expression of the curve distribution model, (7a) and (7b), a nonlinear equation set similar to (8a) and (8b) is obtained. Solving this set yields and .
Next, the sections whose distances to Section C are less than 0.2 are studied. The highly nonlinear distribution model is applicable to these sections. In this model, the compressed area of each section is assumed to be constant and the -axis coordinates of the neutral axis are equal to the coordinate for the section whose distance to the crack is equal to 0.2, which is calculated by using the abovementioned method.
Although the neutral axis is known, there are still two unknowns in this model: the top stress and the stress of the inflection point. By substituting (6), the stress expression of the highly nonlinear distribution model, into (7a) and (7b), nonlinear equations similar to (8a) and (8b) are obtained. Solving these equations yields the stresses and .
In Section D, the stress is linear to the height, as expressed in (3). If the top stress and neutral axis are known, then based on (3) the inertia moment can be expressed as
In other sections of the effect region, the stress is also linear from the beam top to the neutral axis according to the three stress distribution models. In addition, the top stress and neutral axis are determined by solving the nonlinear equations. Similar to the inertia moment of Section D, the inertia moment in the effect region can be expressed using the known top stress and neutral axis, and its formulation is the same as (10).
Although the inertia moment of each section in the effect region can be obtained by using the above method, it is unnecessary to calculate the moments of all the sections, because the variation curve of the moments along the region length can be fitted using the moments of some given sections. If the number of the given sections is sufficient, the fitted curve can represent the real change in the moments.
2.3. Formulation of Stiffness Matrix
If the finite-element method is adopted to analyze the cracked beam, the structure can be divided into beam elements, including () common elements and two cracked elements, as shown in Figure 3. Each element has four degrees (node displacements at each end, vertical displacement, and rotation). The common element can be viewed as a beam element with constant cross sections. So its stiffness is not affected by the crack and easily obtained by a variational method [27].

However, the cracked elements are located near the crack, and their stiffness is significantly affected by the crack. For example, the region from Section D to Section C can be taken as a cracked element with variable sections, whose local stiffness matrix is expressed as [27] where denotes the elastic modulus of the beam material, is the fitted curve of the inertia moments along the element length, is equal to , where is the -axis coordinate of the left node of the macro element, and is the Hermite interpolation function, expressed as [27]
The global stiffness matrices for cracked elements are developed by using transformation matrices and then assembled to form the stiffness matrix of the beam, which can be written as where is the stiffness matrix of the common elements, denotes the stiffness matrix of the th cracked element, is the transformation matrix of the th cracked element, and denotes the serial number of the element. can be written as
When the beam is subjected to dynamic loads, the governing equations of the beam can be expressed as where and are, respectively, the mass and damp matrices of the beam and denotes the vector of the moving load. Explicit expressions for these characteristic matrices and the vector are provided in [27]. is the displacement vector of the beam, where and are the vertical displacement and sectional rotation at the section centroid of the th node, respectively.
When the beam is under static loads, the inertial and damping forces vanish in the governing equations. By solving these equations, the static or dynamic displacements of the beam can be obtained.
3. Examples
3.1. Cracked Beams with a Rectangular Section
To test the proposed method, the cracked beam considered in this study is a simply supported beam with a rectangular section used by Mahmoud et al. [8]. The following beam and crack parameters are adopted: L = 50 m, E = 2.1 × 1011 N m−2, m = 3930 kg m−1, h = 1 m, b = 0.5 m, = 0.5 m, and the distance from the crack to the left-hand support = 25 m. The beam is divided into 26 elements, including 2 cracked elements and 24 common beam elements. The intersection of the cracked elements is the cracked section.
The beam finite-element method (BFEM) proposed in this paper is adopted to calculate the responses of the beam, and the obtained results are compared with the results of the following two methods.
The first one is the solid finite-element method (SFEM) using the commercial finite-element program ANSYS. The beam is simulated with three-dimensional solid elements, and the crack in the program is taken as a shot which is formed by disjointing the nodes of the two adjacent elements at the crack interface. The element size is set to 0.1 m, and the smaller element size requirements near the crack tip are provided by the refined meshing technique, as shown in Figure 4.

The second approach is an iterative modal analysis approach (IMAA) proposed by Mahmoud et al. [8], where the crack is modeled as a rotational spring connecting two undamaged beam segments. The stiffness of the rotational spring is determined using fracture mechanics and derived using the results of a beam with a rectangular section.
To obtain the inertia moments of the sections in the cracked elements, a couple of static moments with a value of 1226 kN m are imposed on the ends of the beam to analyze the stress distribution in the cracked elements. According to the analysis in Section 2.2, the length of a cracked element is equal to the height of the beam, which is equal to 1 m. Then, thirteen sections are selected, which divide the element into 12 segments and the stresses of these sections are analyzed.
The bottom stresses of the sections are estimated using (2), and the top stresses and the neutral axis are calculated based on the stress distribution models and force equilibria. These stress results are compared with the results obtained by the SFEM, as shown in Figure 5.

(a)

(b)
From Figure 5, it is observed that the two results almost completely overlap, except in the small region close to the crack. This indicates that the assumed stress distributions in this study are in accordance with the solid finite-element simulations for the sections that are not close to the crack. Although there is a difference between the two results in the region close to the crack because it is difficult to accurately model the stress concentration phenomenon, the maximum difference is less than 4% of the stress peak and the length of the small region is less than 4% of the length of the cracked elements. Because this difference is small, it may have a slight effect on the calculation of the inertia moments.
Based on the calculated stresses and neutral axis, the inertia moments of the selected sections are estimated considering (10) and used to fit the curve of the inertia moment along the element length, as shown in Figure 6. It is found that the inertia moment reaches its minimum value at the crack and the farther away the section is from the crack, the larger the inertia moment is. The inertia moment of the farthest section in the cracked element is equal to the moment of the uncracked section.

The local stiffness matrices of the cracked elements are calculated by substituting the fitted inertia moments into (11), and then the global stiffness matrices are obtained using the transformation matrix presented in (13). Finally, the stiffness matrix of the beam is formed by combining the global stiffness matrices of the cracked elements and common beam elements.
Suppose that a vertical load with a value of 385.5 kN statically acts on the beam and it is located at the midspan of the beam. The static displacement of the beam is calculated by the proposed method and compared with the two results obtained by the SFEM and IMAA, as shown in Figure 7. In this paper, all displacements are normalized relative to the value , which is the static displacement of the beam with no crack due to at the midspan. From Figure 6, it is observed that the three results are almost coincident. This agreement indicates that the proposed method can yield an accurate prediction of the static responses of the cracked beam.

Then suppose that the load moves into the beam from the left-hand support at time and its velocity is 40 m s−1. The vibrational displacements of the beam obtained by different methods are compared, and their results are shown in Figure 8. In this figure, is the total time that the load requires for one pass the beam. It is seen that the results from different methods coincide with each other, so the vibrational responses are accurately calculated.

These accurate results of the static and vibrational displacements indicate that the stiffness of the cracked beam estimated by this paper is precise enough to predict the displacements and accurately describes the change in the beam characteristics caused by the crack.
Figure 9 investigates the effects of the crack depth on the beam stiffness. It is seen that all the displacements obtained by the three methods increase with the depth and are almost coincident when the depth ratio is less than 0.6. When the ratio is larger than 0.6, the displacement obtained by the IMAA is smaller than the result simulated by the SFEM. However, the result obtained by the BFEM is identical to the simulated result even if the ratio reaches 0.8. Therefore, it is concluded that the proposed method can accurately estimate the stiffness of the cracked beam even when the crack is deep.

(a)

(b)
The effect of crack location on the displacement is investigated in Figure 10. It seems that when the crack is not in the middle of the beam, the displacement predicted by the BFEM still corresponds with the result simulated by the SFEM, which indicates that the crack location does not affect the accuracy of the predicted displacement.

3.2. Cracked Beams with a T-Shaped Section
A cracked beam with a T-shaped section considered here is a simply supported beam, whose crack and beam parameters are the same as those of the beam with a rectangular section mentioned in Section 3.1, except the cross-section. The T-shaped section is shown in Figure 11, and its height is identical to the height of the rectangular section. The crack depth is set as 0.5 m, and the beam is also divided into 26 elements, including two cracked elements.

To calculate the inertia moment along the length of the cracked elements, two static moments of 1720 kN m are imposed on the ends of the beam and the stress distribution of the cracked elements is analyzed. According to the analysis in Section 2.2, the element length is equal to the height of the T-shaped section, which is equal to 1 m. Then, 13 sections are selected, and the stresses of these sections are analyzed.
The bottom stresses of the sections are estimated using (2) and, based on the stress distribution models and force equilibria, the top stress and neutral axis are calculated. These stress results are shown in Figure 12. It is seen that the stresses obtained by the proposed approach and the SFEM are similar both at the bottom and top of the sections. This proves the validity of the stress distribution models used in this study.

(a)

(b)
Based on the calculated stresses and neutral axis, the inertia moments of the selected sections are estimated based on (10) and used to fit the curve of the inertia moment along the element length, as shown in Figure 13. The trend of the fitted curve is similar to that for the beam with a rectangular section.

The fitted inertia moments are used for the formation of the stiffness matrix of the beam subjected to the load mentioned in Section 3.1. When the load statically acts on the midspan, the static displacements are shown in Figure 14 and it is observed that the displacements obtained by the BFEM are in good agreement with the results simulated by the SFEM, but there are some differences between the displacements obtained by the IMAA and the simulated results.

When the load moves into the beam at a velocity of 40 m s−1, the vibrational displacements of the beam are shown in Figure 15. It is seen that the vibrational displacements obtained by the BFEM are closer to the results simulated by the SFEM, compared with the displacements obtained by the IMAA.

These phenomena of static and vibrational displacements are attributed to the different assumptions used for these methods. The IMAA derives the stiffness of the cracked region using fracture mechanics and assumes that the beam section was rectangular, whereas the BFEM assumes that the stress distribution of the region and does not limit the shape of the cross-section. As the beam section here is T-shaped, the BFEM is more applicable and will be used for the following study.
Figure 16 illustrates the effects of the crack depth on the displacements of the beam with a T-shaped section. It is seen that the displacements obtained by the BFEM agree with the simulated results if the ratio of the crack depth to the section height is not larger than 0.7. When the ratio reaches 0.8, the stress concentration in the beam becomes extremely high, so that the stress distributions along the section height could not be simulated by the proposed models. Therefore, it is difficult to accurately calculate the stiffness of the sections near the crack and, consequently, the response of the beam using the BFEM.

(a)

(b)
4. Conclusions
In this paper, a method for the stiffness estimation of cracked beams is proposed based on nonlinear stress distributions near the crack. First, regions whose stresses were affected by the crack were analyzed and, according to the distance from the studied sections to the crack, different nonlinear stress distributions were modeled for the effect regions. The inertia moments of the sections were evaluated by substituting the stress distributions into the internal force equilibrium of the sections. Then, the finite-element technique was adopted to estimate the stiffness matrix of the cracked beam.
The estimated stiffness was used to predict the displacements of cracked beams, and the results show that both static and vibrational displacements were accurately predicted, which indicated that the estimated stiffness was precise enough. Besides, the method could be applicable not only to the stiffness estimation of cracked beams with a rectangular section, but also to that of the beams with a T-shaped section if the crack depth ratio was not larger than 0.7, as the section shape of the beam was not limited in the process of modeling the stress distributions.
In practical applications, the proposed method is attractive for use, as the beam stiffness is established by using the nonlinear stress distributions, in place of fracture mechanics, and the computation becomes simple. Although this method analyzes the beam with a crack, it still works for beams with multiple cracks. But if the cracks are closely located, the stresses of the region in a crack are affected by other cracks, and they do not follow the models of nonlinear stress distributions in this paper. So this method may not be applicable for beams with closely located cracks, whose stress distributions and stiffness should be studied in futures works.
Data Availability
The data analyzed during the current study include the responses and stiffness characteristics of cracked beams. They were derived from the following public domain resources: https://pan.baidu.com/s/1ktwUqubg9ycuXnsLmtRcKg.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is sponsored by the National Natural Science Foundation of China (51508155) and the Fundamental Research Funds for the Center Universities (2017B12314).