Implementation of Virtual Crack Closure Technique for Damaged Composite Plates Using Higher-Order Layerwise Model

A higher-order layerwise model is proposed to determine stress intensity factors using virtual crack closure technique for singleedge-crack aluminum plates with patch repairs. The present method is based on p-convergent approach and adopts the concept of subparametric elements. In assumed displacement fields, strain-displacement relations and three-dimensional constitutive equations of layers are obtained by combination of twoand one-dimensional shape functions. Thus, it allows independent implementation of p-refinement for in-plane and transversal displacements. In the proposed elements, the integrals of Legendre polynomials andGauss-Lobatto technique are employed to interpolate displacement fields and to implement numerical quadrature, respectively. For verification of the present model, not only single-edge-crack plates but also V-notch aluminum plates are first analyzed. For patched aluminum plate with behavior of complexity, the accuracy and simplicity of the present model are shown with comparison of the results with previously published papers using the conventional three-dimensional finite elements based on h-refinement.


Introduction
Adhesively bonded repairs cause minimum stress concentration of damaged structural components as well as reducing stress intensity factor (SIF) at the crack tip [1,2].The crack patching technique may alter the load path that induces efficient load transfer from cracked component to reinforcement.Thus, the reduction of SIF caused by bonded patch repair retards reinitiation of crack or further crack growth.The properly designed patching leads to the asymptotic value of SIF at the crack tip.It means that SIF at the crack tip with patch repair would become independent to crack length [3].However, the single-side patching may considerably lower the repair efficiency because of out-of-plane bending caused by shift of the neutral plane away from that of a parent plate.Thus, it has been suggested that the SIF in a single-side repair would increase indefinitely with increase of crack length due to influence of bending in the component [4].
From a numerical consideration, conventional finite element analysis has largely been utilized to evaluate the effectiveness of bonded patch repair technology.Although the finite element analysis of components with adhesively bonded patch is a three-dimensional (3-D) problem, twodimensional (2-D) finite element models have sometimes been adopted currently because of cost of analysis, computational time, and difficulties of modeling.The 2D models have usually utilized 2D plane stress element or Mindlin plate element and have sometimes adopted shear springs for modeling of adhesive layer [5].For the analysis of patch repairs, Naboulsi and Mall [6] introduced a 2D finite element using three-layer technique composed of a parent plate, adhesive layer, and patch layer.The effects of disbands were investigated by Ouinas et al. [7] using 2D plane stress element.In spite of that, these 2D modeling approaches still have difficulty of representing variation of each stress 2 Advances in Materials Science and Engineering component through the thickness in a structure stacked with different materials.Thus the models cannot show influence of variation of the SIF through the thickness in cracked plate with asymmetric patch.3D modeling using conventional elements such as isoparametric 8-node or 20-node elements has recently been implemented on some studies by Barut et al. [8] and Ellyin et al. [9].One of weak point in 3D modeling is to require large number of degrees of freedom (NDF).In particular, because thickness of adhesive is relatively too thin as compared to that of patch or parent plate, the mesh configuration without consideration of adhesive thickness would result in large aspect ratios such that might lead to erroneous results.
In this study, higher-order layerwise model is proposed to determine SIF using virtual crack closure technique (VCCT) for single-edge-crack aluminum plates with patch repairs.The higher-order layerwise model based on integrals of Legendre polynomials is proposed to overcome some weak points in 3D modeling such as requirement of large NDF and large aspect ratio due to very thin adhesive.For implementation of numerical quadrature in the proposed model, Gauss-Lobatto quadrature is employed.For verification of the present model, V-notch aluminum plates are first analyzed.Next, for patched aluminum plate with behavior of complexity, the accuracy and simplicity of the present model are represented with comparison of the results with previously published papers using the conventional 3D finite elements based on ℎ-refinement.What is more, the reduction of SIF at crack tip is calculated with respect to the crack length and then the influence of patch size is examined as well.

