Approximate analytical solution of buckling of multi-damaged column-like structures using a continuous diffused crack model by variational iteration method

Compressive structural members can be locally damaged by overloading, corrosion, car crash and fire. In this work, a continuous diffused crack model is proposed to study the static stability of Euler–Bernoulli rectangular column-like structures under different boundary conditions. The governing differential equation is formulated by adopting a diffused crack model. The powerful variational iteration method is implemented to find the approximate analytical buckling modes and buckling loads based on the buckling response of the intact column. A novel generalized Lagrange multiplier is derived. The proposed method incorporates the effects of the crack width into consideration when deriving the buckling modes. The stability equation allows addressing the influences of multiple damages and can be applied to both concentrated and distributed cracks. The famous Rayleigh–Ritz method is utilized to verify the computed buckling loads. The proposed diffused crack model and the application with VIM is efficient and accurate for handling buckling problems of cracked columns under different boundary conditions.


Introduction
Euler-Bernoulli columns play an important role in civil, mechanical, aeronautical, nuclear and offshore structures and the related static, dynamic and stability analysis is an important topic for column design and maintenance. Damages to the column-like structures may cause catastrophic disasters, which is one of the concerns of structural health monitoring [1] and it is crucial to determine the buckling loads of locally weakened columns. In engineering practice, local damages are frequently associated with cracks [2,3] and the influences of cracks on the buckling of columns is a fundamental research topic of stability analysis. Extensive studies on the buckling of columns have been conducted and there are several crack models for studying the stability of cracked columns, namely, the discrete spring crack model, the singularity crack model, the stiffness reduction model and the finite element model.
In the discrete spring crack model [4], the column is discretized and separated by rotational springs at each of the cracked cross-sections. To analyze the mechanical behavior of cracked columns, the cracked cross-section is simulated by a rotational spring whose rotational stiffness is determined by the crack size and the crack geometry as proposed by Freund and Herrmann [5]. In this model, the continuity of deflections, bending moment and shear force at each of the cracked cross-sections are required while the slope jump is allowed [6]. The theory of fracture mechanics is utilized to link the rotational stiffness of the spring to the size and geometry of a crack [7]. The discrete spring model was adopted by Wang et al. [8] to analytically investigate the buckling of internally weaken columns. Recently, Planinc and Schnabl [9] applied the discrete spring model to study the buckling loads of twolayer columns with transverse cracks and partial delamination analytically.
In the singularity crack model, the flexural stiffness reduction due to a concentrated crack is simulated by the Dirac's delta function. The singularity crack model can be applied to simulate the rotational discontinuity by a macroscopic approach, which is equivalent to a rotational spring [10]. The representative works were carried out by Caddemi et al. [11][12][13][14][15]. Biondi and Caddemi [14,15] applied the singularity crack model to study the Euler-Bernoulli beams with concentrated cracks and concluded that the singularity crack model is equivalent to the discrete spring model. In a follow-up study, Caddemi and Calio [11] investigated the buckling load of a Euler-Bernoulli column with multiple cracks using the Dirac's delta function. The singularity crack model was also used to derive the exact dynamic stiffness matrix of a Euler-Bernoulli column [16] with multiple concentrated cracks. The exact solution of the buckling load and buckling mode was obtained. The singularity crack model was adopted to study a Timoshenko beam with multiple cracks by Caddemi et al. [13]. The deflection and rotation discontinuities of the beam with multiple external concentrated supports were investigated under static load. Caddemi et al. [12] addressed the effects of elastic end support on the dynamic stability of Beck's column in the presence of weak sections by the singularity model. The exact solution of the buckling modes and buckling load of the weakened cantilever column was obtained.
When the length of the stiffness reduction zone is relatively small compared with the total span of the beam, the discrete spring crack model and the singularity model is acceptable, since the stiffness reduction near the crack is neglected [17]. Besides, it is inconvenient to relate the crack depth to the stiffness of the rotational spring directly [18].
As was noticed by scholars, the flexural stiffness near the vicinity of a crack is reducing gradually and several diffused stiffness reduction crack models were proposed to consider this effect [18,19]. Assuming that the stress decays exponentially from the tip of the crack within the length of stiffness reduction zone, Christides and Barr [19] studied the dynamic response of beams with pairs of symmetrical open cracks on the top and bottom. An exponential function was proposed by Christides and Barr [19] to account for the flexural stiffness reduction. Based on the crack model of Christides and Barr [19], Wahab et al. [20] proposed a stiffness reduction function in terms of the cosine function and the reduction in elastic modulus was considered to represent the stiffness loss near the crack. Based on the studies of Christides and Barr [19] and Wahab et al. [20], Sinha et al. [18] proposed a simplified form of stiffness reduction in the vicinity of a crack, in which the flexibility varies linearly from the cracked section to the intact beam section. According to this model, the stiffness reduction length on one side of the crack is l c .
The Christides and Barr crack model, Wahab crack model and Sinha crack model are only piece-wise continuous in terms of the stiffness of the total length of the beam. It is inconvenient to work with them to obtain an analytical solution of the four-order differential equation that governs the buckling of a column. Most of the crack models in the literature can be applied in a finite element model, the interested readers can refer to the references [21] for more information.
In the practice of civil engineering, the approximate analytical solution is desirable. In the present work, closedform expressions of buckling loads and buckling mode shapes are derived and illustrated for Euler-Bernoulli columns with multiple cracks by applying the diffused stiffness reduction crack model. In this model, a crack is described by three parameters: the crack ratio, crack position and nominal crack width. The diffused crack model can simulate concentrated cracks and local damage with a certain length, such as local delamination and can accommodate multiple cracks without increasing the rank of the coefficient matrix. The powerful variational iteration method is implemented to derive the buckling loads and buckling mode shape and a new generalized Lagrange multiplier is also derived.
The rest of the paper is organized as follows: The differential equation governing the buckling of a cracked column is reformulated with a diffused crack model in Sect. 2. In Sect. 3, the buckling mode solution of an intact column is rederived with a novel Lagrange Multiplier. With the obtained solution in Sect. 3, the buckling solution of a multi-cracked column is developed in Sect. 4. In Sect. 5, the Rayleigh-Ritz method for a cracked column is introduced. Numerical examples are illustrated and discussed in Sect. 6 and conclusions are made in Sect. 7.

