July 2013 

NASA/CR–2013-218022 

Floating Node Method and Virtual Crack 

Closure Technique for Modeling Matrix 

Cracking-Delamination Interaction 

N. V. De Carvalho 

National Institute of Aerospace, Hampton, Virginia 

B. Y. Chen, S. T. Pinho, and P. M. Baiz 

Imperial College London, London, United Kingdom 

J. G. Ratcliffe 

National Institute of Aerospace, Hampton, Virginia 

T. E. Tay 

National University of Singapore, Singapore 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

Access the NASA STI program home page 

at http://www.sti.nasa.gov

E-mail your question to help@sti.nasa.gov

Fax your question to the NASA STI 

Information Desk at 443-757-5803 

Phone the NASA STI Information Desk at 

443-757-5802 

Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320 
National Aeronautics and 

Space Administration 

Langley Research Center Prepared for Langley Research Center 

Hampton, Virginia 23681-2199 under Cooperative Agreement NNL09AA00A 

 

July 2013 

NASA/CR–2013-218022 

Floating Node Method and Virtual Crack 

Closure Techique for Modeling Matrix 

Cracking-Delamination Interaction 

N. V. De Carvalho 

National Institute of Aerospace, Hampton, Virginia 

B. Y. Chen, S. T. Pinho, and P. M. Baiz 

Imperial College London, London, United Kingdom 

J. G. Ratcliffe 

National Institute of Aerospace, Hampton, Virginia 

T. E. Tay 

National University of Singapore, Singapore 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

Trade names and trademarks are used in this report for identification only. Their usage does not constitute 

an official endorsement, either expressed or implied, by the National Aeronautics and Space 

Administration. 
1 

Abstract 

A novel approach is proposed for high-fidelity modeling of progressive damage and failure in 

composite materials that combines the Floating Node Method (FNM) and the Virtual Crack 

Closure Technique (VCCT) to represent multiple interacting failure mechanisms in a meshindependent fashion. In this study, the approach is applied to the modeling of delamination 

migration in cross-ply tape laminates. Delamination, matrix cracking, and migration are all 

modeled using fracture mechanics based failure and migration criteria. The methodology 

proposed shows very good qualitative and quantitative agreement with experiments. 

1 Introduction 

Damage in composite materials generally occurs as a combination of different and interacting 

failure mechanisms, e.g. delamination and matrix cracking. Capturing these interactions 

accurately is essential to confidently model and predict progressive damage and failure. Many 

approaches have recently been proposed that explicitly model different failure mechanisms and 

attempt to capture their interaction [1-3]. However, no single strategy has yet surfaced as totally 

satisfactory. The present approach combines the Floating Node Method (FNM) [4] with the 

Virtual Crack Closure Technique (VCCT) [5,6] to explicitly account for different failure 

mechanisms and their interaction. Delamination, matrix cracking, and migration events are all 

modeled with the same FNM element using fracture-mechanics-based failure and migration 

criteria. The approach is validated using recent experimental results in which a setup capable of 

isolating a single complete migration event was developed [7]. Migration is defined as the 

transition of a delamination from its current ply interface to a new interface via transverse ply 

cracking. The results obtained offer a simple, yet very challenging, benchmark for the validation 

of numerical and analytical approaches aimed at high-fidelity modeling of composite materials, 

such as the one proposed in the present paper. 

2 Floating Node Method (FNM) 

2.1 Introduction 

The FNM is a recently proposed numerical method, capable of representing multiple evolving 

discontinuities in solids [4]. Compared to other methods used for the same purpose, such as 

eXtended Finite Element Method (XFEM) or the Phantom Node Method (PNM), it provides the 

following key advantages. First, FNM does not introduce an error on the crack geometry when 

mapping to natural coordinates. Second, it does not require numerical integration over only part 

of a domain. Third, it is ideally suited for the representation and incorporation of multiple and 

complex networks of weak, strong and cohesive, discontinuities. Fourth, FNM yields the same 

solution as a finite element mesh where the discontinuity is represented explicitly, and finally it 

is conceptually simpler than the PNM or XFEM. 
2 

2.2 Element formulation 

The static equilibrium of a body with volume under body forces with density (acting on ) 

and traction acting on the boundary can be expressed in the weak form as: 





where is the displacement vector; is the test function; is the strain tensor related to 

through the differential operator relative to Cartesian coordinates as ; and is the 

