A Cell Centered Finite Volume Formulation for the Calculation of Stress Intensity Factors in Mindlin-Reissner Cracked Plates

In fracture analysis, the stress intensity factor (SIF) is an important parameter which is needed for describing the stress state at crack tip. In this paper a finite volume formulation is developed for calculating the stress intensity factor (SIF) in Mindlin-Reissner plates with a through-the-thickness crack (through crack). For approximating the field variables and its derivatives the moving least square (MLS) technique is utilized. The problem domain is discretized into a mesh of elements where each element is considered as a control volume (CV). The center of CVs are considered as computational points where the unknown variables are associated with. The equilibrium equations of each CV are written based on the stress resultant forces acting on the boundaries of CV where the first order shear deformation theory (FSDT) is implemented in the formulation. Some benchmark problems of plate with through cracks are solved by the present method and the obtained results are compared with the results of analytical and XFEM numerical methods in order to demonstrate the accuracy of the present formulation. These comparisons illustrate the accuracy of predictions of the present solution method. Nevertheless, it is found that the formulation is free of shear locking property which greatly facilitates the cracked plates analysis due to its dual capabilities of analyzing both thin and moderately thick cracked plates.


Introduction
Plate and shell elements are used in the construction of large number of structures such as airplane fuselages, storage tanks, and ship hulls.These structures are subjected to cyclic pressure which may develop defects inside the plate material and lead to arising cracks through the thickness of plate.Due to cyclic nature of applied pressure on these elements, evolving the through crack, results in the sudden failure; so-called as fatigue fracture.Investigation of this type of fracture is an important issue in the design of these structures.On the other hand, the plate element, due to its simplicity is used more than the shell element for modeling of the above mentioned structures.Thus there is a fracture problem that contains a plate element with through crack under the lateral loading in which internal moments and shear forces can be yielded.Similar to cracked plates under the in-plane loads, the main parameter which determines the crack tip field of the cracked plate under the edge moments, is the stress intensity factor.To study the cracked plate under the lateral loading, indeed three dimensional models represent the stress and displacement fields near the crack tip more accurate than the simplified two dimensional theories but accompanied with difficulties in approaching to solution.Therefore, two dimensional plate theories such as the Mindlin-Reissner plate theory -so called the Reissner plate theory-and Kirchhoff plate theory have been developed by researchers to simplify obtaining the solution.The discrepancy between the Kirchhoff and Reissner models is the effect of transverse shear strains through the thickness of plate which is considered in Reissner model while it is ignored in the Kirchhoff model Considering this effect, allows the boundary conditions on the crack surface to be satisfied in the same manner as introduced in the case of three dimensional analysis.However, many researchers analyzed the cracked plate based on the Kirchhoff theory.Besides the analytical methods the numerical approaches have been presented by many researchers to study the bending cracked plate.In the following parts, some of the important works are reviewed.
The primarily works on the cracked bending plate related to deriving the analytical solutions for stress and displacement fields near the crack tip in a cracked plate under the edge moments.For example, in both Refs.[1,2] Kirchhoff theory has been used for a cracked plate under the edge bending moment.In the first work, Williams asymptotic fields near the crack tip have been derived by an eigenfunction expansion, in which the results of the principle stresses in the crack tip were not the same as those of plane stress state, in the second work, Keer and Sve have derived the expressions for the SIF and strain energy calculations at crack tip for edge and central cracked finite plates.Knowles and Wang in 1960 [3] have investigated bending of Reissner cracked plate with infinite dimensions under edge moments.It was found that contemporary to Kirchhoff theory, the principle stress near the crack tip is the same as that value expected from the plane stress state Also the shear stress near the crack tip varies asymptotically as  −1/2 when  → 0, different to that predicted by Kirchhoff theory.Hartranft and Sih in 1968 [4] investigated the effect of plate thickness on the stress distribution field in a very closely region around the crack tip for a Reissner cracked plate in both cases of a crack of finite length and of infinite length.Authors reported the bending stress as much as sixty-two present greater that obtained by Knowles and Wang [3] for very thin plate.Boduroglu and Erdogan in 1983 and also Joseph Erdogan in 1991 [5,6] have used a method based on the Mindlin-Reissner theory for calculating the SIFs in cracked plate under the edge bending moments.In the first work a cracked plate solved by a perturbation method and the SIFs have been calculated for the various crack length to plate width ratio for cases of central crack and symmetric edge cracks [5].In the second work this problem solved by considering the plate thickness to crack length ratio as the effective parameter and the obtained results were differed considerably with those values obtained by Kirchhoff theory [6].Zucchini et al. in 2000 [7] implemented a three dimensional finite element analysis for thin cracked plates in which the applied loading includes the edge bending moments, edge twisting moments and also out of plane shear forces.In That work, the results were compared with the analytical results presented in literature and obtained using the Reissner or Kirchhoff theories.Reported results denote that the stress field approximated by Reissner theory is conformed to three dimensional elasticity field in distances less than 0.1 of plate thickness from crack tip only, and is diverged away from these distances.On the other hand, the stress fields approximated by Kirchhoff theory is conformed to ones obtained using three dimensional elasticity theory only in distances more than thickness of plate from the crack tip.However authors reported that for thin plate, the value of J integral calculated using Reissner theory has been determined as the average values of this integral obtained by the three dimensional theory.Also authors reported that energy release rates computed based on the Kirchhoff and Reissner theories are the same as values computed using the three dimensional finite element analyses.Hui and Zehnder in 1993 [8] have derived the stress and displacement fields near the crack tip for a Reissner cracked plate using an asymptotic approach.Also the authors calculated the SIFs for Reissner and Kirchhoff thin cracked plate and produced.Also the authors suggested, for very small plate thickness in the order of 0.02 crack length, as long as small scale plastic yielding tip is prevailed on plate thickness, the Kirchhoff approach may be a better choice for solving the problem of cracked plate.But Zehnder and Viz in 2005 [9] proposed the analysis of thick cracked plate under the edge moments and torsional loading must using the Reissner theory or even using the three dimensional elasticity theories.Wang et al. in 2003 [10] have derived the complex stress function for a single crack in an infinite plate, the moment stress intensity for cracked Kirchhoff plate were obtained by satisfying the boundary conditions of crack surface and plate edge in the mentioned function using of the collocation method.
Analytical solutions for finding SIFs for cracked plate under moment and shear loading have been limited to very simple geometries and loadings.Even in Reissner cracked plates, the most results have concerned for infinite plates and a few solutions have been presented for finite ones.Therefore considerable attempts have been made based on the numerical methods in order to fill this gap.Among the numerical methods, the finite element method as a prevalent method in fracture mechanics has been used to study the bending of cracked plate under the lateral loading and edge moments by many researchers.Alwar and Nambissan in 1983 [11] by using the isoparametric elements and a three dimensional modeling performed a fracture analysis of thick plate and found that for this type of plate, the SIF is varied nonlinearly across the plate thickness.Barsoum in 1976 [12] for simulating the singularity in the crack tip, have used a degenerated singular isoparametric element at crack tip.The author calculated the SIF for a central cracked plate with different thickness covering both thin and thick plates.Accurate values obtained for the thick and moderately thick plates but for thin plates a substantial error has been reported compared to Reissner theory.Ahmad and Loo in 1979 [13] proposed a triangle singular element for crack tip to calculate the bending and shearing SIFs in thin plates.In this method, a refined mesh was also used near the crack tip and several problems of cracked plate with edge moments and lateral pressure loading for central and edge cracks have been analyzed.The results obtained by using of a displacement type formulations which deduced from the Williams expansions [1].Alwar and Ramachandran in 1983 [14] by using the three dimensional finite element analysis in which the isoperimetric singular elements have been suited in crack tip, calculated the SIF values in a central cracked plate for various dimensions of plate and crack length.The obtained results found to vary nonlinearly across the plate thickness and are as much five to ten percent more than the results obtained from the Reissner theory.Rhee and Atluri in 1982 [15] have presented a hybrid stress finite element procedure wherein the singular quadrilateral element was used in crack tip and the SIF are calculated for cracked plate based on Reissner theory.Sosa and Eischens in 1986 [16] have derived the path independent integral socalled as J-integral for the cracked plate based on the Reissner theory.The authors by using a simple finite element mesh conformed to crack surface and also using the path independent integral, the SIFs for the cracked plate computed based on Reissner theory and various crack length to plate dimension ratios were considered.Dolbow used the extended finite element method (XFEM) [17] for the fracture analysis of Reissner plate in which the enrichment functions obtained from the analytical solution for crack tip fields used for the crack tip elements.In the above work, a plate with an inclined centered crack with various angles under uniform bending was analyzed successfully; however the shear locking problem for very thin plates has been reported in this Reference.Su and Leung in 2001 [18] by expansion the displacement field near the crack tip and also employing the fractal two-level finite element, both moment intensity and twisting intensity factors have calculated for a bending mode and mixed mode loadings of a center cracked plate based on Reissner plate theory.The results obtained were in agreement with those calculated by other numerical methods presented in literatures.Recently Lasry et al. have used XFEM to calculate the SIFs for Kirchhoff plates [19].They computed SIFs for a very thin cracked plate using two strategies, the direct estimation and J-integral.They also obtained the optimum value of the radius of enrichment domain which have been derived in terms of crack length.Bhardaj et al. calculated SIFs for the Reissner cracked plate with different loads and boundaries using the extended isogeometric analysis [20].Their obtained results were in agreement with those have been obtained from XFEM even were more accurate in some cases.Tanaka et al. in 2015 [21] have been developed an effective meshfree formulation for cracked bending plate based on Reissner theory in which the stabilized conforming nodal integration and subdomain conforming integration techniques have been used to integrate the stiffness matrix.For several through the thickness cracked plate problems, the moment stress intensity factors evaluated using the J-integral approach, have been reported and compared to that obtained from the analytical solutions.Also Tanaka et al. in 2017 [22] the meshfree formulation employed and mixed moment intensity factors have been evaluated for an inclined cracked plate based on Reissner theory via a path-independent J-integral and accurate results have been reported compared to analytical solutions.It should me mention that in the modeling of the cracked plate a very refine mesh is used in a small region around the crack tip.
According to the above works there are issues in the application of finite element method for the analysis of cracked plate.In order to simulate the crack tip singularity in the framework of finite element, one must use the singular elements such as the degenerated isoparametric element with a very refined mesh around the crack tip.Also, XFEM the well-known finite element method for fracture analysis of structures, cannot be assessed as a reliable method for this type of fracture problems due to encountering with the shear locking problem in the analysis of very thin plates based on the Reissner model.So one may seek an alternative method that can overcome these issues.
The finite volume method (FVM), which is a well-known numerical method in the field of fluid dynamics, has been also used for the stress analysis of structures [23].Wheel [24] presented a finite volume formulation for the analysis of bending plates based on Reissner theory.Author declared that no shear locking observed in the analysis of very thin plate by this method.Also Fallah [25] formulated the cell-centered and cell-vertex finite volume methods for the Reissner plate analysis under the lateral loading where no shear locking has been reported.Ivankovic et al used the FVM to study the crack propagation phenomenon for the analysis of pressure pipe test [26].Also Stylianou et al [27,28] applied a FVM based method for studying the dynamic fracture in 2D media by using the methodology of release node technique and cohesive zone concept.According to the published research, no investigation for the application of FVM for the fracture analysis of cracked bending plate has been reported so far.Therefore, in this work, a finite volume method based formulation is presented to study the bending of a cracked plate.Simplicity in implementation and the inherent high stability and also being shear-locking free in the Reissner based plate formulation are those of the attracting feature of the FVM method [23,24].In this work, a local weak-form equilibrium equation of isotropic moderately thick plate, based on the Mindlin-Reissner plate theory is obtained.Moving Least Squares (MLS) approximation which has been used in many of the numerical methods such as, the meshless local Petrov Galerkin method [29] and the element free Galerkin method (EFGM) [30] has been used in this work.MLS technique is able to approximate the displacement field with a desired order of accuracy by using appropriate polynomials in its basis function [31].A regular mesh consist of quadrilateral control volumes (CVs) over the plate domain is constructed and the visibility criterion [30] is used to treat the crack line presence in interpolation of the field variables of nodes around the crack line.Also to enforce the boundary condition at plate edge and crack surface, a set of nodes suited on the boundary of plate and also suited adjacent to crack line are used.In the mentioned nodes the boundary condition of plate and the condition of crack surface are imposed by the collocation technique [32].This treatment of the discontinuity due to crack is more easy to accommodate than what is done in XFEM in which a Heaviside function imported into the finite element mesh and unknown variables defined for nodes located nearing the crack geometry.Also in this research the interaction integral is employed for evaluating the moment intensity factor in which appropriate auxiliary terms of stress and displacement are needed and these terms can be selected from the expansions of stress and displacement fields near the crack tip presented in literatures.Several problems of through the thickness cracked plates with thickness ranging from thin to thick and under different loadings are solved using presented finite volume based method.The results obtained are more accurate than those obtained by other numerical methods.The moment intensity factor calculated for varying ratios of crack length to plate width.Versace ratios of crack length to plate width.
The layout of the present paper is as follows: in the next section, a plate formulation according to the Mindlinreissner model is presented briefly.Then by applying the weak form of differential equations governing the plate equilibrium and considering the boundary conditions especially boundary conditions of cracked faces, the FV based discretized equilibrium equation is obtained for all the CVs representing the plate.In section 3, crack tip fields for cracked Reissner plate and also the interaction integral for calculating the SIFs values of crack tip are presented.In section 4, by solving some benchmark problems, the results of the present study are compared with the results obtained by other numerical and analytical methods.And finally the conclusions are drawn in section 5.

