An efficient analytical model to evaluate the first two local buckling modes of finite cracked plate under tension 1

The analytical approach is presented for both symmetric and anti-symmetric local buckling of the thin-plate in finite sizes and with a center crack under tension. An efficient classical solution based on the principle of minimum total potential energy was provided using only 2 and 1 degrees of freedom for symmetric and anti-symmetric modes and the linear elastic buckling loads are evaluated by the means of RayleighRitz method. In the pre-buckling state, a correction factor for the peak compressive stress in the finite cracked plates is defined with an empirical formula and used in the analytical solution of the buckling. To verify the analytical approach, a wide range of numerical results by aid of finite element method are provided herein and a comparison between theoretical results with the experimental work of other researchers has been done. Both numerical and experimental results accept the accuracy and validity of the presented analytical model.


INTRODUCTION
In a cracked plate subjected to the tension load, perpendicular to the crack orientation, the developed compressive stress near the crack leads to buckling of crack edges.This phenomenon which can be called as local buckling of cracked plates under tension, has been studied by numerous researchers.
The early works were the experimental studies such as the Air Force Flight Dynamics Laboratory (AFFDL) results in AFFDL-TR-65-146 (1965); Zielsdorff and Carlson (1972).Similarly, Dyshel (1990;1999;2002) has presented some variety of experimental results on the stability and fracture of plate with a central crack referred to as a non-classical problem.Structures 12 (2015) 2078-2093 Besides, the numerical techniques such as the finite element method (FEM) have been carried out by many investigators such as Markstrom and Storakers (1980); Sih and Lee (1986); Shaw and Huang (1990); Brighenti (2005a;2009) to evaluate the critical load in various geometry and boundary conditions.Furthermore, Riks et al. (1992) studied the buckling and the post-buckling behavior from the viewpoints of stability and fracture mechanics by the aid of finite element method.

Latin American Journal of Solids and
As Cherepanov (1963) mentioned, an analytical solution for this problem is engaged to some mathematical and fundamental difficulties because the exact solution for in-plane fields of stress and displacement exists just for the infinite plates and mathematically it is complex as well as the local buckling shapes.Guz et al. (2004) presented an analytical solution for the buckling of infinite cracked plate under tension.For the finite cracked plate, Brighenti (2005b) used an approximate model with estimated pre-buckling stress functions.In that work, the stress functions are in the form of hyperbolic series with infinite terms.The buckling shape function is considered to be a radial Gauss-like function (same function of x and y ) in which the values of the two unknown parameters were assumed by the author.Both Guz et al. (2004) and Brighenti (2005b) consider the first mode of buckling, which is called the symmetric mode.
Following the previous studies, the most of non-experimental works concern the numerical studies which use many degrees of freedom to solve the problem.Against, the current paper intends to present an efficient classical solution for the non-classical problem of the central cracked plate buckling under tension, which includes only 2 and 1 degrees of freedom for both first (symmetric) and second (antisymmetric) modes of buckling and has answer for the finite plates with various geometric ratios.In this direction, primarily, the pre-buckling stress state is approximated by a simple function and then the proper shape functions with minimum degrees of freedom are used to describe the buckling mode shapes.Finally, the total potential energy is written as a function of the pre-buckling stress and the buckling deflection and subsequently, the well-known Rayleigh-Ritz method is used to evaluate the critical loads for each mode of buckling.All results were validated with several numerical finite element simulations in ABAQUS environment and some existing experimental works.

General considerations
Figure 1 shows the finite cracked plate, with the thickness t , under tension load.The geometric parameters are shown in this figure as well.As shown, the central crack orientation is perpendicular to the load direction which corresponds to the lowest buckling load and in terms of geometry it has two axes of symmetry consistent with the X and Y .
It is assumed that the material is isotropic and linearly elastic with Young's modulus E and Poisson's ratio ν .Since ij σ denotes the pre-buckling stresses, from Von Karman's linearized theory, the total potential energy Π is obtained based on the transverse displacement ( ) , w x y as the following equation: in which, ( ) is the bending stiffness and subscript of w denotes the derivation.As can be seen, the maximum deflection occurs at the center of the crack for both modes and generally, the deformations are insignificant in the regions away from the crack, because the instability of the crack edges is a local buckling phenomenon due to the compressive stress near the crack.Thus, providing a mathematical solution to this problem requires reviewing and determining the status of pre-buckling stress and especially compressive stress along the crack edges.