stress tensor related to the strains through Hooke’s law as , with being the constitutive 

tensor. In the Floating Node Method, each real node “” is characterized by its nodal coordinates 

and associated Degrees of Freedom (DoF) . In addition, an FNM element contains a suitable 

number of floating DoF without pre-defined associated nodal position vectors. These additional 

floating nodes are used to represent discontinuities. Their number varies with the number and 

type of discontinuities, weak or strong, modeled within each element. A strong discontinuity is 

defined as a jump in a field quantity (e.g. displacements), while a weak discontinuity is defined 

as a jump in the gradient of a field quantity (e.g. strains). Figure 1 shows an example of such an 

element with four nodes and four additional sets of floating DoF required to represent a strong 

discontinuity.

2.2.1 Element formulation without weak/strong discontinuity 

If there is no discontinuity to be modeled by the element (e.g. before failure), the formulation is 

the same as in the standard Finite Element Method (FEM). The vector of nodal coordinates is 

given as . In the case of Figure 1, x is given by: 





The Jacobian of the transformation from physical coordinates x to natural coordinates is: 







where is a matrix of shape functions. Assuming an isoparametric formulation, the 

displacement field u is related to the real DoF q through: 

	
3 

In the case of Figure 1, . The strains can also be expressed in terms of the 

real DoF as: 




where , leading to the stiffness matrix: 





 

where is the integration domain (in natural coordinates) corresponding to (in physical 

coordinates), and to the vector of nodal forces: 





 



 


where is similarly the boundary corresponding to . The weak form of equilibrium, Equation 

1, becomes: 

 

2.2.2 Element formulation with weak/strong discontinuity 

Once a discontinuity in the element is defined, the element is split in two or more partitions 

(depending on the discontinuity). Without loss of generality, a case in which the element is split 

in two partitions, and , is illustrated in Figure 1. For each partition, and , a vector of 

nodal coordinates, and , is defined. For the case in Figure 1: 









Each partition has a separate Jacobian: 









The displacements and , in partitions and , respectively, are interpolated separately 

from the DoFs and associated with each partition: 



These DoFs, and , did not have an associated position beforehand, and therefore were 

considered to be “floating”. As the discontinuity is defined, they are linked to the respective 
4 

position vector. In the case of Figure 1, which represents a strong discontinuity,



and . Note that there are four different sets of 

floating DoF: is different from and is different from . If a weak discontinuity were 

to be modeled, only two sets of floating DoF would be included in the element, and the DoF with 

a prime would coincide with those without a prime. 

The strains then become: 



	






The stiffness matrices for partitions and are: 



 







 






and the force vectors: 

 

 

 

 

Finally, the equations of equilibrium can be written as: 



2.2.3 Element topology and assembly 

To illustrate this approach, a floating node element capable of modeling two weak interfaces 

with an arbitrary crack between them was implemented, as shown in Figure 2(a). The nodal 

position of the initial integration domain, and the real and floating DoF are indicated in the 

figure. Each floating DoF is associated with either an edge or the inner domain of the element. 

Once weak/strong discontinuities are detected, the floating DoF at the edges are used to 

determine the solution for the nodal positions given by the intersection of these discontinuities 

with the element edges (Figs. 2(b) to 2(e)). The additional inner floating DoF are used to 

determine the solution for the nodal positions originated by the intersection of multiple cracks 

within the element (Fig. 2(f)). In Figures 2(b) to (f) only the floating DoF used are represented. 

In each case, the floating DoF that are not used have no assigned position, only a topological 

relation to an edge or to the inner domain. During assembly, inner floating DoF can be removed 

from the analysis through static condensation. The floating DoF associated with the edges are 

assembled with the corresponding DoF of neighboring elements, and will therefore have a 
5 

unique position in the global DoF vector. All floating DoF are then used as required to model 

weak/strong discontinuities as they evolve throughout the simulation. 

The approach can be extended to model an arbitrary number of weak and strong discontinuities 

within an element, and therefore an arbitrary number of plies and interacting cracks, provided the 

topology is adequately defined, and a sufficient number of floating DoF are used. 

3 Virtual Crack Closure Technique (VCCT) 

In the present work, the FNM method has been coupled to the Virtual Crack Closure Technique 

(VCCT) [5]. An extensive review of VCCT is provided in [6]. Traditionally, VCCT is used to 

obtain energy release rates in cases where the crack path is known beforehand and can be aligned 

