Analysis of Ultimate Compressive Strength of Cracked Plates with the Use of DoE Techniques

Abstract The objective of this work is to investigate the structural compressive response of plates with locked cracks accounting for all relevant factors and correlation between them. The nonlinear FE model considering both geometric and material nonlinearities is employed herein, and the FE model of the structural response of intact plates is validated with the available experimental data. In the common studies, based on One Factor at a Time analysis, some of the parameters and interactions between them are excluded. In the present study, the numerical investigations are conducted with the use of the Design of Experiments techniques, where all essential parameters and their interactions are adequately considered. With a total of 32 numerical analyses, the most influential factors and their interactions are identified. As a study outcome, empirical formulations, which allow for a fast estimation of the ultimate compressive strength of intact plates, plates with locked cracks, and repaired cracked plates, are derived. The developed formulations represent a fast and practical tool for estimating the ultimate compressive strength of intact, cracked, and repaired plates, which can be easily employed in the reliability analysis.results followed by their discussion.


INTRODUCTION
Ships and offshore structures are subjected to different degradation effects, such as corrosion and cracks and deformations [1,2]. The latter could be the result of fatigue damage or impact loads. The fatigue cracks are often undetected, and they can have an impact on the strength reduction of structural components.
The effect of locked cracks and their influence on the loadcarrying capacity has been investigated experimentally and numerically for some years now. Some of the first attempts to assess the importance of a transverse crack to the buckling capacity of plates were made in [3][4][5][6], showing the importance of this problem. In all cases, the buckling strength was significantly lower compared to non-cracked plates. Further studies were performed analysing the elastoplastic collapse, considering the geometry, and material nonlinearities. One of the more advanced studies was conducted by Paik et al. [7], where both numerical and experimental analyses were carried out. They found that the presence of cracks can significantly reduce the ultimate strength of plates (by up to 50%).
Additionally, the conservative estimation of strength reduction is related to the cross-sectional area reduction with the presence of a crack. However, due to the cycling load, which may change from tensile to compressive, in some cases the cracks may close, and then their behaviour is somehow similar to a non-cracked plate. Additionally, nonlinear FE analysis was found to be a useful tool in predicting the behaviour of cracked plates subjected to compressive load.
A quite comprehensive review related to the ultimate strength of cracked ship structural elements was presented in [8]. A limited number of experimental studies were conducted with regard to cracked plates [7,9] and cracked stiffened plates [10]. However, only a limited number of cases were investigated. The joint effect of locked cracks and openings was experimentally investigated in [11], showing that the presence of these two phenomena could result in a significant loss of plate capacity. Additionally, it is hard to model the boundary conditions in the experimental domain, and these may differ from the joining conditions of real ship structural components. Recently, more attention has been paid to the numerical analysis. The ultimate strength of different structural elements subjected to cracks was investigated, including plates [12][13][14][15][16], stiffened plates [15,17], stiffened panels [18], box girders [19] and ageing ship hulls [20]. Based on these studies, the governing factors that have the most influence on the structural capacity reduction can be identified. The governing parameters that may be identified with regard to the impact of the crack are the crack length, crack location, crack orientation, and the shape of the crack tip. In all studies, for transverse cracks, the reduction is higher with increase of the crack length. In the case of the longitudinal orientation of locked cracks, no significant effect is observed in most of the studied cases. However, some studies recognise that this parameter is a relevant one. In the case of the crack tip shape, in [21] and [22] it was shown that it will not influence the ultimate strength in general. However, the stress concentration around the crack tip will be influenced by the FE mesh density and tip shape.
Nevertheless, the factors that typically influence the ultimate strength of plates will potentially interact with the crack parameters, such as the plate slenderness ratio, plate aspect ratio, and welding-induced imperfections and residual stresses. The combined effects of cracking and initial imperfections were studied in [16]. It was concluded that for small cracks, the ultimate strength is additionally reduced for the higher levels of initial distortions. However, in the case of long cracks, the crack damage effect plays a dominant role.
Based on the different studies, empirical formulations were developed too. Babazadeh and Khedmati [21] derived empirical formulations for the ultimate strength reduction of transversely cracked plates as a function of the crack length, plate slenderness ratio, and plate aspect ratio. In contrast, Paik et al. [7] derived the formula as a function of the plate slenderness ratio and crack length.
Although in previous studies the importance of different individual factors was investigated, investigations that take into account all relevant factors and their joint effect seem to be lacking. Additionally, the studies usually used the OFAT (One Factor At a Time) [23] technique during the sensitivity analysis, i.e. they changed only one parameter, leaving the rest of them unchanged. This technique is very good for initial studies to examine the relevant parameters. However, the possible interaction between factors is lacking. To investigate that, the DoE (design of experiments) [24] methodology seems to be suitable. Another benefit of the DoE is the significantly reduced number of observations compared to OFAT analysis. Based on the DoE analysis, the response surface takes into consideration the relevant factors, and the interactions between them can be established with a significantly lower number of observations compared to classical regression analysis.
Depending on the considered experimental plan, the response surface may be linear or nonlinear.
The objective of the present study is to assess the ultimate strength of cracked plates subjected to compressive loads with the use of the nonlinear FE method and DoE techniques. At the beginning, all relevant parameters related to the locked crack, as well as the plate, are taken into consideration. However, based on the OFAT analysis some of them are not considered further due to their lower impact. The present study also investigates the importance of the variables related to the drilling holes, which are commonly used to stop crack propagation in emerging cases. Furthermore, once the most important governing parameters are identified, the ultimate strength assessment with the predefined DoE plan is conducted. As a result, the most influential factors and interactions are found. Additionally, the response surfaces, allowing for fast estimation of the ultimate strength of cracked plates as well as intact and repaired plates, are established, and the results are compared with other existing empirical formulations and experimental data from the literature.

