Variable Singularity Power Wedge Element for Multilayered Composites

Multilayered composites can be weakened by local failures, which can give rise to singular stress fields. Singularities can be removed or integrated to obtain bounded quantities which may be used for the analysis. In any case, the participation of singular stresses should be identified. Serving this purpose, a new mixed singular wedge element based on interpolation functions with a variable singularity power is developed, which can adapt to the problem and can represent the non-singular fields as a particular case. The displacements and the interlaminar stresses are assumed as nodal d.o.f. to a priori fulfill the kinematic and stress contact conditions at the material interfaces and the stress boundary conditions in point form. The application to sample cases taken from literature shows that the element can be successfully employed for the analysis of singular fields. The variable singular representation and the mixed formulation give more accurate results than their displacement-based, non-singular, or singular counterparts with a fixed singularity. The failure loads of initially delaminated specimens are accurately predicted either using the virtual crack closure technique or a mesoscale model by running in a reasonable time on a personal laptop computer.


Introduction
Multilayered composites are increasingly finding use, owing to excellent strength and stiffness-to-weight ratios, long fatigue life and advantageous energy absorption properties.As their performances can be weakened by local failures (e.g., fibre breakage, intralaminar and interlaminar matrix cracking, fibre/matrix debonding, delamination), an accurate prediction of their residual strength and stiffness is mandatory.A recent comprehensive review of the mechanisms of damage formation and evolution and of their modeling is given by Ivañez and Sanchez [1], Hassan et al. [2], Tay et al. [3], Càrdenas [4], Liu and Zheng [5] and Hongkarnjanakul et al. [6].
Damaged areas require a full three-dimensional analysis, since the out-of-plane stresses, usually negligible elsewhere, have a significant bearing on keeping equilibrium in these regions and determine the stress singularities rise at matrix and delamination crack tips.As extensively discussed by Ebrahimi et al. [7] and Koguchi and da Costa [8], a different order of the singularity power is determined by the geometry and the constituent materials.The singular stress fields can be (i) removed or (ii) integrated to arrive at bounded quantities which can be physically interpreted and used for the analysis of progressive damage evolution and failure (see Berthelot [9] and Liu and Zheng [10]).Case (i) is that of continuum damage mesomechanics (CDM), while case (ii) is that of fracture mechanics (FM).Approach (i) is easy to implement into the finite elements and treats crack initiation and growth within a unified model.Approach (ii) is nowadays widely employed to predict the delamination growth using the virtual crack closure technique (VCCT) to solve the energy release rates (ERR).As examples, the papers by Okada et al. [11], Leski [12] and Liu et al. [13] are cited.If not experimentally determined, the critical ERR are computed by criteria such as those by Benzeggagh and Kenane [14], Wu and Reuter [15] and Reeder et al. [16].
Usually, finite element models are used as structural models, since closed form approaches should be limited to the analysis of certain range of geometries and loading conditions, due to the simplifying assumptions adopted, while shear-lag approaches, variational methods and stress transfer mechanics use parameters that should be determined by experiments.
However, the presence of singularities should be recognized, their nature appreciated and their stress fields accurately and efficiently predicted, as extensively discussed by Sinclair [17] and Hu and Yao [18], since the participation of the singular stresses must be identified for any real use.Accordingly, special finite elements with singular stress fields should be considered.Elements of this kind have been developed via nodes shifting [19], singular topological transformation [20]- [29], embedding a known solution function of a singular problem [30]- [33], or using interpolation functions with singular derivatives [34]- [36].Progressive refinements have been brought to these techniques, in order to ensure compatibility of the crack-tip singular elements with the standard elements without using transition bridging elements.In this way it is possible to weaken the consistence requirements for the field functions, allowing to use certain discontinuous displacements and to develop more efficient ways to construct the stiffness matrix.
For a crack lying close to a bi-material interface like in laminated composites, the stress singularity is oscillating in nature, i.e. σ ij ≈r -1/2+ic while for a crack perpendicularly terminating at the interface the stress field is of the form σ ij ≈r λ-1 (0<λ<1).In the previous formulas r is the radial distance from the crack tip, λ is the lowest root of the characteristic equation, r ic is the oscillatory term due to the mismatch in elastic properties of the two bonded material.Therefore, an arbitrary order of singularity should be embedded in the elements.Elements with variable singularity power have been recently developed using different approaches and successfully applied for solving different problems in the papers by Zhou et al. [37] and Chen et al. [38].In [37] a three-dimensional variable-order boundary element with variable singularity order and different kinds of displacement formulations was suggested, while in [38] a special five-node singular triangular crack-tip element with extra basis functions in the G space and strain smoothing was proposed.
In the present paper, a solid singular, mixed wedge element for analysis of composite laminates is developed.It uses simple singular interpolation functions with variable singularity power that can degenerate into a non-singular standard representation as a particular case.This element is developed assuming the displacements and the out-of-plane stresses (i.e., the interlaminar stresses) as nodal d.o.f, since in this way, it a priori fulfills kinematic and stress contact conditions at the material interfaces, which are essential for keeping equilibrium, and it easily satisfy the stress boundary conditions in a point form.This element, which can be particularized as a displacement based element, is compatible with the standard ones (solid trilinear displacement based elements and non-singular mixed elements).The mixed formulation is chosen because it can provide more accurate results than displacement-based models with the same meshing and with a comparable processing time.The element is developed considering displacement and stress interpolating functions all of the same order, therefore the equilibrium is satisfied in an approximate integral form.The interpolation functions by Hughes and Akin [39] are used for the displacements, giving rise to singular derivatives with a variable singularity power that can be set as desired and can represent the regular (i.e.non-singular) behavior as a particular case.The interpolation functions for the nodal stresses are obtained deriving those for the displacements with respect to the coordinate that lies in a plane parallel to the layers where delamination/debonding cracks can rise and consequently a singularity line rises at the crack tip.Satisfaction of equilibrium in approximate integral form assuming a unique set indifferently for the displacement or the stress d.o.f.components is a simplification that does not result in an accuracy loss, as shown by Nakazawa [40].More accurate stresses than with displacement-based elements are obtained with the same meshing, along with an improved convergence rate.The present element is employed around the delamination crack singularity, while to ensure compatibility with no need of transition elements, the non-singular brick elements with the same d.o.f and standard three-linear interpolating functions developed by Icardi and Atzori [41] and extensively applied to analysis of composites are used elsewhere.Due to the interpolation scheme with quasi-standard features, the element can be easily implemented into the computer code [41] and used to account for the effects of singularities.The present solid modeling will result in an increase of data preparation and processing time with respect to plate models customarily adopted for the analysis of composite laminates, but simplifying assumptions will not be introduced across the thickness that could limit accuracy, as discussed, e.g., by Weißgraeber and Becker [42].So the present modeling can be more accurate with arbitrary geometries, lay-ups, loading and boundary conditions and the crack growth mechanism can be easily simulated by splitting the interfaces.The accuracy of the present modeling will be assessed considering sample test cases with singular stress fields taken form the literature and initially delaminated specimens, for which experimental results are available, whose failure loads will be predicted using fracture mechanics and a mesomechanic approach.In the former case, the analysis of delamination will be carried out using the singular representation, as in the papers by Ariza et al. [43] and Yao and Hu [44].In the latter case, the regular behavior will be set by removing the effects of the singularity, similarly to McQuigg [45].
The numerical results will show that the failure loads of initially delaminated specimens are accurately predicted with an affordable effort for a personal laptop computer either using the virtual crack closure technique, or a mesoscale model.Correct results are obtained with a relatively coarse meshing, thanks to the possibility to adapt the singularity power and owing to the improved accuracy of the mixed formulation.Of course, more accurate results are obtained when the singularity power is appropriately set and the regular behavior can be correctly predicted by setting the singularity power to zero.The displacement-based elements particularized from the mixed elements give less accurate predictions.