Formulation
In the finite element method the mid-plane of plate is discretized to a set of elements where by interpolating the displacement fields inside each individual element and applying the energy principles the discretized equations are obtained.On the other hand, in the finite volume method the CVs constructed around the field nodes scattered over mid-plane of plate, and the equilibrium equation is written for each CV.By assembling the established equations of all the CVs, the system equations governing the problem domain can be resulted.The distinguished feature of the finite volume method is satisfaction of the equilibrium equation for each local domain of CV when the domain integral is eliminated, hence researchers found this method as an interesting method in the computational mechanics [33,34].In the finite volume method, the CVs can be taken overlapped or non-overlapped [35].On the other hand to approximate the displacement function at an interest point of the domain one can employ the well-known techniques for approximation of displacement which are used in the meshless methods, like Radial point interpolation method (RPIM) and moving least square method.These techniques already have been applied to the analysis of plates successfully [36][37][38].Therefore, by employing these techniques in the finite volume method the displacement of an interest point on the faces of CVs can be interpolated in terms of displacement of nodes fallen inside the 'support domain' defined for that point.Among the later techniques, the MLS is a well-known method for approximating the unknown variables in any point of domain where by constructing a smoothed function over a set of nodes the nodal displacements and rotations can be resulted [31].
In this section, at first a finite volume based formulation for the Reissner plate with a through crack is presented and then the MLS technique for the approximation of field variables at the interest point of mid-plane of plate is presented.