FE MODELLING
The compressive structural behaviour of a cracked plate is modelled with the use of the nonlinear finite element method (FEM), considering both geometrical and material nonlinearities. The implicit static solver (Newton-Raphson iterative procedure) is used, employing the commercial software ANSYS [25]. The material is modelled as bilinear with hardening, and in the following the material properties of normal strength steel are considered, such as a yield strength (R e ) of 235 MPa, Young's modulus (E) of 206 GPa and hardening stiffness (E h ) of 430 MPa. This type of steel is commonly used in shipbuilding, and the material properties are considered as required by Classification Societies Rules, such as Common Structural Rules for Bulk Carriers and Oil Tankers [26]. The plate was modelled with the use of SHELL181 elements, and the contact between the edges was modelled with the use of CONTAC52 elements.
The initial imperfections are modelled considering one half-sine wave as suggested by Smith [27] with an average level of 0.1β 2 t, where β is the plate slenderness ratio, defined as: where w is the plate width, and t is the plate thickness.
At the beginning of the analysis, the optimum mesh density needs to be established. In this order, the element size convergence studies are performed considering the initial parameters of the analysis, as presented in Table 1. At this stage, the drilling holes are not considered, and the central transverse crack is taken into account. The gap between cracked edges is 2 mm, and the cracked tips are modelled as circular ones.
A typical plate being a part of the ship hull structure is spanned between longitudinal stiffeners and transverse girders. Based on that, the unloaded edges can be considered as simply supported. However, in the case of the transverse edges of the plate, since these are supported by very rigid girders, their behaviour can be assumed as something between simply supported and clamped conditions. In this study, the clamped boundary conditions are considered on the loaded edges.
To simulate the strain distribution appropriately around crack tips, the element size of 0.5 mm was found to be the right solution. The additional transition area around the crack is distinguished to provide the proper mesh distribution. The element size along the cracked edges is considered as 3 mm to avoid an excessive aspect ratio of the elements in that region. An example of an element mesh distribution around the crack tip is presented in Fig. 1.
The results of the mesh convergence studies are presented in Fig. 2. The normalised ultimate capacity of the plate is shown in the vertical axis, which is the ratio between the ultimate stress in the plate and the yield stress of the considered material. It can be seen that, with the mesh refinement, the ultimate strength tends to arrive at the actual value. Based on these results, the 20 mm element size is chosen for further analyses, as it provides both accurate results and quite low computation times.
The plate FE model of the considered element size of 20 mm and initial analysis parameters from Table 1 are presented in Fig. 3, where the mesh distribution around the crack edge is also visible. In Fig. 4, the von Mises stress distribution together with the deflection shape of the considered model are presented. As can be seen, the stress concentration around crack tips is significant. Additionally, in the crack position, the deflections are very high, leading to different shapes of plate and levels of deflections compared to the intact plate.
Welding-induced residual stresses may also reduce the plate strength, and are often taken into account [28]. However, based on different studies, it is evident that after multiple cycles of loading the residual stresses are significantly reduced, which is called a shakedown effect [29]. Thus, when the crack initiates and starts to propagate, the influence of residual stresses is already released. For that reason, the residual stresses are not considered in this study, which may lead to an overestimation of the plate strength.