The wedge Element
Figure 1 represents the wedge element in its natural volume of coordinates (η, θ, ζ) and in physical coordinates (x, y, z).Here x, y and η, θ are assumed to lie in a plane Variable Singularity Power Wedge Element for Multilayered Composites parallel to that of the constituent layers, with x, η in the spanwise direction and y, θ in the transverse direction, as a consequence, z and ζ lie in the thickness direction.The displacements and the interlaminar stress components are indicated as u, v, w, σ xz , σ yz , σ zz , respectively.The axes x and y are required to lie in a plane parallel to the constituent layers, while z is required to lie in the thickness direction, so that the singular line from nodes 1 to 4 lies on an interface.Distorted elements with a reduced distance from the nodes 1 and 4 can be constructed in order to approximate a point singularity.The symbols N e disp , N e stress represent the interpolation functions for the displacements and the stresses, which are here indicated by D and S, respectively: It is desirable that the interpolation functions have a variable stress singularity power that can adapt to the problem, in order to treat all constituent materials, geometries, loading and boundary conditions.The interpolation functions should also be able to represent the regular behavior, in order to obtain a conventional element as a particular case.Here, the shape functions with singular derivatives developed by Hughes and Akin [39] are used for interpolating the displacements: The stresses are instead interpolated by differentiating the previous functions with respect to the coordinate ζ, which is supposed to lie in a plane parallel to the layers where delamination/debonding can rise: In this way, the stresses become singular at η =0, i.e. along the edge connecting nodes 1 and 4, with a singularity power μ≡ β -1 that can adapt to the problem under investigation.The regular behavior is obtained by setting β= 1 in (4) to (9) and assuming μ= 1 in (10) to (15).
As shown by Kaneko and Padilla [46], elements based on the shape functions defined above could have a convergence rate depending on the location of the interpolation points.Namely, locations can exist in which they could not work.However neither the present numerical applications, nor other authors have shown the occurrence of such prejudicial effects in the practical cases examined.
It could be noticed that equilibrium is fulfilled in an integral form, because (18) does not always produce intra-element self-equilibrating stresses.As mentioned in the introductory section, this representation is chosen because it simplifies the development of the element and does not have detrimental effects on accuracy.
As all the stresses should have the same singularity power, though they can have different singularity strength, the interpolation functions (10) to (15) are used for representing the inner variation also of the in-plane stress components.The element matrix and the vector of equivalent nodal forces are derived in a straightforward way by the Hellinger-Reissner mixed variational principle: by substituting the discretized representation of displacements and stresses given above.Please note that, the symbol b i represents the volume forces, u i represents the displacements of the points where these forces are applied, t i are the surface tractions, ε u ij =½(u i,j + u j,i ) represent the strains derived from the strain-displacement relations (the differentiation with respect to the physical coordinates is indicated by a comma), while σ kl are their counterparts obtained with the stress-strain relations: The stresses are represented as: the matrix [ D ] being obtained from the relation that defines the in-plane stress components by inserting zeros in the appropriate positions.Since the interlaminar stresses σ xz , σ yz , σ zz are assumed as nodal d.o.f., the interfacial stress contact conditions are a priori fulfilled when the nodes lie over a material interface.Though the continuity of the transverse normal stress gradient σ zz,z , which is also prescribed by the elasticity theory, and the boundary conditions σ zz,z = 0 at the upper and lower bounding faces cannot be enforced with the chosen d.o.f., they are spontaneously met, as observed by Icardi and Atzori [41] and by the present results.The present wedge element can be connected at nodes 2,3,5,6 to the non-singular hexahedral element by Icardi and Atzori [21], which has the same nodal d.o.f. and compatible tri-linear interpolation functions.
Singular and non-singular displacement-based counterpart elements can be particularized from the present wedge element, which is referred as W6-MX-SG, according to the scheme of Table 1.In the numerical applications, the results by the present element will be compared to those of its counterpart elements defined in Table 1, in order to show the advantages of the mixed formulation and the variable singularity power representation.
As customarily, a topological transformation is carried out from the physical volume (x, y, z) to the natural volume (η, θ, ζ) to uniform the computation of the volume integrals.This operation transforms any wedge element into the prismatic element of Figure 1, i.e. the nodal displacements e u i , e v i , e w i and the stresses e σ i xz , e σ i yz , e σ i zz into their counterparts e η u i , e θ v i , e ζ w i , e σ i ηζ , e σ i θζ , e σ i ζζ in the natural volume.To this purpose, the coordinates of the internal points are expressed in terms of the nodal coordinates through the standard tri-linear interpolation functions N j c obtained from the N j disp of ( 4) to ( 9) by setting β= 1: The derivatives with respect to the physical coordinates are computed as: 1 , , , , | ℑ | being the Jacobian of the transformation  In the numerical applications, the element matrix and the vector of the nodal forces will be derived using symbolic calculus, as in this way they are generated once at a time for all the elements.With this approach, the integrals involved in Eq. ( 16) will be computed in a closed form, thus avoiding the possible numerical instabilities due to Gaussian integration.
The readers are referred to Hoa and Feng [47] for a comprehensive discussion of stability and solvability of mixed elements.Summing up, for being admissible, these elements should have a number of displacement d.o.f. that is larger or equal to that of stress d.o.f., but certain choices of the interpolation functions could not yield meaningful results.An eigenvalue test performed over a single finite element with unit length edges, which is free of boundary conditions, should be carried out to ensure admissibility.Instead, a sufficient condition for solvability requires that the number of zero eigenvalues of the stiffness matrix should be equal to the number of rigid body modes, that the number of positive eigenvalues should be equal to the number of stress d.o.f., while the total number of zero and negative eigenvalues should correspond to the number of displacement d.o.f.An isotropic material with unit elastic moduli and a Poisson's ratio of 0.3, should be used in order to have eigenvalues with the same magnitude.The successful results of this test over the element W6-MX-SG and its particularizations (because the results coincide in the test case) is shown in Table 2, while those for the element B8-MX are shown in the paper by Icardi and Atzori [41].Hereafter, the virtual crack closure technique and a damage mesomodel used for predicting the delamination failure loads will be briefly reviewed.