Diffused crack model
Assuming the approximate remaining area of the locally damaged column follows the function where A 0 is the area of the intact cross-section, ci is the crack ratio of the ith crack, x 0i is the position of the ith crack, 0i is the nominal crack width of the ith crack, M is the total number of cracks, A(x) is the remaining area along the column and 0 < x < L is the spatial position along the column starting from the bottom of the column. The total length of the column is assumed to be L and the crack ratio ci is defined as the crack depth over the total depth of the intact column.
For a rectangular column, the residue flexural stiffness takes the form of where E 0 I 0 denotes the flexural stiffness of the intact column and E(x)I(x) is the residue flexural stiffness along the column. A similar form of residue stiffness can be obtained for other shapes of cross-sections.
In the present model, the severity of the crack is represented by the crack ratio parameter ci , the position of the crack is denoted by the position parameter x 0i and the crack width is described by the nominal crack width parameter 0i . This model can simulate cracks with varying severity, different positions and widths. As an example, the remaining total height and the residue flexural stiffness along long the column as a result of cracking is schematically graphed in Fig. 1. On both sides of the crack, the total height and the stiffness are the same as the intact beam. Besides, this mathematical model is arbitrarily differentiable and convenient to deal with analytically.
Under the impact of loads, such as wind, earthquake, vehicle vibration and rain and the effects of corrosion, temperature and material decay, old column-like structures (for example, high piers of high bridges in mountain areas) tend to develop open cracks on both sides of the columns, which are approximately symmetrical on both sides of the element. An exaggerated schematical crack model in a segment of a column-like structure is shown in Fig. 2. It is critical to evaluate the residual buckling load to ensure the reliability and serviceability of the structures. Besides, in this work, the classical assumptions for the Euler-Bernoulli beams are assumed.
It is worth noting that the crack width parameter oi is the nominal crack width. It is assumed that the crack shape approximately follows the famous Gaussian "bell-curve". In practice, assuming that the damaged area covered by the surface crack width is 99.9% of the total "bell-curve", the nominal crack width is defined as the width covering 50% of the "bell-curve". The surface crack width is measurable. Therefore, the nominal crack width can be obtained according to the following equation.

Differential equation governing the buckling load
The study of buckling of a cracked column is explained in this section considering the general case of an elastic column under an axial force. For this purpose, the moment and shear force are assumed at both of the ends of the column and the free body diagram of the lower segment of the column is schematically graphed in Fig. 3.
Assuming there is no transverse force along the column, then the shear force is invariant with respect to the spatial position x, where V (x) is the shear force along the column.
Taking moment equilibrium with respect to the bottom end of the column gives In this paper, the relationship between the moment and curvature of a cross-section is assumed to follow the Euler-Bernoulli beam theory, Substituting Eq. 6 into Eq. 5 and differentiating Eq. 5 twice with respect to x , the fourth-order ordinary differential equation governing the buckling of a column of rectangular cross-section with finite number of cracks is written as [22] Substituting Eq. 2 into Eq. 7 yields the following differential equation governing the buckling of the column, In this paper, the non-dimensional form of the foregoing equation is studied. Introducing the non-dimensional coordinate = x L , the following non-dimensional parameters are defined And the dimensionless buckling mode is defined as Substituting the newly defined non-dimensional relationships into Eq. 8, the reformulated equation is expressed as denotes the non-dimensional eigenvalues and the primes indicate taking the derivatives with respect to .
After the buckling eigenvalues are obtained, the buckling load can be calculated as To simplify the expression Eq. (12) is rewritten as where In the diffused crack model, the continuity of deflection, slope, moment and shear force are automatically satisfied. In this paper, it is assumed that the cracks keep open during the loading process.