Pre-buckling stress state
The exact solution for the stress field of an infinite plate with a crack under the in-plane loads exists in many references such as Sadd (2005) but for the finite plates, due to the complexity of the problem, it cannot achieve such relations.In Figure 3, the distribution of stresses x σ and y σ on the axes X and Y of a finite cracked plate schematically are compared with the infinite one under the same tension stress o σ .It is necessary to note that due to the symmetry, the tangential stress on the symmetry axes is zero and just normal stress exists.In this model, the normal stress x σ acting on the Y -axis includes the local compressive stress parallel to the crack edges that is the main reason for the local buckling and should be considered in the analysis.
For some finite cracked plates problems, infinite plate relations can be used with some simplifications such as applying a correction factor related to the finite size of the plate.For example, in the field of fracture mechanics, the correction factor β is used for the singular stress field at the crack tip of the finite plates.The values of β can be assessed using numerical methods and for different geometric ratios a W , L W are presented in Isida (1971).In this study, another correction factor β ′ is defined and used for the peak compressive stress.From the exact solution of infinite cracked plate under tension in Sadd (2005), the stress distribution x σ on the axis Y is obtained as follows: ) According to the equation ( 2), the maximum compressive stress is obtained for 0 Y = which the value is equal to the applied tensile stress o σ .So by using the correction factor β ′ , the peak compressive stress for the finite cracked plate is achieved as o β σ ′ × .In section 3.1, the values of β ′ are obtained for different ratios a W , L W using the finite element method.