VALIDATION OF INTACT PLATE MODEL
To validate the FE model, the results of intact plate calculations are compared with the experimental investigations of [30], where plates with different aspect ratios and slenderness ratios were investigated. The material properties considered in that experiment were based on tested coupons, and are equal to E = 206 GPa, Re = 290 MPa. For the boundary conditions, all edges were simply supported. The level of initial imperfections was not studied, so in the FE model, the mean initial imperfections are considered, as mentioned in the previous section. The plate dimensions and results for both numerical and experimental results are presented in Table 2. A total of 9 specimens were considered, with the same plate width of 400 mm. Additionally, the results were compared with the formula of Faulkner [31], where the classical formulation for the critical buckling stress of an infinitely long thin elastic plate was extended in terms of the ultimate strength and considering the mean level of the initial imperfections. The ultimate strength, in that case, is a function of the plate slenderness ratio and yield strength: As can be seen, the deviations between the numerical and experimental results are minimal, and this can be the result of the non-ideal behaviour of the boundary conditions during the experiment and unknown values of initial imperfections. Additionally, the mechanical properties considered in the FE model are the mean value, whereas, in real conditions, they are subject to uncertainties. However, it can be seen that Faulkner's simplified formulation underestimates the ultimate strength.
Based on that study, one can conclude that the FE model can predict the behaviour of a real plate accurately.

INITIAL SENSITIVITY ANALYSIS
Based on the literature review, it is evident that the plate aspect ratio, plate slenderness, crack length, initial imperfections level and crack orientation are undoubtedly the most important governing factors on the ultimate strength. However, in the case of the crack position, the conclusions are divergent. Additionally, the influence of the presence of drilling holes has not been studied previously. Considering these three parameters related to the longitudinal and transverse position of the crack and the drilling holes diameter, an OFAT analysis is performed to determine their importance. The rest of the parameters considered here are presented in Table 1, and the variable ranges are given in Table 3.
The plate with 20 mm drilling holes is presented in Fig. 5. The holes are modelled in such a way that they are not extending the length of the crack. The considered crack positions are presented in Fig. 6.
The results of the sensitivity analysis are presented in Fig. 7. As can be seen, the drilling holes have almost no influence on the ultimate strength. Nevertheless, a bigger drilling hole presents a lower stress concentration level. This will have a significant influence in the case of fatigue strength. In the case of the crack position, some influence is visible. In the case of the transverse position, the ultimate strength is higher for the middle crack location compared to the side crack position. In the case of the longitudinal position, the capacity is higher for a crack located near the loaded edge. It may be observed that the relations are not linear, and for the mean positions of the crack, the results do not deviate much from the results of the middle crack position. The differences between the maximum and minimum values are around 6%.  Apart from some significance of the crack position, this effect will not be analysed further. The influence seems to be much less compared to other governing variables, which can be concluded also from the previously mentioned studies.