Virtual Crack Closure Technique
Efficient ways for computing the energy release rates (ERR), which represents the drive force, since the crack propagates when they exceed a critical value, are based on the VCCT, as discussed in a comprehensive review by Krueger [48].It is assumed that a crack propagates from a node i, which is split into two nodes whose relative displacements in the three directions are indicated as Δ i sx, Δ i sy, Δ i sz.The nodal forces before propagation are indicated by F i x, F i y, F i z.The sole assumption of VCCT is that the energy required for a crack propagation Δα is equal to that required to close two separate crack surfaces of the same length Δα.The total ERR due to the crack propagation G is expressed as the sum of the ERR corresponding to the three modes G I , G II , G III , i.e. the opening (mode I), the forward shearing or sliding (mode II) and the anti-plane tearing or scissoring (mode III) modes: Σ being the area generated by the crack propagation Δα.
Usually, the so called two-step method is employed for computing G, because it can be successfully used with a not extremely refined meshing.This method uses the node forces before crack propagation and the node relative displacements after crack propagation: 1 ( ) As an alternative, the one-step method could be used, which however requires a refined meshing.It considers coupled the displacements of nearest nodes before crack propagation, instead of the relative displacements after crack propagation.
The meshing should be chosen by considering the balance between accuracy and computational efficiency and by doing also convergence rate considerations, because a too refined meshing can result disadvantageous by the viewpoint of convergence, as shown by Krueger [48].Under mode I, usually accuracy is not appreciably improved by refining the discretization, while it is improved of a small percentage under mixed modes (see Liu et al. [13]).Here the meshing size is chosen by stopping the refinement when precision is incremented less than 5 %.The following criteria, widely used for the analysis of delamination via VCCT, are employed for estimating the equivalent critical ERR of mixed-modes: B-K law (Benzeggagh and Kenane [14]) where G cr I , G cr II and G cr III represent the critical ERR for modes I, II and III.The choice of the exponents can determine quite different results, as observed by Liu et al. [13] and in the present numerical applications.The values 0.5, 1 and 2 have been used in literature, but the most accurate results are obtained with the latter two values.