Higher-Order Layerwise Model
In higher-order layerwise model, three displacements (, V, and ) of a layer in a quadrilateral laminate plate can be expressed as follows: where the repeated subscripts  and  imply Einstein summation rule with respect to plane and thickness.In (1),    , V   , and    refer to modal variables. is 2D shape functions with respect to natural coordinates  and , which have vertex, side, and bubble modes, and  refers to one-dimensional shape functions, which have vertex, and thickness modes, as shown in Figure 1.In the one-dimensional (1D) shape functions, two vertex modes adopt linear Lagrange polynomials as follows: For the thickness modes to more correctly approximate displacement fields, integrals of Legendre polynomials are used as follows: 2D shape functions with hierarchic properties can be built from the 1D shape functions defined above.Four vertex modes are identical to linear Lagrange interpolations and each mode on four sides of elements in any order of shape functions (-level) can be defined as where the superscripts refer to side indexes.Bubble modes are valid for  ≥ 4 only and can be obtained by taking the product   ()  () so that  +  =  + 2 and that both  and , positive integers, are not less than 3.The number of modes per element is dependent on -level.In an element, the number of modes of each group in Figure 2 with respect to different -levels can be presented as From (1), the strain matrix in the higher-order layerwise model can be expressed as Considering a state of anisotropy which possesses three mutually orthogonal planes of symmetry, the following stress-strain relationships in any layer  can be defined: Here [] is a 3D elasticity matrix which is composed of engineering material constants with respect to principal material axes (1, 2, 3), and [Λ] is the transformation matrix based on the angle between layer axes (, ) and principal axes of anisotropic material (1, 2).The elasticity matrix based on the principal material axes can be written in the form where Here  1 ,  2 , and  3 are Young's moduli in 1, 2, and 3 material directions, respectively.Moreover, ]  are Poisson's ratios, defined as the ratio of transverse strain in the th direction to the axial strain in the th direction when stressed in the th direction. 12 ,  13 , and  23 are shear moduli in the 2-3, 1-3, and 1-2 planes, respectively.Also, the transformation matrix [Λ] can be shown to be Here , , and  refer to the principal axes (1, 2, and 3) of the layer material.The variable   means direction cosine of positive  direction with respect to the -direction.
To derive the element stiffness matrix of the higher-order layerwise model, the principle of virtual work is considered.The displacement fields {Φ} of a layer with components defined in (1) can be written by the following general form: Here the matrices {} and {} are in terms of approximating functions corresponding to nodal modes {} and nonnodal modes {}, respectively, for the layer with the reference axes (x, y, and z).The nodal modes only have vertex modes of 2D and 1D shape functions aforementioned and the other modes except the nodal modes belong to nonnodal modes.The nodal modes have physical meaning, while nonnodal modes with respect to the increase in the order of the hierarchical approximating functions do not have physical meaning but improve accuracy of analysis.The element equations for a layer can be expressed by using the principle of virtual work as follows: where  is external virtual work expressed in the following form: In this equation, the superscripts , , , and  signify point forces, side forces, surface forces, and body forces, respectively.Using (11), the virtual displacements can be expressed as and the corresponding virtual strain can then be expressed as where [] and [] are the strain-displacement matrix corresponding to nodal and nonnodal modes, respectively.Using these definitions, the virtual work relation shown in (12) can now be written as Based on (16), the element stiffness matrix of a layer can be expressed as It is obvious that the higher-order layerwise model represented by (17) can account for a layer-to-layer variation of material properties.Also, if there are no gaps and empty spaces between interfaces of layers, compatibility conditions can be applied at the layer interfaces.

Virtual Crack Closure Technique
For determination of SIF at crack tip, VCCT can be employed in which the constitutive model is developed in the framework of linear elastic fracture mechanics, neglecting material nonlinearity.It is assumed that a crack extension of Δ within interval from  to  + Δ does not significantly alter the state at the crack tip.Furthermore, the assumption that when a crack extends by a small amount, the energy dissipated in the process, Δ, is equal to the external work  required to close the crack to its original length is represented as follows: Figure 3 illustrates the modeling before analysis and the deformation after analysis to calculate the energy release rates using VCCT.After analysis, it is noticed that crack extension has occurred.At that time, the work required to close the crack back to its original position by nodal force  is defined by Therefore the total energy release rate can be computed from the nodal forces (   ,    ) and relative displacements ( and V) as follows: where  is the width of the specimen.In 2D crack problems, the total energy release rate is calculated from the individual mode components as Therefore, one can obtain the energy release rate with respect to individual mode with

Numerical Examples
4.1.V-Notched Plates.For verification of proposed model, first, V-notch specimens with 1/ 1− type singular stress fields, in which  refers to distance from corner tip and is a parameter depending upon the abrupt change of geometry or material at corner tip of specimen, are analyzed.In the V-notch specimens, stress concentration near sharp notch is very high.In particular, the peak stress at the notch tip is singular according to theory of elasticity.In general, SIF is very important quantity to be determined in structures Here  and  I (/) are a nominal stress and nondimensional SIF that depends on the considered geometry, respectively.Gómez and Elices [11] report the values of  that is dependent on the V-notch angle and that is well known as the strength of stress singularity.The values of  with respect to opening angle  are tabulated in Table 1.Because the plate shown in Figure 4(a) is symmetry with respect to -axis about geometric configuration and loading conditions, the half of plate is considered for computational domain.Figure 4(b) shows configuration of mesh with 12 elements based on higher-order layerwise model.Generally, the converged SIFs are obtained at fifth or sixth order with respect to -plane and second or third order with respect to thickness.So, in current analysis, seventh order is chosen as the order of shape functions with respect to -plane and third order as the order of shape functions with respect to thickness.The used NDF is 3,996.Numerical predictions in terms of nondimensional SIF  I (/), when /  = 0.4, are shown in Table 2 as the V-notch angle  varies from 0 ∘ to 90 ∘ .Also, the nondimensional SIFs with respect to depth of V-notch are shown in Table 3 when  = 90 ∘ .Agreement between the present values obtained by the proposed model and other solutions in the literatures [10,12,13] is good in both cases.