with the elements’ edges, such as the case for delamination. 

In VCCT, the Mode I and II energy release rates are obtained from [6]: 





























where and are the normal and tangential components of internal force vector at the crack 

tip; and are the normal and tangential components of the displacement jump 

between the nodes immediately behind the crack tip; and and are the lengths of the crack 

in the elements behind and ahead of the crack tip, respectively (Fig. 3). When combining VCCT 

with FNM, the forces and displacements needed to obtain the energy release rates, Eqs (21) and 

(22), are computed at the floating DoF as the crack develops (Fig. 3). Additionally, floating DoF 

are also added along a virtual crack plane up to a distance , named the ‘enrichment radius’, . 

The associated floating nodes provide additional degrees of freedom along this path, which 

enable the accurate modeling of the deformations near the crack tip [4]. At the end of the 

enrichment radius, the displacements of the floating DoF are interpolated from the displacements 

of the real DoF. In the present work, the enrichment radius was chosen to be equal to the largest 

in-plane dimensions of the numerical model. 

4. Delamination propagation 

In this paper, delamination is modeled using VCCT to compute the energy release rates and 

at each delamination front position. These are then used in a failure criterion: 







where , and is the critical energy release rate given by [8]: 

  







	
6 

For delamination growth according to this criterion, is assumed to equal the critical energy 

release rate of the interface, 

. The delamination front is advanced by one element along the 

projected crack path when 

.

5 Matrix crack propagation 

In the present study, matrix crack propagation is also modeled with VCCT. However, unlike 

delamination, the matrix crack is assumed to propagate following a Mode I path in the throughthickness direction, as supported by experimental evidence [7]. Therefore, the failure criterion is 

assumed to be given by: 








where the critical energy release rate is given by  . 

At each crack growth increment, the Mode I crack path orientation is approximately determined 

using a maximum tangential stress criterion. This criterion, typically written in terms of the 

Mode I and II stress intensity factors [9], can also be written in terms of energy release rates 

and , where the angle that maximizes the tangential stress is obtained by evaluating: 



















and for the two solutions obtained choosing the one that maximizes given by: 













  

 


where is the tangential component of the internal force vector at the crack tip. The tangential 

stress relates to by: 











where, 



 



 



is the Young’s modulus; and the Poisson’s ratio. Once the angle is determined, the criterion 

given by Eq. (25) is assessed. If the criterion is met, the crack advances by one element along the 

projected crack path. 
7 

6 Migration 

6.1 Delamination to matrix crack – migration criterion 

Extensive research has been done on migration of cracks from a bi-material interface. He and 

Hutchinson [10] have proposed a necessary condition for migration to occur: 















 

This condition states that for a crack to migrate from an interface into a given Material B, as 

shown in Figure 4, the ratio between the energy release rate associated with its growth into that 

material  , (at an angle  that maximizes its value) and the critical energy release rate of 

that material 

, has to be higher than the ratio between the energy release rate associated with 

the growth of the crack along the interface 

, and the critical energy release rate 

 of that 

interface. 

With the approach proposed, evaluating Eq. (31) would require more than one FE analysis at 

each interface crack position. In each of these analyses, a branched crack would need to be 

initiated at different tentative angles in order to determine the angle, , that would 

maximize . 

In the present work, an alternative criterion is proposed. With this criterion, all variables needed 

to determine whether migration occurs, can be directly obtained from the FNM results at the 

each interface crack position. The criterion states that: (i) the sign of the tangential component of 

the internal force vector , associated with Mode II shear displacement, dictates the direction 

and therefore the material into which the interface crack migrates and (ii) migration can only be 

completed if it is energetically favorable. The criterion is given by 

as: 



















 

  















 

  



If 

the interface crack does not migrate, and if or the crack migrates 

into material A or B, respectively. In the present form, Eq. (32) assumes that internal force vector 

is computed at the lower surface (Material B). If the internal force vector is computed at the 

upper surface the signs of the inequalities in Eq. (32) need to be reversed. Depending on the sign 

of , the ratio between the energy release rate determined at the interface, and the critical energy 

release rate of Material A or B, 





or 







, is compared to the ratio between the energy release 

rate determined at the interface and the critical energy release rate of the interface, 





. If the 
8 

ratio 





or 







is greater than 





, the necessary condition for migration is met, and 



or 

. With this criterion, all variables needed to determine whether migration occurs, 

and 