Damage Mesoscale Model
This model considers an intermediate scale between that of the laminate and the micro-scale of the mechanisms of damage.Matrix microcracking, local delamination, diffuse damage and diffuse delamination are efficiently described by replacing the discretely damaged medium with an equivalent continuous homogeneous medium, which is equivalent to the damage micromodel from an energy standpoint and can simulate structures of industrial complexity.The behavior is described by means of two mesoconstituents, the single layer and the interface, which represents the thin layer of matrix between two adjacent plies.The damage state of each mesoconstituent is quantified by damage variables linked to the loss of stiffness of the mesoconstituent, whose evolution depends on damage forces that represent the variation of the energy with respect to the damage.These damage variables are assumed constant across the thickness of each mesoconstituent, but can vary changing the mesoconstituent.
The damage mesoscale model (DMM) by Ladevèze et al. [49] is employed in this paper.In summary, the displacement, strain and stress fields representing the exact micro-solution are represented as S M =S͂ + S̅ , S͂ being the solution of the problem with the damage removed and S̅ the solution of a residual problem where each cracked area is loaded by the residual stress σ̴ ͂ .Damage indicators are introduced, which are defined as the integrals of the strain energy for any basic residual problem.The strain energy E j p of the interface γ j is expressed as: where k͂ 1 , k͂ 2 and k͂ 3 are the elastic stiffness coefficients of the interface and I̅ 1 , I̅ 2 and I̅ 3 the damage indicators, which depend just on the material properties of the mesoconstituent layer.They can be conveniently computed by a 3-D finite element analysis, as shown by Ladevèze and Lubineau [50], since the damage state of a ply depends on the state of damage of the adjacent plies.To this purpose, in this paper a cell subjected to the residual stresses is discretized by non-singular wedge elements, assuming the stiffness of the interface equal to the shear modulus of the matrix divided by its thickness, which is assumed as the thickness of a ply divided by 5.The progressive failure analysis by DMM is carried out at each increment of displacement or loading by evaluating the homogenized stresses at any point and comparing them to the strength.Damage is extended to the points where the ultimate strength is reached.