Finite Volume Formulation for Cracked Reissner Plate
The Mindlin-Reissner plate is a two dimensional model for the three dimensional body of the plate where one of its dimension, i.e. thickness, is small rather than the two other dimensions.In this model, different to Kirchhoff model the main assumption is that the plane perpendicular to the mid-plane of the un-deformed plate is not perpendicular to the deformed mid-plane after deformation due to the presence of the lateral shear stresses [39].Accordingly, at every point of interest, P, the plate section rotations about the x and y axes of Cartesian coordinates, i.e.   and   deduced as follows: where w is the transverse displacement, β x is rotation about y axes and   is the rotation about the x axes corresponding to the considered point P, Figure 1.In this way, the displacement components of interest point, P, of the mid-plane can be presented as follows: Where u and v are the in plate displacement along the x and y axes respectively.
As displayed in Figure 2 to analyze the cracked plate using the finite volume method one can discretize the midplane of plate to a set of CVs which are corresponded to the field nodes distributed over the domain.In the present work, a uniform node distribution, hence identical rectangular or square CVs, depends to the plate geometry, can be constructed.Although one can use the irregular node distribution where results in CVs with different geometries.It should be noted that in the cracked plate, the CVs are considered on both sides of the crack line as shown in Figure 2. Also it should be mentioned that in the present work, we use non-overlapping CVs.In FVM, the governing equations for cracked plate can be obtained by writing equilibrium equations of acting external transverse forces and internal stress resultant forces on faces of a CV, see Figures 3 and 4. Therefore the equilibrium relations for a CV placed on the mid-plane of plate lied on the xy plane, can be written as follows: ∑   = 0 (5) ∑   = 0 Where Equations 5 and 6 show the balance of moments about the y and x axes respectively and Eq.7 represents the balance of forces acting along the z axes which is perpendicular to the mid-plane.According to Figure 3 in which a typical CV with an arbitrary number of faces is shown, the force and moments corresponding to the faces of CV which are resulted from the normal and shear stresses distributed over the plate thickness can be integrated on the CV faces, These integrals can be approximated by using Gauss quadrature rule, which is applied in the present work.Hence, the above equilibrium equations can be approximated as follows: The parameters presented in Equations 8 to 10 are defined as follows: n k is the number of faces of CV, n g is the number of gauss point considered on each face, X Qj k is the coordinates for gauss point j-th on face k-th, α x k and α y k are cosine directions of outward normal to face k, w ̂j and J qk are weight and Jacobin values of the integration in the j-th gauss point of the k-th face, respectively.It should be noted that if the face of CV laid on the crack path it is free of traction so no Gauss point is defined there.Equations 8-10 expresses the relations of internal moments and shear forces acting on CV faces and the external loading on the plate.One can use the constitutive equation of Reissner plate [39] to express the moments and shear forces in terms of strains components, this equation is presented as follows: . [ Where  is the flexural rigidity and defined by: In which t is the plate thickness, E, G, and ν are elastic modulus, shear modulus of plate material and Poisson's ratio, respectively.In Eq.11, the strain vector can be divided to two parts as, bending strains, ε b , and shearing strains, ε s , which are expressed in terms of the displacement and rotations derivatives as follows [39]: Where β x and β y are rotations at the interest point on the faces of CV.Hence Equations 8-10 can be expressed in terms of unknown displacement and rotations by using the above equations.

Moving Least Squares (MLS) Approximation
In this section the MLS technique for approximating the field variables of the Reissner plate, namely rotation components and transverse displacement, is explained.Since in MLS approximation for field variables of the plate, similar interpolation function can be used for all the variables, hence the technique is presented for a general field variable say u.
In the MLS technique the unknown function u(X) at a favorite point X of the domain of problem is approximated continuously as follows [31]: Which m is the number of the monomial basis used in the MLS approximation, and p T (X) is the vector of monomial basis which is a linear combination of variables 'x and y' with the desired orders.The first-order and second-order of p T (X) for a two dimensional domain are stated as follows: For linear basis, And for quadratic basis, we have: () = {1, , ,  2 , ,  2 } (16-b) Note that usually the quadratic basis is selected as basis function in application of MLS to plate problems [31].In Eq. ( 15), a(x) is called as the associated unknown coefficients and x = [x, y] T is a vector of space coordinates.The a(x) coefficients are obtained by minimizing a weighted discrete L2 norm as follows: Where n is the number of nodes fallen into the support domain of interest point X.Support domain may be assumed in any shape, for example, circle, rectangle, or elliptical, where the rectangular shape centered at the interest point X is used as the support domain in this work, Figure 5.Note that if crack geometry crossed the support domain of interest point, some nodes fallout from the support domain.This issue can be treated using techniques such as visibility criterion or diffraction method [40].
Which α s is a parameter depending on the average field node spacing and its appropriate value for each problem is determined by numerical experiments [31].
In Equation 17, w is the weight function, which is equal to unit in the vicinity of point X I and equal to zero on the edges of support domain.The conventional forms can be used for weight function such as the cubic spline function, quartic spline function, or exponential function.In this work the quadratic spline function is used for the weight function which is defined as follows: In Equation 20, |X − X I | is the distance from node I to the interest point X, and r s is the size of support domain for node I which is equal to d sx and d sy in x and y directions, respectively.
In Equation 18, u ̂ is the vector of unknown parameters associated to field nodes, also A(X) and B(X) called as coefficients moment matrices which are expressed by the following equations: By solving the Equation 18, a(X) can be obtained as follows: Substituting Equation 23in Equation 15 results in the following equation for u (X): Where ϕ I is the interpolation function which is expressed as follows: () = ∑   ()  =1 ( −1 ()())  = p T (X)( −1 ()())  (25) According to Equation 25, displacement and rotation components of plate i.e. w, β x , and β y are interpolated using MLS technique as follows: Where superscript h denotes the approximated values of variables associated with the displacement and rotation components of field nodes.

Essential Conditions
In the finite volume method the essential boundary conditions can be enforced by using of line CVs which are placed on the boundaries of plate [25].As shown in Figure 6a. to perform this task, the nodes associated to line CVs are placed on the essential boundaries of the plate.For example for plate with simple support edges or fixed support edges, the displacement value of boundary node, say node B in Figure 6 considered as a line CV should be set as follows: For fixed support: And for simply support: Where as seen in Figure 6b,  (36) For example, for any field node on the free edges of a plate, the moments and shear force can be expressed as follows: = 0 ,   = 0 ,  = 0 (37) It should be mentioned that in case of a mixed boundary conditions in which the essential conditions and the natural conditions are applied together, e.g. for simply supported plate, a set of three equations is selected from Equations 31 to 33 and 34 to 36.
Another issue that can be raised is the boundary conditions on crack surface.Since the surface of crack is assumed to be free of tractions, so for the crack surface, the natural conditions corresponding to the free edge are imposed.To accomplish this need, line CVs are considered along the crack line, Figure 6c, and the natural boundary conditions, Equations 34 to 36 relevant to crack surface are imposed to the field nodes considered at the center of line CVs.Also it should be noted that by substituting Equations 11 to 14 in Equations 34 to 36, the natural boundary conditions can be expressed in terms of the field variables β xB , β yB , and w B .By substituting the approximated values for field variables from Equations 26 to 28, the equations of natural boundary conditions are obtained in terms of unknown parameters of neighboring field nodes.
By assembling the equilibrium Equations 8 to 10 for all the CVs together with Equations 30 to 41 for boundary conditions, a system of linear equations is obtained in terms of unknown parameters associated to transverse displacement and components of rotations corresponding to the field nodes.

Crack Tip Fields for a Through Crack in Reissner Plate
Because the plate model is a simplified model of a three-dimensional solid, one can derive the asymptotic fields for stress and displacement near the crack tip of a plate by modification of the fields corresponding to three dimensional model.Hartnaft and Sih [41] derived an Eigen function expansion for general stress and displacement fields in three dimensional cracked problems.Sosa and Eischen [16] picked the singular terms of these functions to employ in the Jintegral for calculating the SIF in the cracked Reissner plate.Dolbow [17] extracted all of the singular terms in the mentioned expansion and also the asymptotic stress and displacement fields near the crack tip in the cracked Reissner plate.The above mentioned findings corresponding to near-tip fields for kinematic variables and also for the resultant moments and shear forces on transverse section of the plate are presented as follows: (38-a) And for moments and shearing force fields we have: )) (39-a) )) (39-c) In Equations 39a to 39e K 1 is the SIF of crack tip corresponding to the case in which the cracked plate loaded by the edge symmetric moments and also K 2 , and K 3 are SIFs corresponding to the cases in which the cracked plate loaded by the edge twisting moments and the edge out of plane shearing force, respectively, see Figure 7.In fact the SIFs for a cracked plate are defined as the limiting values of forces and moments near the tip of crack as follows:

Interaction Integral
Different numerical methods have been used for calculating the SIF value of plate with a through crack under the edge moments.Wilson and Thompson [42] calculated the SIF value for a central cracked plate under the edge moments by substituting the calculated nodal displacement from the finite element method in the Kirchhoff solution for the cracked plate.They used a refined mesh around the crack tip in order to evaluate the SIF value with more accuracy.Ahmad and Loo [13] used a typical singular element adjacent to the crack tip and applied the Kirchhoff solution for calculating the bending and shearing SIFs of cracked plate.Rhee and Atluri [15] presented a hybrid finite element method to analyze the Reissner plate in which a singular element at crack tip has been used and the SIF directly obtained from the singular function corresponding to the crack tip.Sosa and Eichsen [16] applied the path independent J integral for calculating the SIF value in cracked Reissner plate under the edge moments.In order to achieve an exact solution, the authors used a fine mesh around the crack tip.
In this study the domain form of the interaction integral which has been configured for Reissner plate in Ref. [17] has been used.In that reference the domain form of the interaction integral for cracked plate has been concluded from the J-integral.The authors of that work recommended the interaction integral as an efficient technique for the accurate prediction of SIF value.Also recently, the formulation of interaction integral has been applied to Kirchhoff plate by Lasry et al [19].
In Ref [17] the relation between the SIFs, i.e.K 1 , K 2 , and K 3 for a Reissner cracked plate and interaction integral has been presented as follows: Where I is the value of interaction integral and E, μ and t are elastic modules of plate material, shearing modules of plate material and thickness of plate respectively.Also, the interaction integral for a cracked plate based on Reissner theory has been concluded as follows [17]: Where indices i and j are selected from the variables x or y, and M and Q denotes internal moments and shearing forces acting on lateral plane of the plate respectively.In Equation 42, β and w are rotation and transverse deflection of the mid-plane respectively.Also in the above equations, q is the unit-step function which is defined in such a way that, takes the unit value over the central region of the integral domain and takes zero value in borders of the domain, see Figure 8a.Also q ,j denotes the derivative of function q relative to j, δ 1j is the Kronecker delta which is defined as follows: And W int is called the interaction strain energy which is calculated as follows: int =  ∶   aux +  .  aux =  aux ∶   +  aux ∶   (44) Where ϵ b and ϵ s are the portions of bending strain and shearing strain from the total strain at the interest point of plate respectively, which are defined as follows: Where ∇ denotes the derivative operator defined as follows: The domain of interaction integral encloses the crack tip and can be assumed in any shape, where in the present work this domain is assumed as a square area centered at the crack tip as shown in Figure 8 (b).
In Equation 45 variables superscripted by 'aux' are related to the auxiliary state and depends to the considered cracking mode, and they are calculated from the corresponding crack tip fields using Equations 38a to 39e.To evaluate the SIFs values for a cracked plate in the case of mixed modes loadings one can use the proper terms for the interest mode from the Equations 38a to 39e and substitutes the auxiliary superscripted terms in Equation 42.For example for calculating K 1 , auxiliary terms in Eq. ( 42) are replaced by the coefficients of K 1 picked from Equations 38a to 39e, and in Equation 44 the SIFs should set to values , K 1 aux = 1, K 2 aux = 0, and K 3 aux = 0, which in turn the following equation is derived for K 1 : A similar procedure can be followed for evaluating two other SIFs, i.e.K 2 and K 3 .