Simplified analytical model
A simple analytical model that is used to estimate the buckling load of the finite cracked plate is shown in Figure 4.It is clear that one quarter of the plate is considered and the distribution of the stress x σ on the Y -axis is represented as a simple bilinear function.For both finite and infinite cracked plates, the stress x σ along Y -axis changes from a negative value to a positive one and at a special point (e.g.Y d = ) it is zero (see Figure 3).From equation ( 2), the zero-value location d for the infinite plate can be found by setting 0 x σ = which results in It is assumed that the value of d depends only on the crack length so that the constant 0.8 d a δ = ≈ can be defined and used for the finite plate with arbitrary dimensions as shown in Figure 4.This assumption is further verified in section 4.1.Now, using dimensionless parameters y Y a = and , the simple bi-linear function x σ with the value of σ at y = ℓ and zero for y δ = can be written as In the bi-linear function, c σ is the peak compressive stress and can be achieved as c o σ η σ = × , where η equals β ′ multiplied by a reduction coefficient C related to the linearization of x σ .In this study, roughly C = 0.9 is used for all analytical cases but the exact value can be set by a regression analysis on the results.Apart from that, using the force equilibrium along the X direction, t σ value is obtained in terms of the variable c σ according to the following equation.
As shown in Figure 4, the analytical one quarter model of the cracked plate is divided into cracked ( 0 X a ≤ ≤ ) and un-cracked (a X W ≤ ≤ ) rectangular panels which are named panels 1 and 2, respectively.Since in the local buckling mode, the buckling deformation must be vanished with the distance from the crack, the buckling deformation in panel 2 is limited and panel 1 will experience the most bulging.To define the buckling shape functions ( ) , w x y for each mode of buckling, the conditions presented in Table 1 are considered.For the symmetric mode, the following shape functions are defined in terms of the dimensionless variables x X a = , y Y a = and w W a = .
Latin American Journal of Solids and Structures 12 (2015) 2078-2093   ( is the maximum deflection and where A is an arbitrary parameter and a continuing condition at . In this case, the division line 1 x = acts as a simple support.It is necessary to note that these functions are obtained through multiplying an exponential decay that exerts locally the nature of buckling, to polynomial functions so that the overall satisfies all conditions in Table 1.
Table 1: Boundary conditions of buckling shape modes at panels' edges.
From Table 1, it is clear that for the anti-symmetric mode, the panel 2 has four edges with zero deflection and it is not unreasonable to assume ( ) , 0 w x y = for 1 x ω ≤ ≤ .In other words, the division line x = 1 acts as a clamped support.This assumption simplifies the solution for anti-symmetric mode so that the shape function can be lessen as To calculate the potential of applied stresses x σ and o σ shown in Figure 4, the first assumption is that the stress x σ acts throughout the horizontal shortening (in X -direction) due to the buckling deflection of panel 1.For the stress o σ which acts in the vertical direction (Y -direction), since the crack edge is free to move, the vertical shortening due to the buckling deflection appears as crack mouse opening.Consequently it can be assumed that the stress o σ is relatively without moving so that the related potential can be neglected in the analysis.Now, using equation ( 1) and the well-known Rayleigh-Ritz method, the critical loads are obtained for both symmetric and anti-symmetric modes.It should be noted that, for the symmetric and antisymmetric shape modes, just 2 and 1 degrees of freedom are used respectively which can greatly reduce the computation.

EVALUATIONS USING FEM
A wide range of models were prepared in ABAQUS environment in order to determine the correction factor β ′ and numerical evaluation of the buckling loads in which the effect of parameter L W , which so far has received less attention, is also included.
After some mesh size sensitivity analysis with relation to the buckling load, the proper finite element mesh with adequate refinement near the crack tips is illustrated in Figure 5 in which, 4-node doubly curved general-purpose shell elements S4R have been used.

Peak compressive stress correction factor
In section 2.2, the correction factor β ′ was defined for the peak compressive stress in the finite plates under tension as that amount can be calculated using the finite element models.Given that the prebuckling stress field is independent of the material properties, the β ′ is only related to the geometry of the plate.Therefore, several models have been developed with different aspect ratios a W = 0.2 0.5 ∼ and 1 2 L W = ∼ , the results are presented in Table 2. Like the infinite cracked plates, the value of β ′ for small cracks in comparison with plate dimensions should approach to 1 which is inserted in the second row of the above Table 2.
Latin American Journal of Solids and Structures 12 (2015) 2078-2093 Due to the use of β ′ in the proposed solution, it will be useful to determine an approximate formula for β ′ based on the numerical results.Using the method of separation of variables, the general form of β ′ can be considered as follows.
( ) ( ) Each of the functions F and G has certain limit values which are useful to determine their overall form.For example, by changing a W from 0 to 1, ( ) F a W must change from 0 to infinity respectively and similarly ( ) G L W changes from infinity to 1 when L W changes from 0 to infinity.Due to the described conditions, several functions were considered for F and G and they were analyzed using numerical results and applying least squares method.Finally, appropriate forms of the functions F and G are recommended below.

( ) )
( ) Due to the simple form and precision of the defined formula for β ′ , it can very well be used in the analytical solution presented in section 2.

Buckling loads for the symmetric and anti-symmetric modes
Using equation ( 1) with the Rayleigh-Ritz method, it can be easily shown that the buckling load is directly related to the elastic modulus E and inversely with 2 t (the thickness to the power of 2).Thus, the other parameters involved in the problem are the crack length a, the dimensions of the plate L , W and the Poisson's ratio ν that their influence must be evaluated on the buckling load.For this purpose, finite element models have been developed with different Poisson's ratios ν = 0.1 0.4 ∼ and the aspect ratios 0.2 0.5 a W = ∼ , 1 2 L W = ∼ and the values of buckling load are presented in Table 3 as the dimensionless quantities.It should be noted that in these models, the constant parameters are E = 200 GPa, t = 1 mm, W = 200 mm and all exterior edges are considered to be simply supported.These results are used to validate the analytical results and are discussed in section 4.