Unpatched or Patched Single-Edge-Crack Plates.
A plate of thin aluminum sheet with a single-edge crack of length  = 15 mm is considered as shown in Figure 5;   = 200 mm, To check the accuracy of the finite element model used, comparisons of values of nondimensional SIF of the unpatched plate with an edge crack are made with analytical solutions by Gross & Brown, and Tada obtained from the crack handbook [14].It is noted from Figure 6 that the nondimensional SIFs by the proposed higher-order layerwise model show excellent agreement with those in references as the ratio of /  (crack length to width of plate) is increased from 0.1 to 0.8 when the ratio of height to width in the plate (  /  ) is fixed as 5.
Meanwhile, for the values of SIF in cracked plates, the external loadings which most references have dealt with are remotely applied stresses.The condition of remotely applied stress means that the ratio of height to width (  /  ) shown in Figure 5 is not less than 2. When ratio of height to width in cracked plates is less than 2, however, the available values about SIFs of cracked plates have been considerably limited.Thus, Table 5 shows the values of nondimensional SIF obtained from the present model and from Yan's work [15] which used analytical solution, when the ratios (  /  ) are smaller than 2. It is noted from the results that even when external loadings are not remote applied stresses, the present finite element modeling approach also yields accurate results for the unpatch plates with single-edge crack.
Next, the behavior of patched plates is investigated by the proposed higher-order layerwise model.The results of conventional finite elements based on linear analysis by Umamaheswar and Singh [16] are considered as reference values for comparing with the results of the present models.shows finite element meshes based on the proposed higherorder layerwise model.The computational domain takes the half of the entire system by taking advantage of symmetry conditions identically for both meshes.For the computational model of the patch repair configuration, the adhesive is treated as cracked along with the plate.The crack is modeled as a free face.The mesh from Umamaheswar and Singh [16] adopted different levels of discretization for different subregions of the computational domain which are usually used in modeling technique based on ℎ-refinement and then 10,920 eight-node solid elements.The modeling technique is adopted in order to reduce computational resources.The mesh of the present model shown in Figure 7 which employs -refinement is discretized by totally 36 elements of which 20 elements are assigned to aluminum plate, 8 elements are used for patch material, and the number of elements for adhesive film is identical with that of patch material.It is noted that simplicity of modeling can be improved in appearance by the modeling technique based on -refinement.Concerning solution economy with respect to computational resources, the comparison of the reference and the present model is represented in Table 6.It is seen from Table 6 that the present model requires much less computational resources than that presented by Umamaheswar and Singh [16] based on ℎ-refinement.To see solution accuracy of the present model, convergence test about SIF is first implemented for a crack length of  = 15 mm. Figure 8 shows convergence characteristic of the present model with variation of -levels as compared with reference values.It can be told that the values of SIF start to converge from -level = 4 or 5.It is shown that the converged value is similar with reference value.Also, Table 7 shows the comparison of SIFs between reference and the present model for crack length of  = 15, 25, and 30 mm.It is shown that the present model with -level = 5 (NDF = 5,718) yields results that are within ±3% of those reported by Umamaheswar and Singh [16].In addition, the SIFs of the present model are all calculated by VCCT.As we aware of it, the element size Δ aforementioned in Figure 3 is usually nonsensitive to yield SIF values when Δ varies from 1 to 5% of the crack length, and the crack length is less than about 60% of the width of the plate.Also, it is observed that smaller Δ should be used as the crack length becomes considerably large.
In order to see characteristic of behavior of patch repair, first, the patching effect is investigated in terms of nondimensional SIF of patched and nonpatched configurations as the crack length is increased from /  = 0.1 to 0.8 as shown in Figure 9.The characteristic of patched plate with asymptotic values of SIF as a function of increased crack length compares well with the results of previous researches [17,18].It is noted that the values of the nondimensional SIF are virtually unchanged even though the crack length is increased under the condition of suitable patch size.For this behavior, the influence of patch size denoted by ratios of half width of patch to crack length (  /a) has been investigated.Figure 10 illustrates that the increase of patch size beyond   /a = 1.0 has no meaning because nondimensional SIFs are stable.It is recommended that suitable patch size (2  ) requires to be two times the crack length  and the reduction in nondimensional SIF can be approximately 54% when the ratio of /  is 0.375.It is also important to know the influence of patch thickness on the values of nondimensional SIF.When the patch thickness   is the same as the thickness of aluminum plate, that is   /  = 1.0, the reduction in nondimensional SIF is 57% as shown in Figure 11.Also, the nondimensional SIF can be reduced by 69% for   /  = 2.0 in Figure 11.However, the thickness beyond   /  = 1.0 seems to be not practical.
In general, with single patch repair, a variation of SIF across the thickness is expected due to the occurrence of out-of-plane bending action along with primary membrane action.For the patched plate with the crack length of  = 15 mm, the maximum deflection of 0.867 mm is calculated by the present model, while the deflection is zero in the case of unpatched plate.As mentioned before, it means that the behavior of patched plate subjected to uniform tensile loading is dominated by bending moment as well as axial force.So, if the magnitude of uniform tensile loading is considerably large, large transverse deflection would be induced.An increase in the number of elements across thickness within aluminum part is required for obtaining some SIFs across thickness by VCCT. Figure 12 shows the effect on the variation of SIFs across the thickness by the out-of-plane displacement.It is also noticed that the SIFs get larger away from the patched side and the variation across thickness is approximately linear.Also, shown in Figure 12 are the results of the nonlinear analysis in Umamaheswar and Singh's work [16].It is seen that the results of linear analysis of the present model are larger than those of nonlinear analysis, while there exists little difference between the present model and linear analysis by Umamaheswar and Singh [16].However, the difference between linear and nonlinear analysis is not much.Nevertheless, it can be told that membrane effect sometimes needs to be considered for more correct SIFs, especially, in using thick patch.Figures 13(a plots of the in-plane stress   thorough patched, middle, and free sides of the aluminum plate, respectively.They represent well that the reduction in stress concentration around the crack tip is significant in patched plates due to the patching effect, especially, at top surface, that is, patched plane.