DOE MODEL
Having established that five variables are the most influential ones, the proper engineering design experimental plan needs to be chosen. A full factorial analysis (2 k ) seems to be an excellent option to analyse the design space and provide information about the interaction between effects [32]. In DoE, the factors considered are referred to with letters from the alphabet. The factors with the minimum and maximum values are presented in Table 4.
The plate aspect ratio is considered between 1 to 4. In the case of the slenderness ratio, the range is based on the statistical data [33]. The typical values of that parameter in the case of ship structures are between 1.4 and 2.5. The corresponding thickness values for the slenderness ratios from Table 3 are 15 mm and 8 mm, respectively. In the case of the crack length, the maximum length corresponds to 40% of the plate width. The crack length cannot be too long because it will not be possible to stop its propagation [7]. Thus, 40% of the plate width is considered as a critical crack length. In the case of the crack orientation, 0 degrees corresponds to a transverse crack, while 90 degrees corresponds to a longitudinal crack. The level of the initial imperfections is estimated by w = w 0 β 2 t where the range of w 0 is based on the statistical analysis that may be found in [27].
The output of the analysis is the normalised ultimate stress. Each combination of factors considering either their minimum or maximum values is considered in the full factorial analysis. This leads to a 2 k number of cases, where k is the number of factors. In this case, five factors are taken into account, so 2 5 = 32 cases are computed. The analysed points are only those at the extreme values of the variable range, so linear behaviour is assumed in the region between the observations. The upper limit of the variable is considered as +1, whereas the lower limit of the variable has the value of -1. Having designed a test matrix, the FE simulations are carried out, and the values of output normalised ultimate stress are presented in Table 5.  Based on the results from Table 4, further statistical analysis is carried out. First, the effect of each factor and interaction is evaluated to analyse whether it is significant or not by calculating the average responses: (2) where Response + is the response where either the factor or the multiplication of factors is positive, and Responseis the response to the negative value of the factor or their multiplication. Accordingly, n + is the number of cases for positive responses and nis the number of cases for negative responses. In the presented study, n + = n -= 16.
In To evaluate the significant effects, several methodologies may be used. The first of them is based on the statistics that are commonly used in the DoE methodology. The cumulative probability value is calculated for the sorted effects and plotted at a half-normal plot [34], as presented in Fig. 8. The significant effects are those that diverge from the dashed line, which represents the probability of observation noise. That means that all effects that are lying on that line are stochastically originated. The effects that deviate from the dashed line are B, C, A, D, CD, AE, BE, DE, E, and AB.
Another possible way is to perform the t-test for each effect (to find whether the difference between the population means is significant or not). In this way, the t-values are used to measure the size of the difference between means concerning variation in the data. The t-value for each effect is calculated as follows: where MS Res is the mean of the sum of squares of the residuals (insignificant) effects, equal to: where DOF is the number of degrees of freedom, and is equal to the number of insignificant effects, and N is the total number of observations (equal to 32 in this case). The DOF in the present study is 21. MS Res is calculated, and it is 0.0316.
Based on that, the t-values for the analysed effects are shown in Table 5. To quantify whether the effect is significant or not, the critical t-value is calculated, which is equal to a two-tailed t-value for a 0.05-probability and the given value of the degrees of freedom. In this case, the critical t-value is 2.074. All effects with t-values higher than the critical one can be considered as important. To distinguish that, the Pareto chart is plotted, as shown in Fig. 9. The results of this analysis are in line with the half-normal plot shown in Fig. 8.
Based on the importance analysis, the significant and insignificant effects are distinguished, and the significant effects are shown in Table 6. The effects are ranked from the most important one to the least important.  As can be seen, ten effects from the overall number of 31 were found to be significant. The essential factor is the plate slenderness ratio, which is known as a parameter that has the most impact on the ultimate strength in general. The three next factors (crack length, crack orientation, and plate aspect ratio) have a similar impact on the ultimate stress. In general, five main effects and five interaction effects are found to be influential. Additionally, in the case of the interaction effects, only interactions consisting of two factors are considered.
From all the interaction factors, the most important one is the interaction between the crack length and crack orientation. This is obvious because, for a long crack, the structural capacity is severely reduced where it is a transversely oriented crack and slightly reduced when it is a longitudinally oriented crack. Similar observations can be found in [8,22]. What is interesting is that the impact of the initial imperfections level is smaller itself in comparison to its interaction with other factors. The two interaction effects concerning the initial imperfections are influential, with the plate aspect ratio and plate slenderness ratio due to their mutual magnification effect. The interaction between the plate slenderness ratio and crack orientation is also essential. In the case of slender plates, for transversely oriented cracks, the reduction of the structural capacity is not very severe in comparison to stockier plates.
Based on the analysis of the effects, the Response Surface [35] can be established. Since the preliminary plan was based on the Full Factorial Design and considering the main effects and interactions of two factors, the function will take a general form of: (6) where x i corresponds to the model variables and the parameters B i and B ij correspond to linear coefficients, and are equal to half of the considered effect values. It needs to be noted that, in this case, the coefficients may be positive or negative, depending on the effect sign. In the importance analysis, the absolute values of the effects were taken into account. B 0 is calculated as a mean value from all observations.
One can see that, in the DoE analysis, the variables x i were in the range between -1 and 1. However, as can be seen in Table 3, the range of the real variables is different. For that reason, the variables need to be normalised, and their values are estimated to: The response surface gives the estimation of the normalised ultimate strength. As mentioned previously, a linear relationship in the factors range is assumed here, which can be different, considering the FE analysis. To get a higher level of estimation, a Central Composite Design can be considered. In this type of design, additional central and axial points are considered, and the total number of observations is equal to 2 k + 2k + 1. Based on the results of this analysis, one can obtain the polynomial response surface, which is more advanced compared to the linear one. To avoid the more complex polynomial response surface, the formula from Eq. (12) is verified for the central and axial design points.
To expand the Full Factorial Design into a Central Composite Design, additionally 2k + 1 axial and central points need to be calculated. In the presented study, the number of additional points is 11. For these points, both exact FE results and estimations from Eq. (12) are presented in Table 7.
The difference between the exact FE solution and the estimation from Eq. (12) is presented in Fig. 10. One can see that for both the FFD points and the central and axial points, the differences are not very significant. The mean error of the estimation for the FFD points is about 3.6%, whereas for the central and axial points it is about 4.8%. The correlation factor of the linear regression is equal to 0.956. This leads to the conclusion that the response surface from Eq. (12) also predicts the ultimate strength very well for axial and central points. Nevertheless, the derivation of the polynomial response surface can reduce the estimation error slightly.
Based on the presented results, the simple empirical formulation in predicting the ultimate strength of cracked plates considering five variables is established. Additionally, the impact of different factors and their interactions is established. It seems that 32 observations are sufficient even in the case of a physical experiment.