Pre-buckling stress distribution and the correction factor β ′
In Figure 6, the curves of correction factor β ′ from equation ( 10) are drawn for the different ratios of a W and L W and compared with the finite element results from  10) and FEM from Table 2.
Figure 6 shows that, the approximate equation ( 10) for β ′ is in good agreement with the finite element results so that the maximum difference is less than 1.5 % in all cases.Also, by changing any of the parameters a W and L W , the FEM results will follow the expected and defined trend for β ′ function.For example, with decreasing a W , the value of β ′ decreases and approaches to 1 for all ratios of L W . Overall, the curves in Figure 6 show the effect of plate dimensions and the crack length on the peak compressive stress, so that, the peak compressive stress increases with increasing the crack length and decreasing the length of the plate in comparison with plate width.
Apart from the correction factor β ′ , Figure 7 shows the normalized actual distribution of stress  FEM results show that for different crack length the distribution of the stress x σ has been changed, but the relative zero point of stress on the Y -axis is almost constant and about 0.8 Y a = . Such value was obtained in section 2.3 for infinite plates and was used as 0.8 d a δ = ≈ to define the bilinear function.Figure 7 also show that, the curve trend for the infinite cracked plates is descending far from the crack, so that, the positive stress x σ tends to zero next to the plate edge.Against, FEM results indicate that this trend is ascending for the finite cracked plates especially for relatively larger cracks compared with the length of the plate.This represents the effect of finite size of plates on the pre-buckling stress state.Comparing the curves obtained for different ratios of a L implies that proportional to increase in the peak compressive stress, the maximum tensile stress on the Y -axis increases which is consistent with the equilibrium condition in the X direction written in section 2.3.The larger view of the compressive region of stress x σ shows that under section 2.3, the peak compressive stress in the equivalent bilinear function is multiplied by a factor of 0.9.

Symmetric and anti-symmetric buckling loads
In this section, in order to demonstrate the accuracy and ability of the presented analytical model to estimate the buckling load of the cracked plates with different geometric and material properties, fairly complete comparison between analytical model (Theory) and FEM results are depicted graphically in Figures 8 to 10.
Curves depicted in Figure 8 show that with increasing the crack length, buckling capacity of the plates will drop drastically.The reason is that with increasing the crack length, on the one hand, the coefficient β ′ increases and on the other hand, the geometric stiffness of the plate decreases and both factors are simultaneously reducing the buckling capacity of the plate.It is necessary to note that this drop in the anti-symmetric mode is greater than the symmetric mode, so that with increasing crack length, the difference of the two modes decreases but never reaches zero.Therefore, the buckling load in the symmetric mode is often less than the anti-symmetric mode and the symmetric mode can be introduced as the first and prior mode.In contrast, with decreasing the crack length, the buckling load is rising with a steep gradient and it is expected that it tends to infinity at 0 a W = because the plates without crack never buckle under tensile load.
Overall, the results of the presented analytical model for different values of crack length have an acceptable agreement with the FEM results in terms of quality and quantity which indicate that the simplifying assumptions of the analytical model are reasonable and do not affect the accuracy.
The curves plotted in Figure 9 show that increasing the ratio L W or, in other words, increasing the length of the plate will always increase the buckling capacity.But, the sensitivity of the buckling load to the plate length is reduced by increasing the ratio L W , so that, it looks for 2 L W > , the buckling load is almost constant.
By comparing Figure 9 with Figure 8, it can be concluded that the buckling load is less sensible to L W than to a W , however, the analytical model is able to show that slight variations of buckling load as well and is in good agreement with the FEM results.
According to the Figure 10, the results of both the presented analytical model and FEM show that the sensitivity of the buckling load to Poisson's ratio is less than its sensitivity to a W and L W .Moreover, in the Figure 11, the presented approach is compared with some existing experimental works on the symmetric buckling mode, which will be described.The first experimental data was reported by Air Force Flight Dynamics Laboratory (AFFDL) in AFFDL-TR-65-146 (1965) as load versus buckling deflection curves and in the present paper, the extended Southwell technique from NACA TN-658 (1938) is used to estimate the buckling loads from those curves.In the work of Dixon and Strannigan (1969), the upper and lower limits of buckling load are given and here the upper limit is considered.The other experiment was done by Zielsdorff and Carlson (1972) in which the buckling loads were presented and used directly in this paper.The characteristics of these samples are listed in Table 4 and ν = 0.33 is considered for all samples.
From Figure 11, it can be seen that the results of AFFDL-TR- 65-146 (1965), with a small distance, are at the top and bottom of the present theoretical curve, so that the theoretical results are well consistent with the mean of the experimental results.Dixon and Strannigan (1969) experiments include a wide range of relative crack lengths and are in very good agreement with the theoretical curve.In the work of Zielsdorff and Carlson (1972), a fewer range of cracks were considered and the difference between the experimental results and the theory increases with decreasing the relative crack length.
Looking at three experimental works we can see that for the small cracks, the experimental results show lower values than the analytical results or in other words, the analytical results are somewhat overestimated for small cracks.In sum, given that the experiments were carried out under different conditions and on different materials, Figure 11, can confirm the accuracy of the analytical model based on the experimental results.