Numerical Studies
In this section, the calculation of SIFs factors for some of the benchmark cracked plate problems by the use of the presented finite volume method is performed.The solutions are compared with the results of the analytical and other numerical methods presented in references.In all of the considered problems the crack is assumed through the thickness of plate and crack surface is assumed free of any traction.Also the plate material is assumed to be homogeneous and isotropic.

Plate with a Central Crack Under Edge Moments
In this problem, the SIF values of a plate with a central crack under the uniform bending moments applied on two edges which are parallel to the crack line are investigated by using the present method, Figure 9.This problem has been also studied in Refs.[6] and [17] by applying analytical method and XFEM numerical technique, respectively.The geometry of plate is assumed similar to that presented in the later reference which reads, the length of plate, l, equal to 10 mm, the width of plate, b, equal to 10 mm and plate thickness is t where several values are considered for the thickness to represent thin to thick plates in the analysis.The plate material reads, the modules of elasticity equal to 1000 Mpa, and the Poisson's ratio equal to 0.3.The crack half-length, a, reads, 1 mm.To solve the problem by the present finite volume method the plate domain is discretized to a set of uniform square CVs.For enforcing the boundary conditions on plate edges and also on crack surface further nodes distributed aligned with the edges of plate and crack line, see Figure 10.For the nodes located on crack surfaces the free surface conditions have been imposed and the nodes located on edges of plate according to state of the plate edges, the corresponding boundary conditions are applied.For construction the shape function in the MLS approximation a quadratic monomial basis is used.The support domain is supposed as a rectangular shape with dimension equal to 2.4 times the nodal spacing.For evaluating the internal forces corresponding to the faces of CVs, two quadrature points are used.
The domain size of the interaction integral varies in the range of 0.3 to 0.45 of the half-crack length depending to the thickness of the plate.Several plates of different thicknesses with the above mentioned properties are considered in the analysis for calculating their SIF values.In Figure 11 the obtained values of K 1 versus the ratio of plate thickness to half crack-length are displayed.It is observed that for small values of the thickness ratio the results of the present method are in close agreement with the exact values governed by the analytical method presented in Ref. [6] where XFEM predictions are different for thin plates.The lack of accuracy in XFEM predictions for thin plates is due to the shear locking phenomena which arises in the thin cracked Reissner plates analysis using XFEM [17] where such locking phenomena is not appeared in the present method.The capability of the present method for the analysis of cracked Reissner plates with small thickness ratios is also illustrated in Figure 12 where the normalized value of K ̃1 is presented in which K ̃1 = K 1 M 0√ a .As can be seen, the FVM predictions are in close agreement with the analytical results where XFEM numerical method has not shown such accuracy.It can be observed that for both of the present FVM and XFEM, for plate thickness of about 0.1 of the crack length, the accuracy of results is in good agreement.However, when the thickness of plate decreases to values smaller than the 0.1 of the crack length, the accuracy of SIF obtained from XFEM decreases considerably whereas the accuracy of FVM results is not affected.
To investigate the accuracy of the FVM in the prediction of K 1 of the cracked plates with different crack length, several tests have been conducted.In each test, by assuming a plate with constant thickness ratio, t b ⁄ and considering different crack length ratio, a ̃= a b ⁄ , the normalized value of K ̃1 is obtained by using the present method.The obtained results are compared with the predictions of XFEM and with the analytical predictions as shown in Figures 13 to 16.As can be observed, corresponding to the cracked plates with the considered thicknesses, the results of the present method are comparable to the analytical predictions due to the shear locking free feature of the method.
Noticeably, the results of normalized moment intensity factors obtained using the presented FVM are more accurate than those obtained from XFEM when the plate thickness decrease to values less than 0.1 of plate width.Also in Figures 15 and 16 it is observed that the error of SIF values obtained from XFEM exceeds 7 and 9 percent of the analytical values respectively.The poor values of SIF obtained from XFEM as reported in Ref. [17] is occurred due to the shear locking phenomenon which have been observed in all of the finite element based methods.