COMPARISON WITH OTHER EMPIRICAL FORMULATIONS
A comparative analysis is carried out to investigate the differences between the proposed formulation and other existing formulations as proposed by Paik et al. [7] and Babazadeh and Khedmati [21].
The formula presented in [7] seems to be the most conservative one, which means that all estimates based on it are higher than expected. Firstly, the ultimate strength of the intact plate is calculated based on the formula of Paik et al. as a function of the plate slenderness ratio [36]: (13) Then, based on the formula from [7], the capacity of a cracked plate is estimated as: (14) where A c is a reduced plate cross-sectional area due to the presence of a crack and A 0 is the cross-sectional area of the intact plate.
The second formula as proposed in [21] takes into account the plate slenderness ratio, plate aspect ratio, and crack length. The formula was derived for transversely oriented cracks only. For that reason, in the case of longitudinally oriented cracks, the crack length is taken as equal to 0. In the case of inclined cracks, the crack length is calculated as a projection of the crack in the transverse direction. The formula considers the medium level of the initial imperfections. The ultimate strength of the intact plate is calculated as: (15) Then, the ultimate strength of a cracked plate is calculated based on: (16) where CL is the crack length.
Comparisons between the proposed formula in [7] and other empirical formulas are presented in Figs. 11 and 12. The straight line shows the situation when the estimations from both formulations are equal to each other, whereas particular points are the values of the ultimate strength for the specified values of the design variables. When the selected point is above that line, it means that the estimation from the proposed formula is less conservative compared to the other formulation. In the case of Fig. 11, it is visible that the estimation is very conservative. For almost all cases, the ultimate strength is higher based on the formulation developed in the current study. The mean difference between the two formulations is about 27%.
In the case of the comparison presented in Fig. 12, it can be seen that the deviations are significantly less between the formula developed in the current study and the one proposed in [21]. In general, the formulation from Eq. (12) seems to be slightly less conservative compared to Babazadeh and Khedmati's [21] formula. The mean difference between the two formulations is about 12%. The differences may originate from various sources, such as differences in the Formula developed from the current study vs one proposed in [7] boundary conditions. Additionally, in the case of the second formulation, only the mean level of the initial imperfections is considered.