Two-Material Wedge
As a first test, the sample case of a two-material wedge is considered (Figure 2).Many studies have been published since the 60's for this or similar cases, where analytical solutions have been presented.The exact elasticity solution by Hein and Erdogan [51] and the eigenvalue analysis by Yamada and Okumura [52] show that a singularity due to the dissimilar properties of the interfaced materials rise at the corner interface for certain side angles.The singularity power for this case considerably varies with geometry and material properties.With a mild variation of the properties and side angles of 60° or less there is no stress singularity, while for larger angles and dissimilar materials a strong singularity rises.The largest singularity power corresponds to the case with an interface angle of 90° when one of the materials is rigid.This latter case was simulated using wedge elements around the corner singularity and brick elements elsewhere, Variable Singularity Power Wedge Element for Multilayered Composites as this represent the most severe case.The problem was simulated subdividing the domain into 600 elements, 275 of which in the rigid sector and 325 in the elastic sector, with a progressive refinement along the interface and at the left corner (Figure 2), where the results are compared to the reference solution.Around this corner, 30 singular wedge elements were used.A Young's modulus of 7.3 GPa and a Poisson's ratio of 0.3 were assumed for the deformable sector, while a modulus larger by a factor 100 was assumed for simulating the rigid sector.The variation of the ratio of the elastic moduli of the two sectors produces the variation of the singularity power shown in Figure 3.A singularity power α=0.22 corresponds to the case with an interface angle of 90° and one of the materials rigid.This value was found by a log-linear procedure, because a straight line is shown in the stress vs. distance log-log plot at the right α for each case with different elastic ratios and side angles.The exponent μ that appears in the interpolation functions is related to the singularity power by: then the singular stress field is represented as: ) ( (30) ℜ being the radial distance from the singularity, κ ij the singularity strength and O( ℜ -α+1 ) negligible higher-order terms.Figure 4 shows the curves obtained in the log-log plots when α ≠ 0.22 and the straight line corresponding to α=0.22.Table 3 reports the results when the singularity power is incorrectly set to 0.5 (symbol *) or it is set to 0 and thus represents the non-singular case (symbol °).It is evident from these results that the interpolation functions with variable singularity power give more accurate predictions.Forward and in the next Section, other cases with different singularity powers are considered, in order to assess the accuracy of the present finite element simulation.Figure 5 shows the results with a crack at the tip.The variation of the stress singularity power with the ratio of the elastic moduli in the two sub-domains is shown for a Poisson's ratios of 0.2 of the elastic sector, this case being also considered by Hein and Erdogan [51] and Yamada and Okumura [52].To enable a comparison with the results reported in [52], the characteristic equation used to solve the eigenvalue problem was reconstructed in the finite element model simulation by work-conjugating the nodal stresses to the strains obtained from the nodal values of displacements and the interpolation functions, in a region around the corner interface.In this way, it was possible to compute the eigenvalue λ, which is related to the singularity power α since α=|1-λ|, (the symbol | | was used to indicate the modulus).The results of Figures 3 and 5 and Table 3 show that the solution by the present finite element model is free from spurious modes; it gives accurate predictions of the singular behavior and can correctly determine the real and imaginary parts of the eigenvalue λ.

L-Shaped Region
As a further test, it is now considered the L-shaped plane geometry represented in Figure 6, where a stress singularity is generated by the opening corner (point S).The reference solution for this case is represented by the finite element analysis with non-singular displacement-based triangular plate elements by Zienkiewicz and Taylor [53].In the present analysis, the problem was discretized considering a single layer of solid elements.Four singular wedge elements were used around the opening corner (point S), while 213 regular wedge elements were used elsewhere.The reference solution [53] shows that, approaching the corner singularity along the line B-C, σ xx becomes about eight times greater than far from it.A similar result is obtained by the present finite element model when the regular behavior is set in the wedge elements around the opening corner, non-singular elements being used in the reference solution.In Figure 6 also the results obtained using displacement based elements B8 are reported.The comparison of these results to those by the mixed elements W6-MX-SG and W6-MX shows that even with refining meshing, elements B8 always give less accurate results.This result confirms the well-known fact that hybrid and mixed elements are more accurate than displacement-based elements with the same meshing.Figure 7 shows that the results by the present finite element simulation are still in a very good agreement with the exact solution by Hein and Erdogan [51] and the eigenvalue analysis by Yamada and Okumura [52] for this case.Forward, applications will be presented to laminates.