Simply Supported Plate with a Central Crack Under Uniform Pressure
In this test problem the K 1 value for a square plate with a central crack which is under the transverse uniform pressure is investigated by using the present method.The plate's edges are assumed simply supports as displayed in Figure 17.This problem also has been studied by using of FEM and XFEM numerical methods in Refs.[16] and [43], respectively.The plate geometrical characteristics are taken as, the length equal to 20 units and the thickness equal to 1 unit.The material properties of the plate are taken as, the elasticity modulus of material equal to 1000 units, and Poisson's ratio equal to 0.3.Pressure loading is taken equal to 1 unit.A finite volume model constructed using a 40 × 40 mesh of elements where each element is considered as a CV and a computational node is allocated to the center of each CV.In addition to this, in order to apply the boundary conditions some nodes are placed on the plate boundaries.By doing this a total of 1681 nodes are uniformly distributed over the plate domain.The Parameters needed in the analysis of this problem by the present method are justified as follows: the number of monomial basis used in MLS approximation equal to 6.The support domain assumed as rectangular shape with the same length and width values equal to 2.4 times of nodal spacing.To evaluate the internal forces on faces of CVs, two Gauss points are used.Also the value of the domain size used to evaluate the M-integral is dependent to crack length and varied in the range of 0.05 to 0.65 times the crack length.
In Figures 18 and 19 the contours of bending moment M x and lateral deflection w of a plate in which the crack length is equal to 8 units are displayed.In Figure 20 the normalized values of K I versus the crack half-length to plate width ratios are shown.The results are compared with the predicted results presented in Refs.[16,20,43] where the finite element method and XFEM are used, respectively.It can be observed that the results obtained from the present method have a good agreement with the values predicted by XFEM and FEM.