, can be directly obtained from the FNM results. Note that although this criterion is 

necessary, it is not sufficient for migration to occur unless the ratio being assessed is greater than 

one. 

In the present work, the interfacial crack, or delamination, propagates between a  and a 

ply. Assuming Material A is the  ply, and Material B is the  ply, the critical energy release 

rate to propagate a crack into a  ply is much greater than the critical energy release rate of the 

interface, since fiber fracture is required for that process to occur, i.e. 

 [11]. 

Thus, 





is always smaller than 





. Therefore, even if the sign of favors migration into 

Material A, migration will not occur, since it is not energetically favorable. Typically for 

composite laminates, the critical energy release rate of an interface crack (delamination) and of a 

matrix crack are of the same order of magnitude, i.e. 

 [12]. Hence,  is assumed 

to be given by Eq. (24) and 



 to be equal to  . Therefore, assuming Material B is the 

ply, provided the sign of favors migration into Material B and it is energetically favorable, 

migration can occur. 

6.2 Matrix crack to delamination – migration criterion 

When a matrix crack reaches a weak interface, it is assumed to trigger delamination. This 

delamination is at first contained within the FNM element, as illustrated in Figure 2(f). In the 

following steps, the delamination will propagate, or migrate, as detailed in Sections 4 and 6.1, 

respectively. 

7 Simulation of delamination migration specimens 

7.1 Delamination migration test configuration 

In reference [7] a test aimed at investigating delamination migration was proposed. Its 

configuration is illustrated in Figure 5. The specimen possesses three main features that permit 

the controlled observation of delamination growth followed by migration to another ply 

interface. First, the specimen geometry is in the form of a beam with the intent of promoting 

uniform delamination growth and migration across the specimen width. Second, the specimen 

contains a PTFE insert (acting as an artificial delamination) at an interface between a  ply 

(specimen span direction) and a stack of four  plies (specimen width direction) (Fig. 5). This 

provides an opportunity for the delamination to migrate to another ply interface by kinking 

through the  ply stack. Third, the specimen can be loaded in a manner to cause delamination 

growth from the PTFE insert that eventually migrates to another ply interface. This sequence of 

fracture events is made possible by the way in which specimen loading affects shear stresses 

acting across the delamination front. 
9 

7.2 Finite Element Model 

The experiments presented in [7] were simulated using the finite element solver 

ABAQUS/Standard® (Implicit). Plane-stress and a small-displacement formulation were used. 

The floating node element (Fig. 2) was developed and implemented as a user defined element 

(UEL) in ABAQUS, and applied in the center region of the model, as shown in Figure 6. The 

specimen layup used in the experiments and in the model, is also given in Figure 6. In the layup 

description T represents the PTFE insert. Each block of plies of the same orientation was 

modeled with a separate element CPS4 (through-thickness), except for the highlighted plies at 

the center of the model (Fig. 6). These were modeled within a single FNM element, as illustrated 

in Figure 2. The material properties used are given in Tables 1 and 2. A prescribed displacement 

was applied at different positions, , at a single node, to simulate the displacement-controlled 

test. Four load positions were simulated corresponding to (Fig. 6). 

In the experiments, the specimens were subject to moisture intake from the ambient air leading to 

relaxation of residual thermal stress. Therefore, residual thermal stresses were not considered. 

7.3 Load-displacement response and sequence of events 

The computed load-displacement curves obtained for different values of are given in Figure 7. 

Even accounting for the compliance of the test fixture, measured in [7], all load-displacement 

curves obtained from the simulations show a slightly stiffer response than what was observed 

experimentally. The maximum load is in general very well predicted for the cases to 

, and slightly over-predicted for the case (Figs. 7(a) to 7(d)). For , 

the simulations predict unstable delamination after the peak load (Fig. 7(a)). This unstable 

delamination stops just before migration. After the migration event, delamination continues 

stably along the next  interface. A similar sequence of events was observed 

experimentally [7]. However, in the experiments, after the first unstable event, further loading 

was necessary before migration occurred. For the cases to , the simulations 

predict that (i) the first peak in load is followed by a region of stable delamination growth (Figs 

7(b) to 7(d)) that increases with the load offset, and that (ii) stable delamination is followed by a 

region of unstable delamination, which corresponds to the load-drop in Figures 7(b) to 7(d). The 

migration event occurs within this load-drop. After the migration, unstable delamination 