Assessment of the Fracture Mechanics Model
Now a simply-supported (90°5/0°5/90°5) cross-ply beam under a point centered loading is considered to the purpose of assessing whether the implementation of the fracture mechanics model was correct.The beam has a transverse crack in the matrix and an interfacial delamination at the upper interface of the intermediate layer.The beam is 2.187 mm thick and has a length of 50.8 mm, the constituent materials have the properties reported in Table 4 and the loading has a magnitude P=13.345N.
For this case, the results by Liu et al. [54] and Bui et al. [55] are available for comparisons.The width was chosen large enough to avoid the interaction among the stresses from the free edge and those due to the delamination crack.
In the present finite element simulation, a width of 25.4 mm was assumed.Owing to the symmetry of the problem, only a quarter of the beam was discretized.Span, width and thickness are subdivided into 9, 12 and 3 rows of elements, respectively.As usual, the wedge elements are used at the crack tip, while the brick elements are used elsewhere.In this case, 15 wedge elements were used at each of the tips of the interfacial delamination crack.The results by Liu et al. [54] and Bui et al. [55] reported in Figure 8 show the variation of the strain energy release rates G I and G II separately and their sum G I +G II as a function of the ratio a/h, a being the half length of the delamination crack and h the thickness of the beam.The results by the present finite element model reported in this figure were obtained either by setting the singularity power to 0.5 to enable a comparison with [54] and [55], or setting it with the log-linear procedure.Either displacement-based elements that can be particularized by the mixed elements, as shown in Table 1, or the wedge and brick mixed elements were considered in the analysis.The results of Figure 8 show that the strain energy release rates are correctly computed by the present finite element simulation.However, also in this case the most accurate predictions are still obtained by the mixed elements and when the singularity power is set by the log-linear procedure.Once shown the finite element model to be accurate and the fracture mechanics model correctly implemented, delaminated specimens will be analyzed.These new tests will also serve as assessments of the damage mesoscale model.