SIMPLE FORMULATIONS FOR INTACT AND REPAIRED PLATE
The objective of this work was to perform a structural reliability assessment based on the experimentally and numerically estimated ultimate strength.
To determine the difference in capacity reduction during the entire life-cycle of the plate, two additional experiments were carried out, following the same methodology as presented in the section on the DoE model above. The first experiment is related to the intact plate, and the second is related to the plate after repair. In both cases, the FFD analysis is carried out.

INTACT PLATE
In the case of the intact plate, only three variables need to be considered, which are the plate aspect ratio, plate slenderness ratio, and level of initial imperfections. This led to 2 3 = 8 cases which need to be calculated. The ranges of the variables are the same as presented in Table 4. The results of the calculations with the DoE matrix are presented in Table 8.
Further, the importance analysis is performed, and the following factors and interactions are found to be influential: A, B, E, AE, BE. It can be seen that the same effects were found in the case of a cracked plate. Based on the analysis results, a simple formulation is derived. The variables x 1 , x 2 and x 5 will be the same as presented in Eqs. (7), (8), and (11), respectively. The estimated ultimate capacity for the intact plate is equal to: Compared to the exact FE results, the mean value of the estimation concerning Eq. (17) is equal to 1.5%.

REPAIRED PLATE
The repair of crack damage is usually done by butt welding of the cracked edges. Sometimes an additional doubler plate is considered. The welding could cause additional residual stresses in the region around the crack, which may have an impact on the plate capacity reduction.
In the presented study, the simplified distribution of residual stresses acting in the longitudinal direction of the crack is considered. As presented in [37], in the weld zone, there are existing tensile stresses with a value equal to the yield stress. The total width of the heat-affected zone is taken as three times the plate thickness. In the region outside the welding zone, there are existing compressive stresses of about 20% of the yield stress. The width of that region is calculated to satisfy the internal force equilibrium in the plate and is equal to 7.5 times the plate thickness on one side of the crack. The resulting typical zones for the inclined crack are presented in Fig. 13. The residual stresses are applied as the initially induced stresses in the FE model.
Similarly to the cracked plate, the five variables with their ranges, as presented in Table 4, are considered in the experiment. The analysis results, together with the factorial test matrix, are presented in Table 9.  As a result of the importance analysis, the factors A, B, E, AE, BE are found to be the most significant. This leads to the conclusion that the crack parameters and associated residual stresses do not have a significant impact on the structural capacity of the plate compared to the intact one. And only a slight reduction of 1-2% is observed in the case of a 90-degree crack orientation concerning the direction of the applied compressive load. The estimated ultimate strength for the repaired plate is defined as: repaired = 0.8628 + 0.06648x 1 -0.1093x 2 --0.03032x 5 + 0.04153x 1 x 5 + 0.03214x 2 x 5 (18)

CONCLUSIONS
The presented study analysed the ultimate strength of cracked plates subjected to compressive loads considering five governing parameters related to the plate aspect ratio, plate slenderness ratio, crack length and orientation, and the level of initial imperfections. Based on the initial sensitivity analysis, it was found that the existence of drilling holes and the crack position do not influence the plate capacity significantly, and they were not considered in the present analysis. The DoE techniques were used considering the Full Factorial Design. When incorporating the DoE techniques into the FE analysis, one can perform the sensitivity analysis in a quite straightforward way. Compared to the OFAT analysis, which was used in previous studies, the interaction between factors is established, and the number of observations needed in the analysis is significantly reduced. Additionally, this allows more variables to be considered with a similar computational effort.
In terms of the ultimate strength of cracked plates, the plate slenderness ratio was found to be the most influential factor. The plate aspect ratio, crack length, and orientation have a similar impact. Four interaction effects were also found.
As a result of the analysis, a linear response surface was found, which allows for a fast estimation of the ultimate strength of cracked plates. The simplified formulation was found to deviate slightly from the exact results. The developed formulation is a fast and practical tool for estimating the ultimate strength of cracked plates and can be successfully used in the reliability analysis. Additionally, the simplified formulations for estimation of the ultimate strength of intact and repaired plates subjected to compressive loads were developed. These three formulations cover the entire life-cycle of the plate element in terms of possible crack damage.
The comparison between the intact and cracked plates reveals that, for the most critical case of a transverse crack, the reduction of the ultimate strength reaches the level of 35%. In the case of repaired plates, their ultimate strength is only slightly reduced compared to intact plates.