continues propagating along the next  interface. Loading the specimen further eventually 

leads to a transition from unstable to stable delamination propagation. This sequence of events is 

very similar to what was reported experimentally [7]. The main difference is that in the 

experiments, for the cases and , after the peak load, delamination 

propagates through a series of unstable events, rather than stably, as predicted by the simulations. 

7.4 Migration location 

Figure 8(a) provides an illustration of the typical crack paths predicted by the simulations and 

obtained experimentally. Figure 8(b) compares numerical and experimental results for the 

distance from the load-application point at which migration onset occurs, , as a function of the 

normalized load offset, . Overall, except for , migration is predicted to occur closer 

to the load application point than what was observed experimentally. The simulations also 
10 

predict a decreasing trend in with , which is similar to what was observed 

experimentally for the cases to . 

7.5 Matrix crack path 

In Figure 9, the simulated matrix crack paths obtained for the different load positions was 

superimposed on an image of a typical crack path observed experimentally, where . The 

predicted matrix crack path agrees well the crack path observed experimentally [7]. The 

simulations do not predict variation in the matrix crack path with a change in load offset. 

7.6. Discussion 

All load-displacement curves obtained numerically show a slightly stiffer response than the 

experimental results. In the simulations, the engineering properties used were obtained from 

uniaxial test data. These are known to be higher than their flexural counterparts [13]. Since the 

specimen is predominantly loaded in flexure, the difference between uniaxial and flexural 

properties is therefore likely to be the cause for the overall stiffer response obtained in the 

simulations. Overall, the simulated load-displacement curves, sequence of events, and failure 

morphology are in very good agreement with the experimental results. For the case , 

experimental results show an increase in load at the end of the first unstable event, before 

migration. This increase in load is likely to be caused by the development of fiber bridging 

during the unstable part of the delamination, which was not accounted for in the model. For the 

load cases to , the only significant difference between simulations and 

experiments is that delamination growth onset and propagation, prior to the final load-drop see 

(Figs. 7(c) and 7(d)), does not occur stably as predicted in the simulations, but rather through a 

series of unstable events. In these experiments the unstable delamination growth onset was likely 

due to the presence of resin rich pockets at the tip of the PTFE insert, while the subsequent 

unstable events, prior to the final load-drop, may be due to features accompanying delamination 

growth, such as fiber bridging. Both of these were not accounted for in the simulations. 

The migration location was slightly under-predicted for all cases except . Nevertheless, 

the migration location trend observed experimentally, for the cases to , was 

well captured. This further supports the validity of the approach, and in particular, of the 

migration criterion proposed. The latter is additionally confirmed by the good agreement 

between the predicted matrix crack path and that observed experimentally. This also suggests 

that the criterion used to determine the matrix crack orientation (as it propagates in the throughthickness direction) is adequate. 

8. Conclusions 

A novel approach that combines the FNM (Floating Node Method) with VCCT (Virtual Crack 

Closure Technique) has been demonstrated to adequately capture delamination migration in 

cross-ply laminates. Delamination, matrix cracking, and migration, are all modeled with the 

same FNM element using fracture mechanics based failure and migration criteria. Results show 

very good agreement with the maximum load, and migration location obtained experimentally. 

In addition, simulations capture the overall failure morphology and typical crack path, observed 

in the experiments. 
11 

References 

[1] D. Ling, Q. Yang, B. Cox, “An augmented finite element method for modeling arbitrary 

discontinuities in composite materials,” International journal of fracture, Vol. 156, No. 1, 

pp. 53-73, 2009. 

[2] E. V. Iarve, M. R. Gurvich, D. H. Mollenhauer, C. A. Rose and C. G. Dávila, “Meshindependent matrix cracking and delamination modeling in laminated composites,” 

International Journal for Numerical Methods in Engineering, Vol. 88, No. 8, pp. 749-773, 

2011. 

[3] F. P. van der Meer, L. J. Sluys, S. R. Hallet and M. R. Wisnom, “Computational modeling of 

complex failure mechanisms in laminates,” Journal of Composite Materials, Vol. 46, No. 5, 

pp. 603-623, 2012. 

[4] B.Y. Chen, S. T. Pinho, P. M. Baiz, N. V. De Carvalho and T. E. Tay, “Floating node method 

for the modeling of multiple discontinuities within an element,” 19th International conference 

on Composite Materials, Montreal, 2013. 