Initially Delaminated Specimens
The following sample cases are considered, whose experimental failure loads are known: Case A. Double cantilever unidirectional beam made of AS4-3K/PEI layers under mode I; Case B. Double cantilever angle-ply beam made of XAS/931C layers under mode I; Case C. Angle-ply end-loaded split beam made of T400/6376C layers under mode II; Case D. End-notch-flexure cross-ply beam made of AS4-3501-6 layers under mode II.Case A. This case refers to a beam made of 16 unidirectional layers aligned in the 0° direction, each 0.22 mm thick, whose properties are given in Table 5.The beam is 146.5 mm long; it has a width of 20.0 mm and a thickness of 3.52 mm.A centered non-stick film 59 mm long and 15 μm thick is interposed at the symmetry plane interface, in order to give rise to an initial delamination.G Icrit was determined by Fassine and Pavan [56] as 1.841 KJ/m 2 and the experimental failure load as 74.94 N. The delamination load for this case is reported in Table 6.The computations were carried out discretizing the beam into 30, 10 and 16 solid elements in the spanwise, width and thickness directions, respectively, with a progressive refinement at the tip where 15 wedge elements were used.
The results show that a more accurate prediction is obtained by VCCT when the proper singularity power is set, than with setting the non-singular behavior or a singularity power of 0.5.Since the failure occurs as a matrix failure under traction, the contributions by G II and G III are negligible, thus the propagation criteria ( 25) and ( 27) are reduced to G cr equiv = G cr I , while (26) becomes G equiv /G cr = (G I /G cr I ), consequently just α m =1 has meaning.The damage mesoscale model (DMM) by Ladevèze et al. [49] also provides results in a good agreement with the experiment carrying out the homogenization process by the finite element model using as strength the matrix strength under traction and setting the non-singular behavior in the wedge elements.
In Table 6, the processing time required to solve the problem by a laptop computer with a dual-core, 64 BIT operating system, 2.20 GHz CPU and a 4 GB RAM is reported.Please notice that the times provided comprise solution of the finite element scheme and computation of the failure load, but not the generation of the element matrices and the homogenization process, which are carried apart and take less than 10 minutes.VCCT appears advantageous over DMM as it does not require additional costs due to the homogenization process, while the singular representation offers the capability to predict presence and nature of singularities at the same costs of the non-singular representation.However, in this case an accurate prediction of the failure load can be obtained also using stress-based criteria (the results by Chang-Springer's [57] and Davila-Camanho's [58] criteria are reported) if the stresses σ͂ ij computed by setting the non-singular behavior are averaged over a distance from the crack tip equal to the thickness of the laminate.
Case B. The beam is made of XAS/931C layers, each 0.125 mm thick, with the material properties reported in Table 5 and a [-45°/0° (45°) 2 /0°/(-45°) 2 /0°/45°] 2 lay-up.The beam has a length of 150 mm, a width of 30 mm and a thickness of 3 mm.A Teflon film with a length a 0 =29 mm and a thickness of 10 μm is interposed at the mid-plane to generate an initial delamination.Due to the 3-D stress field, delamination in this case occurs as a combination of the three modes.Robinson and Song [59] found an experimental failure load of 87 N for this case and numerically estimated G crit = 0.35 KJ/m 2 .The results by the present finite element model are still reported in Table 6.The beam was discretized by 30, 10 and 21 solid elements in the spanwise, width and thickness directions, respectively.The meshing was refined at the crack tip where the wedge elements were used, while the brick elements were used in the other regions.Also in this case 15 wedge elements were used.The mounting bulks bonded to the edges in order to apply the opening loading were discretized, as they represent a local increment of stiffness.The critical ERR for the three single modes were numerically estimated as suggested by Robinson and Song [59].
The results show that the propagation criteria ( 25)-( 27) determine a different failure load with a different choice of the exponents (the law (26) was applied using equal exponents), as well as the important role played by the singularity power assumed in the interpolation functions.As shown by comparison with experiments, the results obtained with the singularity power set to 0 or 0.5 are again less accurate than those obtained setting the proper value with the log-linear procedure.
The results by the DMM model have the same accuracy than those by VCCT also in this case.Since the stress state is no longer generated by a pure transverse traction, the Davila and Camanho's criterion was used with DMM as matrix failure rule.Table 6.Failure loads of initially delaminated specimens [N] and computational time {}; λ1 variable singularity power; λ2 singularity power=0.5;λ3 singularity power=0.DMM: Damage mesoscale model [49]; C-S: Chang and Springer's criterion [57]; D-C: Davila and Camanho's criterion [ The processing time is reported in Table 6 also for this case leaving apart the generation of the element matrices and the homogenization process.
In this case, the Chang and Springer's and Davila and Camanho's criteria predict failure loads that are comparable each other and to the experimental load only if the stresses σ͂ ij are averaged over a distance from the delamination crack tip at least of 4 mm, which is larger than in the previous case.Please notice that numerical tests have shown that these criteria give monotonically increasing failure loads increasing the distance from the crack tip at which the stresses are evaluated, so they can overestimate the failure load at a sufficiently large distance if stress averaging is uncorrected.
Case C. The case studied by Choi et al. [60] of a [(-45°/ 0°/ 45) 2s / (45°/ 0°/ -45) 2s ] angle-ply end-loaded split beam made of T400/6376C layers under mode II is considered.The beam is 100 mm long; it has a width of 24.5 mm and a thickness of 3.4 mm.Each constituent layer, whose material properties are reported in Table 5, is 0.142 mm thick.An initially delaminated area with a length of 60 mm is generated by interposing a Teflon film with a thickness of 12.5 μm at the mid-plane interface.Owing to the primary importance of interlaminar stresses, each ply was subdivided into four layers of elements, i.e.48 elements were used across the thickness, while 30 elements were used in the spanwise direction, which were refined at the tip, and 5 were used in the width direction.Also in this case the mounting bulks were discretized and the critical ERR were still numerically estimated.As usual, wedge elements were used at the crack tip, whereas brick elements were used in the other regions.
The results show that the choice of the exponents of the propagation criteria produce larger variations than in the previous sample case, since the contributions by modes II and III are much more important.As shown by Table 6, the best results by VCCT were obtained by using the exponent 1 with all the propagation criteria, but good results were also obtained by using the exponent 2, while the exponent 0.5 gave inaccurate predictions.Again, the most accurate results were obtained by setting the singularity power with the log-linear procedure.
DMM still predicted a failure load in good agreement with VCCT and with the experiment, computing again the matrix failure through the rule by Davila and Camanho.The stress-based criteria still can provide accurate results considering the stresses computed using non-singular elements if these stresses are correctly averaged.The Davila and Camanho's and Chang and Springer's criteria predicted failure loads in a well agreement with the experiment when the stresses σ͂ ij were averaged over a distance of about 6 mm from the crack tip, while the Choi-Chang's delamination criterion used for a comparison required just a distance equivalent to the thickness of a ply.
The computational time again shows that the present analysis is affordable by a personal computer.However, it is presumed that a refined meshing should be used to more accurately determine the participation of the crack tip stresses for this case.
Case D. The latter case examined is that of an end-notch-flexure cross-ply beam under mode II, with a [(0°1 0 /90°/0°1 0 )/0°/90°/(0°1 0 /90°/0°1 0 )] lay-up, a length of 38.1 mm, overall thickness of 5.1 mm and a width of 16.5 mm.Each constituent layer is 0.116 mm thick and has the properties reported in Table 5 (material AS4-3501-6).A Teflon film is interposed at the mid-plane during manufacturing, in order to simulate an initially delaminated area of length 19.05 mm.Kim and Sun [61] found an experimental failure load of 1192 N and the critical strain energy release G II crit =0.515 KJ/m 2 .Each constituent layer was discretized by a layer of elements, i.e. 44 elements were used across the thickness and the width was discretized by 5 elements, while 30 elements with a progressive refinement at the crack tip were used in the span direction.Table 6 shows that the results by VCCT are now very sensitive to the choice of the exponents of the propagation criteria.Still better results are obtained by setting the proper singularity power of wedge elements at the crack tip with the log-linear procedure.Also in this case accurate results are computed by the DMM model carrying out the homogenization process by finite elements and using the Davila and Camanho's criterion.Application of the Davila and Camanho's and Chang and Springer's criteria using the stresses σ͂ ij provided by the finite element analysis with the regular behavior shows that accurate results are obtained averaging the stresses over a distance of 5 mm from the crack tip.It is shown that also in this case the computational time is affordable by a personal computer.