Conclusions
This paper presents SIF of rectangular plates with single edgecrack under tensile loading, and the SIF of crack emanating from the V-notch by using the proposed -convergent layerwise model.Some strong points of the proposed model are referred to in the point of view of efficient modeling.In addition, for behavior of patch repair, it is noted that a well designed crack patching leads to an asymptotic value of SIF at the crack tip.The SIF at the crack tip becomes  independent to the crack length when the patch length is the order of two times the crack length.Also, the reduction in SIF can be approximately 54% when the ratio of /  is 0.375.In addition to these, it is confirmed that single-side patches lead to out-of-plane displacement away from the patched side and SIF variation across the thickness.

Figure 1 :
Figure 1: Configuration of modes of 1D and 2D shape functions.

Figure 3 :
Figure 3: Basic concept of virtual crack closure technique.

Figure 4 :
Figure 4: A symmetrical V-notched plate under uniaxial tension and higher-order FE model.

Figure 6 :Figure 7 :
Figure 6: Nondimensional SIF with respect to variation of crack length in unpatched plates.

Figure 8 :
Figure 8: Convergence characteristic of nondimensional SIF with respect to different -levels.

Figure 9 :Figure 10 :
Figure 9: SIF of unpatched and patched plates with respect to variation of crack length.

Figure 11 :Figure 12 :
Figure 11: Reduction in nondimensional SIF with respect to influence of patch thickness.

Figure 13 :
Figure 13: Contour of the normal stresses (  ) in the unpatched and patched plates.

Table 1 :
[18]es of  as a function of the V-notch angles[18].
[10] V-notch for evaluating fatigue strength of the structures.Thus, accurate SIF for the notches with various geometries is necessary to the research.A V-notch plate with isotropic material under uniaxial tension as shown in Figure4(a) is considered.The aluminum plate is in plane stress state with   = 200 mm,   = 40 mm, Young's modulus  = 70 GPa, Poisson's ratio ] = 0.32, plate thickness   = 1.5 mm, and tensile load  = 1.0 MPa.The opening angle  and the depth  of the notch are variable.In general, the SIF of V-notched specimen,   I , can be defined by Chen[10]:

Table 3 :
Nondimensional SIF with respect to depth of V-notch when  = 90 ∘ .

Table 5 :
Nondimensional SIFs with respect to crack extension with variation of   /  .=40, and   = 1.5 mm.In the region of patch repair,   = 13 mm,   = 1.4 mm, and   = 0.2 mm.Patch material is glassepoxy.Material properties in patch part including adhesive film are given in Table4.External loading is a uniaxial tensile load of 3.5 kN giving a remote stress state of  = 58.33MPa.

Table 6 :
Number of degree of freedom used in ℎ-refinement and -refinement.

Table 7 :
Comparison of nondimensional SIF with respect to crack length.