[5] E. F. Rybicki and M. F. Kanninen, “A finite element calculation of stress intensity factors by 

a modified crack closure integral,” Eng. Fract. Mech. Vol. 9, pp. 931-938, 1977. 

[6] R. Krueger, “Virtual crack closure technique: history, approach and applications,” Applied 

Mechanics Reviews, Vol. 57, pp. 109-143, 2004. 

[7] J. G. Ratcliffe, M. W. Czabaj, T. K. O’Brien, “Characterizing delamination migration in 

carbon/epoxy tape laminates,” Proceedings of the American Society for Composites – 

Twenty-Seventh Technical Conference, Arlington, TX, USA, paper 313, 2012. 

[8] M.L. Benzeggagh and M. Kenane, “Measuring of mixed-mode delamination fracture 

toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus,” 

Composites Science and Technology, Vol. 56, pp. 439-449, 1996. 

[9] N. Moës, J. Dolbow and T. Belytschko, “A finite element method for crack growth without 

remeshing,” Int. J. Numer. Meth. Engng. Vol. 46, pp. 131-150, 1999. 

[10] He, M-Y., and J. W. Hutchinson, “Kinking of a crack out of an interface,” Journal of 

Applied Mechanics, Vol 56, pp. 270-278, 1989. 

[11] M.J. Laffan, S.T. Pinho, P. Robinson and A.J. McMillan, “Translaminar fracture 

thougness: The critical notch tip radius of 0° plies in CFRP,” Composites Science and 

Technology, Vol. 72, No 1, pp. 97-102, 2011. 

[12] M. W. Czabaj and J. G. Ratcliffe, “Comparison of intralaminar and interlaminar mode-I 

fracture thoughness of unidirectional IM7/8552 graphite/epoxy composite,” Proceedings of 

the American Society for Composites – Twenty-Seventh Technical Conference, Arlington, 

TX, USA, paper 114, 2012 

[13] T. K. O'Brien, G. B. Murri and S. A. Salpekar, “Interlaminar shear fracture thoughness 

and fatigue thresholds for composite materials,” Composite Materials: Fatigue and Fracture, 

Vol. 2, ASTM STP 1012, Paul A. Lagace, Ed., American Society for Testing and Materials, 

Philadelphia, pp. 222-250, 1989. 

[14] T. K. O'Brien and R. Krueger, “Analysis of flexure tests for transverse tensile strength 

characterization of unidirectional composites,” Journal of Composites Technology and 

Research, Vol 25, pp. 50-68, 2003. 

[15] R. Krueger, “Development of and application of benchmark examples for mixed-mode 

I/II quasi-static delamination propagation predictions,” NASA/CR-2012-217562, 2012. 
12

Table 1. Elastic properties, IM7-8552 [14].









    















161.0 11.38 0.32 0.44 5.17 3.98 

Table 2. Critical energy release rates for delamination of IM7-8552 [15]. 











0.21 0.77 2.1 

Figure 1. Overview of the Floating Node Method (FNM). 

1 2

4 3

1 2

4 3 4 3

1 2

5 6

4 3 5

1 2

3

5

6

x

y

1 2

4 3

ξ

η

ξ

η η

= + ΩA

ΩB

x

y

x

y

x

Ω y

Ξ

5

6

5

6

Ξ Ξ

Ω

3

ΩA

x

y

6

4

ΩB

x

5

6

Without discontinuity With weak/strong discontinuity 

Nodal coordinates of 

initial integration domain 

Nodal coordinates of 

crack position 

Degree of freedom 

Floating degree of freedom

Key: 
13



(a) topology of the floating node element (b) two weak interfaces, e.g. due 

to plies with different 

orientations 

(c) delamination

(d) delamination with migration to matrix crack (e) matrix crack (f) matrix crack with 

migration to delamination 

Figure 2. Floating node element used showing different weak/strong discontinuities scenarios 



Figure 3. Virtual Crack Closure Technique and Floating Node Method for arbitrary crack propagation.

Figure 4. Bi-material interface. 

x

y

A 

B B 

Ω 

90° 

0° ΩA

90

ΩB

ΩC 0° 

ΩA

ΩB

ΩC

` 

Nodal coordinates of 

initial integration 

domain 

`Nodal coordinates of 

crack position 

Degree of freedom 

Floating degree of 

freedom 

Noda

Key: 

ΩA

ΩC

ΩB2

ΩB1

ΩA