Concluding Remarks
A new mixed singular wedge element with a variable singularity power that can adapt to the problem and can represent the non-singular behavior as a particular case was developed for analysis of damaged composites.The accuracy of the element was initially assessed considering sample cases with singular stress fields due to dissimilar properties of materials, geometry or cracks, for which analytical solutions are available for comparisons.The numerical results proved that the present finite element simulation with mixed singular wedge elements around singularities and solid mixed elements elsewhere accurately predicts the presence and the nature of singularities, the related stress fields, the participation of the singular stresses and the regular behavior with an affordable computational effort for a laptop computer.The variable singular representation and the mixed formulation gave better results than their displacement-based, non-singular, or singular counterparts with a fixed singularity, confirming the results in literature.The failure loads of initially delaminated beams were computed using the virtual crack closure technique, a damage mesoscale model and comparisons were also made using stress-based criteria.The virtual crack closure technique was carried out setting the singular behavior with the appropriate singularity power in the wedge elements, while for the damage mesoscale model and the stress-based criteria the non-singular behavior was set.The comparison with experiments proved that accurate predictions can be obtained using fracture mechanics, but nearly equally accurate results can be obtained using damage mesomechanics and even stress-based criteria with an appropriate averaging of stress, as claimed in literature.The homogenization process required by the damage mesomechanics model, resulted in an additional cost for the damage mesoscale model.The results of the fracture mechanics model were sensitive to the exponents employed for estimating the critical energy release rates by the propagation criteria.The stress-based criteria required the stresses to be averaged over an appropriate distance from the tip, which appeared case dependent The present finite element simulation is accurate and efficient with a relatively coarse meshing and low computational costs.In all the examined cases, the computing time and the memory occupation were affordable by a personal computer, thus the present element can be successfully employed with the three classes of models currently used for the delamination analysis.

Figure 1 .
Figure 1.The wedge element in its natural volume and numbering of nodes Their corresponding nodal values, which represent the nodal d.o.f. of the present element, are indicated as: {D e | S e } )Power law(Wu, and Reuter [15]) Universal Journal of Engineering Science 2(1):16-

Figure 2 .
Figure 2. Geometry, loading and meshing for the two-material wedge problem

Figure 3 .Figure 4 .
Figure 3. Variation of the stress singularity power α with the ratio of the

Figure 5 .
Figure 5. Two material wedge problem: real and imaginary parts of the eigenvalue λ with a crack

Figure 6 .
Figure 6.L-shaped region: variation of the in-plane stress σxx with the distance from the opening corner S

Figure 7 .
Figure 7. Singularity eigenvalue analysis for the L-shaped region

Table 4 .Figure 8 .
Figure 8. Variation of the energy release rates GI , GII and GI +GII with the ratio a/h of the delamination crack half-length with respect to the thickness 12 13 14 15 16 ;

Table 1 .
[41]ents that can be particularized from the present element and from[41]

Table 2 .
Solvability test for the present singular wedge element

Table 3 .
Singularity power of the two material wedge problem by setting a different singularity power in the interpolation functions (+)eigenvalue analysis

Table 5 .
Material properties of initially delaminated specimens