Conclusion
In this work, a finite volume based formulation has been developed for calculating stress intensity factor of bending plates with through cracks.The capability of the finite volume method for the analysis of solid mechanics problems especially for the bending analysis of plates has already been discovered by the researchers where some of these studies have been presented in the reference section of this paper.In the present formulation Reissner plate theory has been considered which enables the shear effects to be accounted in the analysis.Also, the MLS technique has been employed for approximating the displacement variables and their derivatives at the desired points of the plate using the displacement values corresponding to the field nodes.Using the MLS technique enables to approximate the field variables and their derivatives with a desired order of accuracy which are appeared in the discretized governing equations.Based on the presented solution procedure, some benchmark problems have been solved and it has found that the present finite volume based formulation is able to predict accurate results.Also, it has demonstrated that the method can deal with the both thin and moderately thick cracked Reissner plates.This dual capabilities of the present formulation is due to its shear locking free feature which has been already observed in the un-cracked plate analysis as shown in the works published by the second author.Exploring the capability of the finite volume method for the fracture analysis of bending plates is a new research work.Surely future findings reveals more capabilities of the finite volume method for such challenging engineering problem.

Figure 1 .Figure 2 .
Figure 1.Displacement of the interest point P in the Reissner plate model

Figure 3 .Figure 4 .
Figure 3.A typical CV with through crack aligned with one face