ΩC

ΩB2

ΩB1

0◦ ply

90◦ ply

ΩA

ΩC

ΩB2

ΩB1

q

F 

a1

a2

renr

Interpolation 

Material B 

Material A

ω

y 

x 
14

300 

Clamp 

Clamp 

a0=49 12.7 

12.7 

L 

S = 115 Baseplate 

Specimen 

[0o]3

[90o]4

[0o] 

[90o]4

[0o

]2

Applied 

load 

0

o

h=5.25 

All units in mm. Not drawn to scale 

Specimen width, B=12.7mm 

90o-specimen width direction 

0o

-speciment span direction

Figure 5. Schematic of delamination migration test configuration [7]. 
15

Figure 6. Finite element model, and boundary conditions of migration test specimen [7]. 

12.7 mm 

[904/03/(90/0)2s/02/0/904/ T /0/904/0/0/(90/0)2s/0/0/903 [90 /0 /(90/0) /0 /0 /0/90] 

Stacking 

12.7 mm 

a0 = 49 mm 

5.5 mm 

0/904/ T /0

Plies modeled within the 

FNM element 

FNM 

u2 = 

ΩA

ΩB

ΩC

0°

90°

0°

115 mm 

u L 1 = u2 = 0

u1 = u2 = 0

u1 = u2 = 0

u1 = u2 = 0

= Δ 
16

(a) load offset: , (b) load offset: , 

(c) load offset: , (d) load offset: , 

Figure 7. Experimental and predicted load-displacement plots. 

0 

50 

100 

150 

200 

250 

300 

350 

0 0.5 1 1.5 

Force, 

N 

Displacement, mm 

Experiments 

	


Migration 

0 

50 

100 

150 

200 

250 

300 

350 

0 0.5 1 1.5 

Force, 

N 

Displacement, mm 

Experiments 

Simulation 

Migration 

0 

50 

100 

150 

200 

250 

300 

350 

400 

450 

0 0.5 1 1.5 

Force, 

N 

Displacement, mm 

Experiments 

Simulation 

Migration 

0 

50 

100 

150 

200 

250 

300 

350 

400 

450 

0 0.5 1 1.5 

Force, 

N 

Displacement, mm 

Experiments 

Simulation 

Migration 

s a0 

L=a0

a0 

L=1.1a0

a0

L=1.2a0

a0

L=1.3a0
17

(a) qualitative comparison between simulated and 

experimental crack paths 

(b) variation of with normalized load offset 

Figure 8. Crack path and migration location. 

Figure 9. Comparison of experimentally observed crack path for , and simulated crack 

paths for different load-offsets, . 

0 

2 

4 

6 

8 

10 

12 

14 

16 

1 1.1 1.2 1.3 

ΔK, 

mm 

Load offset, L/a0

Simulations 

Experiments 

0.0 

0.1 

0.2 

0.3 

0.4 

0.0 0.5 1.0 1.5 

mm 

mm 

L=a0 L=1.1a0 

L=1.2a0 L=1.3a0 

a0 

L ΔK 







REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Contractor Report

 4. TITLE AND SUBTITLE

Floating Node Method and Virtual Crack Closure Technique for Modeling 

Matrix Cracking-Delamination Interaction

5a. CONTRACT NUMBER

 6. AUTHOR(S)

De Carvalho, N. V.; Chen, B. Y.; Pinho, S. T.; Baiz, P. M.; Ratcliffe, J. G.; 

Tay, T. E.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center 

Hampton, Virginia 23681

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

Langley Technical Monitor: Thomas K. O'Brien

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 24

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

A novel approach is proposed for high-fidelity modeling of progressive damage and failure in composite materials that 

combines the Floating Node Method (FNM) and the Virtual Crack Closure Technique (VCCT) to represent multiple 

interacting failure mechanisms in a mesh-independent fashion. In this study, the approach is applied to the modeling of 

delamination migration in cross-ply tape laminates. Delamination, matrix cracking, and migration are all modeled using 

fracture mechanics based failure and migration criteria. The methodology proposed shows very good qualitative and 

quantitative agreement with experiments. 

15. SUBJECT TERMS

Composite materials; Crack closure; Delamination; Failure; Matrix cracking; Migration

18. NUMBER

 OF 

 PAGES

22

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

 NNL09AA00A

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

031102.02.07.05.9D46.12

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/CR-2013-218022

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 07 - 2013