CONCLUSIONS
In this paper, an approximate classical solution based on the principle of minimum total potential energy was provided for non-classical local buckling problem of central cracked plate under tension which includes both symmetric and anti-symmetric modes.Given that the mathematical functions for symmetric and anti-symmetric shape modes are written using only 2 and 1 degrees of freedom, the proposed approach mathematically is efficient and has an advantage in comparison with the numerical methods with many degrees of freedom.
After investigating the pre-buckling state, the relative zero point of normal stress x σ acting on the Y -axis (perpendicular to the center of the crack as Figure 3) was defined.Then using the finite element models, it was shown that this point is approximately constant and independent of the plate dimensions so that its theoretical value for the infinite cracked plate can be used for the cracked plates with any sizes.Moreover, the factor β ′ was defined as correction factor for the peak compressive stress in the finite cracked plates; then, in order to mathematically describe β ′ , an empirical formula was presented based on the finite element results.This correction factor was used in the analytical solution of buckling.
To verify the results of the presented analytical model, a wide range of numerical finite element models (FEM) were produced in ABAQUS environment which can evaluate the effects of the geometric characteristics and material properties of the plates on the buckling load.By comparing the results of these models with the results of the analytical model, the accuracy and the ability of the presented analytical model was demonstrated.
Finally, the results of three existing experimental works compared with the presented analytical model that confirm the validity of the analytical model to estimate the buckling load.

Figure 1 :
Figure 1: The finite plate with a central crack under tension.

Figure 3 :
Figure 3: Comparison of stress distribution in the one quarter of a finite cracked plate (solid line) and the infinite one (dash line) under same tension.(+) for tensile and (-) for compressive stresses.

Figure 4 :
Figure 4: Simplified model of the finite cracked plate.

Figure 5 :
Figure 5: Typical finite element mesh and crack tip detail.

xσ
along the normalized Y -axis obtained from FEM for square plates ( 1 L W = ) with different crack lengths ( 0.2 0.5 a L = ∼ ) in comparison with the equivalent bilinear function from equation (6) which is used in the buckling analysis.

Figure 7 :
Figure 7: Normalized distribution of stress x σ from FEM and equivalent bilinear function from equation (6) with a larger view of the compressive region.

Figure 8 :
Figure 8: Normalized buckling stresses versus relative crack length for both symmetric and anti-symmetric buckling modes.

Figure 9 :
Figure 9: Normalized buckling stresses versus relative plate's length for both symmetric and anti-symmetric buckling modes.

Table 2 :
β ′ for different aspect ratios from FEM.

Table 3 :
FEM dimensionless results of buckling load for both symmetric and anti-symmetric modes.
Latin American Journal of Solids and Structures 12 (2015) 2078-2093 Figure 6: Correction factor β ′ as a function of a W and L W from equation (

Table 4 :
Characteristics of experimental samples.