He's variational iteration method
In this section, the powerful Variational Iteration Method (VIM) is reached to obtain the approximate analytical solution of the equation governing the buckling of the column. VIM is a powerful method for the solution of linear and nonlinear ordinary and partial differential equations. It is a good candidate for obtaining the approximate analytical solution for problems where the exact solution is difficult to obtain. ∼ n = 0(n = 0,1, 2, … ) . In the generalized Lagrange's multiplier ( ; , ) , is the variable while and are parameters. In this work, a semicolon after and a coma between and is used to denote the difference.
To solve Eq. 14, the following iterative formula with the correction functional according to the variational iteration method is constructed Taking variation with respect to n on both sides of Eq. (17) gives The extremum condition of n+1 ( ) results in δ n+1 ( ) = 0 , based on which the stationary equations are obtained by comparing the coefficients as and the boundary conditions are where the primes indicate taking derivatives with respect to the variable . The optimal generalized Lagrange multiplier obtained in this paper is more efficient than the cubic form ( − ) 3 6 derived by Coşkun and Atay [23] when investigating the critical buckling loads of elastic columns of constant and variable cross-sections. In this work, the term 2 �� n ( ) is also taken into consideration when finding the generalized Lagrange multiplier. In this paper, the optimal generalized Lagrange multiplier is also dependent on the eigenvalue parameter and the position parameter .
Substituting ( ; , ) into the iterative equation in Eq. 17 and removing the restrictions leads to the following iteration expression with the explicit generalized Lagrange multiplier for the buckling mode shape Assuming the initial approximate analytical solution is 0 ( ) , substituting it into Eq. (25) and applying Leibniz's rule for differentiation under the integral sign and integration by parts, the first-order solution is obtained as The second-order solution is expressed as

The buckling mode of an intact column
In this section, the buckling mode shape of an intact uniform column under different boundary conditions is obtained by the Laplace transform, which will serve as the initial approximate analytical solution. The fourth-order differential equation governing the buckling of a uniform column is Performing the Laplace transform of the above equation gives where Φ(s) = L{ ( )}s is the Laplace transform of the transverse displacement and s is generally a complex parameter.
Performing the inverse Laplace transform, the general solution of the buckling mode shape function can be determined as The buckling mode shape of the intact beam is determined by the four constants according to the boundary conditions concerned and the two boundary conditions at the bottom are automatically considered when applying the Laplace transform. When applying the other two boundary conditions on the top, the buckling shape can be determined. The buckling mode shape function in the above equation serves as the trial function for the Rayleigh-Ritz method and the initial approximate analytical solution of the variational iteration method.
Then, by imposing the two boundary conditions at the top end of the column, two homogeneous linear equations are obtained for different combinations of boundary conditions. In this work, four cases of boundary conditions are considered, namely, Simple-Simple (S-S), Clamped-Clamped (C-C), Clamped-Simple (C-S) and Clamped-Free (C-F). As the buckling mode shape functions and buckling load equations will be utilized in the variational iteration method and the Rayleigh-Ritz method, they are presented below.
For the Pinned-Pinned column, the buckling mode shape function is

The characteristic buckling load equation is
For the Clamped-Clamped column, the buckling mode shape function is (29)

The characteristic buckling load equation is
For the Clamped-Free column, the buckling mode shape function is

The characteristic buckling load equation is
For the Clamped-Pinned column, the buckling mode shape function is The characteristic buckling load equation is

Buckling of cracked columns
Assuming the buckling mode shape function of a cracked column takes the following form where Ψ j ( )(j = 0, 1, 2and3) are four shape functions that are to be determined.
Substituting Eq. 40 into Eq. 27, one obtains the iterative equation for the functions Ψ j ( )(j = 0, 1, 2and3) as where the first letter of the sub-indices in Ψ j,n+1 indicates the number of the coefficient and the second one denotes the order of iteration.
In this paper, the author presents the buckling mode shapes functions and the characteristic buckling load equations of four frequently encountered boundary conditions. The corresponding buckling mode shapes and buckling load equations are given below, selecting the mode shape function of the intact uniform column (34) as the initial approximate analytical solution of the mode shape of the cracked beam. In each case, the two equations satisfying the boundary conditions of the top end of the column ( = 1) can be written into a matrix form as follows. For a nontrivial solution, the determinant of the coefficient matrix should be zero and the characteristic buckling equation is generated, which is an eigenvalue problem. The buckling mode shape functions corresponding to each eigenvalue can then be written out.