Figure 5 .
Figure 5. Definition the support domain for the interest gauss point on the faces of CVThe dimensions of rectangle support domain in x and y directions shown as, 2ds x and 2ds y are defined as follows[31]:d sx = α s d cxand d sy = α s d cy(19)

Figure 7 .
Figure 7. Cracked plate under the different loads

Figure 8 .
Figure 8.(a) Interaction integral domain, (b) the shape of unit-step function used over the interaction integral domain

Figure 17 .Figure 18 .Figure 19 .Figure 20 .
Figure 17.Simply supported plate with a central crack under the uniform pressure =   cos   +   sin   Where β xB and β yB are components of rotation of interest point B corresponding to the global coordinates, x and y.It is obvious that one can use Equations 26 to 28 in Equations 31 to 33 for expressing the field variables, β xB , β yB , and w in terms of the corresponding values of the neighboring field nodes bounded by the support domain until the final equation for essential boundary condition concluded in terms of nodal parameters of displacement and rotations.In this case, the value of tractions on the boundary of plate which are expressed in local coordinate nt can be transformed to the corresponding values in global coordinate xy as follows:   =    2   +    2   + 2      (34)   = −      +       +   ( 2   −  2   ) (35)  =     + nB and β tB are components of rotation of interest node B in local coordinate n-t which n and t denote respectively the outward normal and tangential to the boundary at the interest point B. According to Figure6bone can express the local components in terms of global components using the following equations: , aux +     , +    , aux +  aux  , ] −  int  1 } , dA +