Case I: (P-P)
The buckling mode shape function is written as The characteristic buckling load equation is obtained as

Case II: (C-C)
The buckling mode shape function is written as The characteristic buckling load equation is obtained as

Case III: (C-F)
Considering the boundary conditions, the following homogeneous equations are obtained, The buckling mode shape function is written as The characteristic buckling load equation is obtained as

Case IV: (C-P)
The buckling mode shape function is written as The characteristic buckling load equation is obtained as As is observed from the above equations, one advantage of the diffused crack model and the variational iteration method is that the highest rank of the determinant is two for all four cases of boundary conditions, which is convenient to apply.

Rayleigh-Ritz method
To check the solution by the variational iteration method, the Rayleigh-Ritz technique is implemented to obtain the critical buckling load and transverse bending deflection of the buckling column. The trial buckling mode functions corresponding to varied boundary conditions obtained by the Laplace transform are implemented.
The strain energy stored in the column can be written as The work done by the external load is written as According to Hamilton's principle, the total energy of the buckling column is obtained as Assuming the deflection of the columns takes the form of where N R is the total number of terms that are implemented for the Rayleigh-Ritz method, a k is the unknown participating coefficient for each term and k ( ) is arbitrary admissible buckling shape functions satisfying the boundary conditions. In this paper, sixteen terms are applied for numerical calculations.
Substituting the trial functions obtained by the Laplace transform and minimizing the total energy function, the buckling load can be determined. Taking partial derivatives with respect to the coefficients, a system of linear equations is obtained and the characteristic buckling load can be calculated.

Results and discussion
In this investigation, the second-order solution is adopted to compute the buckling mode shape functions and the buckling loads. Four cases of combinations of boundary conditions (P-P, C-C, C-F and C-P) are considered in the previous section. In this section, the dimensionless form of the equation governing the buckling of a column is implemented and the total length of the column is normalized.

Single crack
In the case of a single crack, the author keeps the nominal width of the crack to be 0.01 and set the crack ratio equal to 0.2. The buckling load corresponding to the varying position of the crack along the column is computed. The calculations of the first two buckling loads suing VIM and the Rayleigh-Ritz method are tabulated in Table 1. It is observed that the symmetry for the Pinned-Pinned and Clamped-Clamped boundary conditions is well preserved. The percentage errors between the results of the variational iteration method and the Rayleigh-Ritz method are below 0.2% for all cases of boundary conditions considered.

Two cracks
In the case of two cracks, the first crack is located at 1/3 and the second at 2/3 of the beam length, respectively. The nominal crack widths for both cracks are set to be 0.01 and the crack ratios are 0.2. The buckling loads of the first five modes calculated according to the variational iteration method and the Rayleigh-Ritz method for four different boundary conditions are compared in Table 2. It can be seen that the differences between the two methods for all five modes are below 0.35% for all cases of boundary conditions. The calculated results in the table show good agreement between the VIM and the Rayleigh-Ritz method.

Influences of nominal crack width
To show the influences of the nominal crack width, a crack ratio of 0.2 is selected. The crack is located at 1/3 of the column from the bottom. The buckling loads for the first two modes calculated by the variational iteration method and the Rayleigh-Ritz method are compared in Table 3 when the nominal crack width is varying from 0.005 to 0.025 with an interval of 0.005. From Table 3 it can be seen that the buckling loads obtained by the second-order VIM solution agree well with those obtained by the Rayleigh-Ritz method and the maximum error between the two methods is 0.4%.

Buckling mode shapes
To illustrate the buckling mode shapes of a cracked column, a crack ratio of 0.2 and a nominal crack width of 0.01 are selected and the crack is located at 1/3 of the column from the bottom. The buckling mode shapes for eigenvalues within 20 are graphed in Figs. 4, 5 , 6 and 7, together with the buckling mode shapes of the intact column. For the Pinned-Pinned boundary conditions, the slope � (0) is set to 1 when computing the buckling mode shapes, while for the clamped-clamped, clamped-free and clampedpinned columns, ��� (0) is made to be 1.

Conclusion
In this study, the diffused crack model is adopted to formulate the differential equation governing the buckling of a column with multiple cracks. The variational iteration method is applied to solve the obtained differential equation based on the solution of the intact column, which is obtained by Laplace transform. The second-order analytical solution is derived and the calculated buckling loads are compared with those from the Rayleigh-Ritz method.
To justify the accuracy of the proposed method, the computed buckling loads for the proposed variational method

Compliance with ethical standards
Conflict of interest The author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.