REFERENCES

Andreotti, S., Franzoni, E., Degli Esposti, M., and Fabbri, P. (2018). Erratum: Andreotti, S.; Franzoni, E.; fabbri, P. Poly(hydroxyalkanoate)s-Based hydrophobic coatings for the protection of stone in cultural heritage. Materials 2018, 11, 165. Materials 11 (3), 389. doi:10.3390/ma11030389

Andreotti, S., Franzoni, E., Ruiz-Agudo, E., Scherer, G. W., Fabbri, P., Sassoni, E., et al. (2019). New polymer-based treatments for the prevention of damage by salt crystallization in stone. *Mater. Struct.* 52 (17), 1-28. doi:10.1617/s11527-018-1309-6

Angeli, M., Jean-Philippe, B., Benavente, D., Menéndez, B., Hébert, R., and David, C. (2007). Salt crystallization in pores: quantification and estimation of damage. *Environ. Geol.* 52 (2), 187-195. doi:10.1007/s00254-006-0474-z

Bracciale, M. P., Sammut, S., Cassar, J., Santarelli, M. L., and Marrocchi, A. (2020). Molecular crystallization inhibitors for salt damage control in porous materials: an overview. *Molecules* 25 (8), 1873. doi:10.3390/molecules25081873

Calia, A., Laurenzi Tabasso, M., Maria Mecchi, A., and Quarta, G. (2013). The study of stone for conservation purposes: Lecce stone (Southern Italy). *Geol. Soci., London, Special Publications* 391 (1), 139-156. doi:10.1144/sp391.8

Cassar, J., Marrocchi, A., Santarelli, M. L., and Muscat, M. (2008). Control of los daños por cristalización en la caliza de Malta mediante inhibidores de sales *Materiales de Costrucción* 58, 289-290. doi:10.3989/mc.2008.v58.i289-290.83

Charola, A. E., and Bläuer, C. (2015). Salts in masonry: an overview of the problem. *RBM.* 21 (4-6), 119-135. doi:10.1515/rbm-2015-1005

Charola, A. E. (2000). Salts in the deterioration of porous materials: an overview. *J. Am. Inst. Conserv.* 39 (3), 327-343. doi:10.1179/019713600806113176

Charola, A. E., and Wendler, E. (2015). An overview of the water-porous building materials interactions. *RBM.* 21 (2-3), 55-65. doi:10.1515/rbm-2015-2006

Doehne, E., and Price, C. A. (2010). *Stone conservation: an overview of current research.* 2nd Edn, Los Angeles, CA, United States: Getty Publications

EN 12370 (1999). *Natural stone test methods – determination of resistance to salt crystallization*, United Kingdom: BSI

EN 15801 (2009). *Conservation of cultural property – test methods – determination of water absorption by capillarity*, Santry Dublin, Ireland: NSAI

Espinosa-Marzal, R. M., and Scherer, G. W. (2010). Advances in understanding damage by salt crystallization. *Acc. Chem. Res.* 43 (6), 897-905. doi:10.1021/ar9002224

Espinosa-Marzal, R. M., and Scherer, G. W. (2009). Crystallization of sodium sulfate salts in limestone. *Environ. Earth Sci.* 56, 605-621. doi:10.1007/s00254-008-1441-7

Flatt, R. J., Caruso, F., Sanchez, A. M., and Scherer, G. W. (2014). Chemo-mechanics of salt damage in stone. *Nat. Commun.* 5, 4823. doi:10.1038/ncomms5823

Flatt, R. J. (2002). Salt damage in porous materials: how high supersaturations are generated. *J. Cryst. Growth* 242 (3-4), 435-454. doi:10.1016/s0022-0248(02)01429-x

Franceschini, M., Broggi, A., Bracciale, M. P., Sommei, L., Santarelli, M. L., and Marrocchi, A. (2015). Effectiveness of phosphocitrate as salt crystallization inhibitor in porous materials: case study of the roman mosaic ofOrpheus and the beasts(Perugia, Italy). *Int. J. Architect. Herit.* 9 (3), 195-200. doi:10.1080/15583058.2012.760121

Franzoni, E. (2014). Rising damp removal from historical masonries: a still open challenge. *Construct. Build. Mater.* 54, 123-136. doi:10.1016/j.conbuildmat.2013.12.054

Franzoni, E., Sassoni, E., and Graziani, G. (2015). Brushing, poulting or immersion? The role of the application technique on the performance of a novel hydroxyapatite-based consolidating treatment for limestone. *J. Cult. Herit.* 16 (2), 173-184. doi:10.1016/j.culher.2014.05.009

Franzoni, E. (2018). State-of-the-art on methods for reducing rising damp in masonry. *J. Cult. Herit.* 31, S3-S9. doi:10.1016/j.culher.2018.04.001

Goudie, A. S. (1986). Laboratory simulation of “the wick effect” in salt weathering of rock. *Earth Surf. Process. Landforms* 11 (3), 275-285. doi:10.1002/esp.3290110305

Granneman, S. J. C., Lubelli, B., and van Hees, R. P. J. (2019). Mitigating salt damage in building materials by the use of crystallization modifiers—a review and outlook. *J. Cult. Herit.* 40, 183-194. doi:10.1016/j.culher.2019.05.004

Granneman, S. J. C., Shahidzadeh, N., Lubelli, B., and van Hees, R. P. J. (2017). Effect of borax on the wetting properties and crystallization behavior of sodium sulfate. *CrystEngComm* 19 (7), 1106-1114. doi:10.1039/c6ce02163h

Graziani, G., Sassoni, E., and Franzoni, E. (2015). Consolidation of porous carbonate stones by an innovative phosphate treatment: mechanical strengthening and

physical-microstructure compatibility in comparison with TEOS-based treatments. *Heritage Science* 3 (1), 1-6. doi:10.1186/s40494-014-0031-0

Graziani, G., Sassoni, E., Franzoni, E., and Scherer, G. W. (2016). Hydroxyapatite coatings for marble protection: optimization of calcite covering and acid resistance. *Appl. Surf. Sci.* 368, 241-257. doi:10.1016/j.apsusc.2016.01.202

Gupta, S., Terheiden, K., Pel, L., and Sawdy, A. (2012). Influence of ferrocyanide inhibitors on the transport and crystallization processes of sodium chloride in porous building materials. *Cryst. Growth Des.* 12 (8), 3888-3898. doi:10.1021/cg3002288

Houck, J., and Scherer, G. W. (2006). Controlling stress from salt crystallization. *Fracture and failure of natural building stones*. Dordrecht, The Nederlands: Springer, 299-312. Chap 5

Lubelli, B., Cnudde, V., Diaz-Goncalves, T., Franzoni, E., van Hees, R. P. J., Ioannou, I., et al. (2018). Towards a more effective and reliable salt crystallization test for porous building materials: state of the art. *Mater. Struct.* 51 (2), 1-21. doi:10.1617/s11527-018-1180-5

Lubelli, B., and van Hees, R. P. J. (2007). Effectiveness of crystallization inhibitors in preventing salt damage in building materials. *J. Cult. Herit.* 8 (3), 223-234. doi:10.1016/j.culher.2007.06.001

Muxika, A., Etxabide, A., Uranga, J., Guerrero, P., and de la Caba, K. (2017). Chitosan as a bioactive polymer: processing, properties and applications. *Int. J. Biol. Macromol.* 105 (2), 1358-1368. doi:10.1016/j.ijbiomac.2017.07.087

Naidu, S., and Scherer, G. W. (2014). Nucleation, growth and evolution of calcium phosphate films on calcite. *J. Colloid Interface Sci.* 435, 128-137. doi:10.1016/j.jcis.2014.08.018

Naidu, S., Sassoni, E., and Scherer, G. W. (2011). New treatment for corrosion-resistant coatings for marble and consolidation of limestone in Jardins de Pierres – Conservation of stone in Parks, Gardens and Cemeteries, Paris, France, June 2011, 289-294

Ocak, Y., Sofuoğlu, A., Tihminlioglu, F., and Böke, H. (2009). Protection of marble surfaces by using biodegradable polymers as coating agent. *Prog. Org. Coating* 66 (3), 213-220. doi:10.1016/j.porgcoat.2009.07.007

Ocak, Y., Sofuoğlu, A., Tihminlioglu, F., and Böke, H. (2015). Sustainable bio-nano composite coatings for the protection of marble surfaces. *J. Cult. Herit.* 16 (3), 299-306. doi:10.1016/j.culher.2014.07.004

Pedna, A., Pinho, L., Frediani, P., and Mosquera, M. J. (2016). Obtaining SiO₂ 2-fluorinated PLA bionanocomposites with application as reversible and highly-hydrophobic coatings of buildings. *Prog. Org. Coating* 90, 91-100. doi:10.1016/j.porgcoat.2015.09.024

Philibert, T., Lee, B. H., and Fabien, N. (2017). Current status and new perspectives on chitin and chitosan as functional biopolymers. *Appl. Biochem. Biotechnol.* 181 (4), 1314-1337. doi:10.1007/s12010-016-2286-2

Rivas, T., Alvarez, E., Mosquera, M. J., Alejano, L., and Taboada, J. (2010). Crystallization modifiers applied in granite desalination: the role of the stone pore structure. *Construct. Build. Mater.* 24 (5), 766-776. doi:10.1016/j.conbuildmat.2009.10.031

Rodriguez-Navarro, C., Doehne, E., and Sebastian, E. (2000). How does sodium sulfate crystallize? Implications for the decay and testing of building materials. *Cement Concr. Res.* 30 (10), 1527-1534. doi:10.1016/s0008-8846(00)00381-1

Rodriguez-Navarro, C., Linares-Fernandez, L., Doehne, E., and Sebastian, E. (2002). Effects of ferrocyanide ions on NaCl crystallization in porous stone. *J. Cryst. Growth* 243 (3-4), 503-516. doi:10.1016/s0022-0248(02)01499-9

Ruiz-Agudo, E., Lubelli, B., Sawdy, A., van Hees, R., Price, C., and Rodriguez-Navarro, C. (2010). An integrated methodology for salt damage assessment and remediation: the case of San Jerónimo Monastery (Granada, Spain). *Environ Earth Sci* 63 (7), 1475-1486. doi:10.1007/s12665-010-0661-9

Ruiz-Agudo, E., Putnis, C. V., Pel, L., and Rodriguez-Navarro, C. (2012). Template-assisted crystallization of sulfates onto calcite: implications for the prevention of salt damage. *Cryst. Growth Des.* 13 (1), 40-51. doi:10.1021/cg300744x

Ruiz-Agudo, E., Rodriguez-Navarro, C., and Sebastian-Pardo, E. (2006). Sodium sulfate crystallization in the presence of phosphonates: implications in ornamental stone conservation. *Cryst. Growth Des.* 6 (7), 1575-1583. doi:10.1021/cg050503m

Sammut, S., Cassar, J., Marrocchi, A., Santarelli, M. L., Russo, C., and Sommei, L. (2014). Investigating a method to limit damage in Globigerina limestone, a soft porous stone widely used in historic buildings. *Proceedings of SWBSS 2014* 3rd international conference on salt weathering of buildings and stone sculptures, Brussels, Belgium, 14-16 October 2014. Royal Institute for Cultural Heritage Brussels, 101-115
Sandrolini, F., and Franzoni, E. (2007). Repair Systems for the Restoration of Ancient Buildings—dampness Rise Problem/Instandsetzungssysteme für das Restaurieren historischer Gebäude—aufsteigende feuchtigkeit. *RBM*. 13 (3), 161–172. doi:10.1515/rbm-2007-6129

Sassoni, E., Franzoni, E., Stefanova, M., Kamenarov, Z., Scopece, P., and Verga Falzacappa, E. (2020). Comparative study between ammonium phosphate and ethyl silicate towards conservation of prehistoric paintings in the magura cave (Bulgaria). *Coatings* 10 (3), 250. doi:10.3390/coatings10030250

Sassoni, E., Graziani, G., and Franzoni, E. (2016a). An innovative phosphate-based consolidant for limestone. Part 1: effectiveness and compatibility in comparison with ethyl silicate. *Construct. Build. Mater.* 102 (1), 918–930. doi:10.1016/j.conbuildmat.2015.04.026

Sassoni, E., Graziani, G., and Franzoni, E. (2016b). An innovative phosphate-based consolidant for limestone. Part 2: durability in comparison with ethyl silicate. *Construct. Build. Mater.* 102 (1), 931–942. doi:10.1016/j.conbuildmat.2015.10.202

Sassoni, E., Graziani, G., Franzoni, E., and Scherer, G. W. (2018). Calcium phosphate coatings for marble conservation: influence of ethanol and isopropanol addition to the precipitation medium on the coating microstructure and performance. *Corrosion Sci.* 136, 255–267. doi:10.1016/j.corsci.2018.03.019

Sassoni, E. (2018). Hydroxyapatite and other calcium phosphates for the conservation of cultural heritage: a review. *Materials* 11 (4), 1–48. doi:10.3390/ma11040557

Sassoni, E., Naidu, S., and Scherer, G. W. (2011). The use of hydroxyapatite as a new inorganic consolidant for damaged carbonate stones. *J. Cult. Herit.* 12 (4), 346–355. doi:10.1016/j.culher.2011.02.005

Scherer, G. W. (1999). Crystallization in pores. *Cement Concr. Res.* 29 (8), 1347–1358. doi:10.1016/s0008-8846(99)00002-2

Scherer, G. W. (2004). Stress from crystallization of salt. *Cement Concr. Res.* 34 (9), 1613–1624. doi:10.1016/j.cemconres.2003.12.034

Scherer, G. W., and Wheeler, G. S. (2008). Silicate consolidants for stone. *Kemi* 391, 1–25. doi:10.4028/www.scientific.net/kem.391.1

Selwitz, C., and Doehne, E. (2002). The evaluation of crystallization modifiers for controlling salt damage to limestone. *J. Cult. Herit.* 3 (3), 205–216. doi:10.1016/s1296-2074(02)01182-2

Smith, B. J., Gomez-Heras, M., Viles, A., and Cassar, J. (2010). *Limestone in the built environment: present-day challenges for the preservation of the past*. London, United Kingdom: Geological Society of London, 331

Steiger, M., and Asmussen, S. (2008). Crystallization of sodium sulfate phases in porous materials: the phase diagram Na2SO4-H2O and the generation of stress. *Geochem. Cosmochim. Acta* 72 (17), 4291–4306. doi:10.1016/j.gca.2008.05.053

Tsui, N., Flatt, R. J., and Scherer, G. W. (2003). Crystallization damage by sodium sulfate. *J. Cult. Herit.* 4 (2), 109–115. doi:10.1016/s1296-2074(03)00022-0

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 Bassi, Sassoni and Franzoni. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Mechanics of Materials

<figure><img src="image_196.png" /><figcaption>Laure Bouquerel</figcaption></figure>

Dr. Laure Bouquerel graduated in 2016 from the Engineering school of Ecole Centrale de Nantes, France, with a double diploma from Cranfield University, United Kingdom, in Advanced Materials. In 2019, she defended her Ph.D. at Mines Saint-Etienne, France, funded by Hexcel Reinforcements and in collaboration with INSA-Lyon, France. During her Ph.D., she focused on modeling the mechanical behavior of dry carbon-fiber reinforcements during forming for primary aircraft structure applications. Since 2019, she has been working at Hexcel Reinforcements as an R&T Engineer and has been involved in dry carbon-fiber reinforcement characterization as well as material development.

Declaration by Dr. Larissa Gorbatikh: "The corresponding author, Dr. Larissa Gorbatikh, would like to instead highlight the contribution of the manuscript's first author, Dr. Man Zhu."

<figure><img src="image_197.png" /><figcaption>Man Zhu</figcaption></figure>

Dr. Man Zhu was born in China and is currently a PostDoc student at Nanjing University of Aeronautics and Astronautics, China. She received her Master's degree in Solid Mechanics in 2015 from Beihang University, China, and her Ph.D. in Materials Engineering in 2019 from KU Leuven, Belgium. Her research interest is in fatigue and fracture analysis of composite materials.

<figure><img src="image_198.png" /><figcaption>Stephanie TerMaath</figcaption></figure>
Dr. Stephanie C. TerMaath is an Associate Professor and a Zeanah Faculty Fellow at the University of Tennessee, Knoxville, in the Mechanical, Aerospace, and Biomedical Engineering Department. Her research group develops numerical models to solve multi-physics problems. Before joining the faculty, she managed a physics-based computing group at Applied Research Associates (ARA) and held positions at Boeing Phantom Works and Lockheed Martin Aeronautics, where she worked on numerous programs, including the High Speed Civil Transport, 737, 777, C-17, Delta IV, and Joint Strike Fighter. She currently leads an afterschool program to teach STEM to girls at socioeconomically disadvantaged elementary schools.

<figure><img src="image_199.png" /><figcaption>Evelyn KF Yim</figcaption></figure>

Dr. Evelyn Yim received her Ph.D. in Biomedical Engineering at the Johns Hopkins University before her post-doctoral training at Duke University. Between 2007 and 2015, Evelyn was a faculty member with joint appointments in the Departments of Biomedical Engineering and Surgery at the National University of Singapore and as a principal investigator in the Mechanobiology Institute Singapore. Evelyn joined the Department of Chemical Engineering at the University of Waterloo in 2016. Experienced with nanofabrication technologies and stem cell culture, Evelyn and her group are interested in applying their knowledge of biomaterial-stem cell interaction to direct stem cell differentiation and tissue engineering applications.
<figure><img src="image_200.png" /></figure>

Check for
updates

# A Contribution to the Study of the Forming of Dry Unidirectional HiTape® Reinforcements for Primary Aircraft Structures

Laure Bouquerel¹*, Nicolas Moulin² and Sylvain Drapier²

¹R&T Reinforcements, Hexcel Reinforcements SASU, Route des Nappes, Les Avenières, France, ²Mines Saint-Etienne, University of Lyon, CNRS, UMR 5307 LGF, Centre SMS, Saint-Etienne, France

In the context of developing competitive liquid composites molding processes for primary aircraft structures, modeling the forming stage of automatically-placed initially flat stacks of dry reinforcements is of great interest. In the case of HiTape®, a dry unidirectional carbon fiber reinforcement designed to achieve performances comparable to state-of-the-art pre-impregnated materials, the presence of a thermoplastic veil on each side of the material for both processing and mechanical purposes should also be considered when modeling forming in hot conditions. As a dry unidirectional reinforcement, HiTape® is expected to exhibit a transversely isotropic behavior. Computation cost and strong characterization challenges led us to model its behavior at the forming process temperature (above the thermoplastic veil melting temperature) through a homogeneous equivalent continuous medium exhibiting four 'classical' deformation modes and a specific structural mode, namely out-of-plane bending. The response of both single plies and stacks of HiTape® to this latter structural mode was characterized at the forming process temperature using a modified Peirce flexometer. Results on single plies showed a non-linear softening moment-curvature behavior and a corresponding flexural stiffness much lower than what can be inferred from continuum mechanics. Moreover, testing stacks revealed that the veil acts as a thin load transfer layer between the plies undergoing relative in-plane displacement, i.e. inter-ply sliding. This inter-ply response was then characterized separately at the forming process temperature thanks to a specific method relying on a pull-through test. Experiments performed at pressures and speeds representative of the forming stage revealed that a hydrodynamic lubricated friction regime predominates, i.e. a linearly increasing relationship between the friction coefficient and the modified Hersey number. From an industrial point of view, high forming pressures and low speeds are therefore recommended to promote inter-ply slip to limit the occurrence of defects such as wrinkles.

**Keywords:** bending, carbon fiber, dry, forming, friction, reinforcement, thermoplastic, unidirectional

**OPEN ACCESS**

**Edited by:**
Patricia Krawczak,
IMT Lille Douai, France

**Reviewed by:**
Damien Soulat,
ENSAIT, France
Frédéric Charles Lebon,
Aix-Marseille Université, France

**\*Correspondence:**
Laure Bouquerel
laure.bouquerel@hexcel.com

**Specialty section:**
This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials

**Received:** 11 June 2020
**Accepted:** 04 November 2020
**Published:** 27 January 2021

**Citation:**
Bouquerel L, Moulin N and Drapier S (2021) A Contribution to the Study of the Forming of Dry Unidirectional HiTape® Reinforcements for Primary Aircraft Structures. Front. Mater. 7:571779.
doi: 10.3389/fmats.2020.571779
# INTRODUCTION

For decades, composite materials have been recognized as optimal solutions for lightweight structures in the industry of transportation. They have widely participated to the development of innovative aircraft structures, including heavily loaded primary structures. Carbon fiber based composites are now used extensively and represent up to more than half of the weight of modern aircraft structures. Classically, the best performances are reached with long fiber semi-products pre-impregnated with non-consolidated organic resin, and processed in autoclave environments. However, such semi-products, made up of the controlled arrangement of tows containing thousands of carbon fibers require heavy storage and processing environments. In the context of continuous cost reduction and production rate increase, liquid composites molding (LCM) processes (Hexcel, 2015a) appear now as ideal routes for manufacturing primary aircraft structures. These Out-of-Autoclave (OOA) processes consist in draping and shaping dry composite semi-products, and then inject (for Resin Transfer Molding for instance) or infuse (in infusion-based processes) liquid resins in this potentially complex fibrous architecture.

In this context, Hexcel Reinforcements is continuously developing specific reinforcements solutions for primary aircraft structures manufactured with LCM processes. HiTape® (Hexcel, 2015b) is a dry unidirectional (UD) tape made up of parallel spread tows of thousands of carbon fibers with a thermoplastic (TP) veil on either side (Figure 1A). The veil is a key feature of this innovative reinforcement; not only it enables automated dry fiber placement (DFP) but it also contributes to yield mechanical properties equivalent to pre-impregnated based composites (Hexcel, 2015b), with high fiber volume fractions of 58–60% reached with OOA processes.

Going a step further in increasing production rates, dry semi-products such as HiTape® can be first deposited with automated placement in 2D stacks of several plies, and then formed in an automated way under mechanical loading to reach a given geometry and fiber arrangement. In this paper, the considered forming process solution is double-membrane hot vacuum forming as illustrated in Figure 1B. It consists firstly in compacting the dry stack between soft membranes by pulling vacuum (①), and secondly in forcing this membranes/stack set to deform (Pickett, 2018) by pulling vacuum between the mold and the membranes/stack set yielding an about 1 bar forming pressure (②). This forming step is performed above the melting temperature of the veil, notably for the veil to maintain the shape of the resulting preform. Preventing defects and imperfections (e.g., wrinkles) occurrence during this crucial forming stage is of prime importance and requires to understand the ply response to thermo-mechanical loading, and how plies move with respect to each other.

The aim of this paper is to characterize the behavior of HiTape® reinforcements in order to feed models of the forming process of stacks. To the knowledge of the authors, literature on the characterization of such dry UD materials is quite sparse. In Section 2, we hence propose a literature review on the behavior of related materials such as dry multi-axial reinforcements and pre-impregnated UDs, in order to define the major deformation modes to be identified for forming modeling, along with the corresponding appropriate experimental methods. Such review allows us to identify the most relevant deformation modes of HiTape® to be characterized: bending and inter-ply sliding. In Section 3, the implementation of corresponding dedicated experimental methods is described. Finally, results of HiTape® characterization in bending and inter-ply sliding are presented and discussed in Section 4.

<figure><img src="image_201.png" /></figure>
<figure><img src="image_202.png" /><figcaption>**FIGURE 2** | The three main scales of observation of composites materials; an intermediate (mesoscopic) ply scale is selected for forming modeling and HiTape® characterization.</figcaption></figure>

## FROM STATE OF THE ART TO CONSIDERATIONS ON HITAPE® REINFORCEMENTS MECHANICAL BEHAVIOR

Extensive work has been performed on the characterization of the forming behavior of composite semi-products with heavy multi-axial complex architectures, either for pre-impregnated (e.g., Ten Thije et al. (2011), Lightfoot et al. (2013)) or dry reinforcements (e.g., Boisse (2004), Cao et al. (2008), Senner et al. (2014)). However, the particular structure of HiTape® is expected to yield a specific intermediate type of behavior. On one hand, an individual ply (intra-ply) is likely to respond like pre-impregnated UD materials (Leutz, 2015), i.e., with a theoretically transversely isotropic behavior, but with even less cohesion and very low stiffness in the isotropy plane orthogonal to the fiber direction. On the other hand when considering several plies forming in industrial conditions (hence at the veil melting temperature), the presence of the melted veil at the interface between plies (inter-ply) may modify HiTape® stacks response compared to multi-axial dry semi-products (Creech and Pickett, 2006). The aim of this section is to perform a selective literature review in order to identify which characterization work should be carried out on HiTape® in order to describe its mechanical behavior to be used as input for forming simulations.

### Intermediate Scale of Observation

Reinforcements are intrinsically multi-scale materials as they are classically composed of thousands of fibers gathered in tows (possibly spread) which are assembled in tapes or plies, then laid up as a stack (Figure 2). The scale of observation for this study must allow for the identification of the main mechanisms which control both intra-ply and inter-ply responses during forming, but should also enable modeling the process with tractable computational efforts.

Modeling the forming at fiber scale would help understanding interactions between fibers, but it is not suited for representing the forming of stacks which contain millions of fibers with specific architectures, nor accounting for the contact between top and bottom plies with the surrounding tooling. On the opposite, the part (or stack) scale is too coarse as heterogeneous layer-wize constitutive responses are expected. The material response model should indeed be able to capture local effects (inter-ply shear, frame rotations...) which control the relative motion between plies as well as the intra-ply response itself. It should also represent the large deformations undergone by the whole stack, which translate into finite strains and rotations that a continuous approach can certainly not handle even with refined kinematics (Dufort et al., 2001). Eventually, an intermediate scale, so-called mesoscopic, is considered for modeling HiTape® forming as well as to characterize the ply and stack responses. Indeed, as illustrated in Figure 2, modeling each ply through an equivalent homogeneous medium (potentially non-linear and orthotropic) ensures acceptable computational efforts (Fagiano, 2010; Dörr et al., 2017) while allowing for a reasonably fine description of the mechanisms acting in the ply and between plies (Bouquerel et al., 2017).

### Deformation Modes at ply Scale During Forming and Corresponding Characterization Methods

As illustrated in Figures 3A,B, reinforcements exhibit both in-plane and out-of-plane (i.e. along thickness direction) intra-ply deformation mechanisms: traction, compression as well as shear. Moreover, like any fibrous medium, HiTape® may also display structural modes like out-of-plane bending due to the reinforcement intrinsic thickness. Finally, these intra-ply modes must be complemented by specific ones associated with the presence of the TP veil on either side of the layers, which come into play when considering HiTape® stacks. These inter-ply
<figure><img src="image_203.png" /><figcaption>**FIGURE 3** | Deformation mechanisms of reinforcements at ply scale during forming, illustrated in the case of a UD with fibers oriented along $E_1$ direction: (A) in-plane modes, (B) out-of-plane modes with specific bending about axis $E_2$: $M = B\kappa_2$, and (C) inter-ply loading with specific inter-ply sliding in plane $(O, E_1, E_2)$, with friction coefficients $\gamma_{1,2}$.</figcaption></figure>

modes correspond to sliding or opening/closing as presented in **Figure 3C**.

In this section let us consider these intra- and inter-ply deformation modes and how the corresponding responses can be identified. We shall not describe all these modes in detail, but focus on the predominant ones during forming.

## Tensions and Compressions

The longitudinal response along fiber direction $E_1$ (**Figure 3A**) is dominated by fibers quasi-inextensibility. Tension behavior can hence be approximated by a simple rule of mixture $E_L \approx V_f E_f$ where $E_L$ is the reinforcement longitudinal modulus, $V_f$ the fiber volume fraction and $E_f$ the fiber longitudinal modulus in tension known to be constant for carbon fibers and available in the fiber data sheet. As for compression in this same longitudinal direction, it is assumed that the external compaction that occurs in double-membrane forming constrains the transverse displacement on either side of the plies, and therefore largely limits the occurrence of potential out-of-plane fiber micro-buckling that may appear if layers are free (Drapier et al., 1996; Drapier and Wisnom, 1999). Moreover, due to the high volume fraction ($V_f \approx 60\%$ under 1 bar compaction as measured by Blais (2016)) and in-plane confinement of plies, in-plane fiber micro-buckling could only appear at extremely high loads (Rosen, 1964; Drapier et al., 1996) which cannot be reached in standard forming conditions. In conclusion, the prescribed transverse compaction during forming limits the potential difference in tension-compression response, and the longitudinal compressive modulus can hence be taken equal to the constant elongation modulus $E_L$.

In transverse directions $E_2$ and $E_3$, the characterization of tension behavior of UDs is quite a challenge due to the lack of cohesion. Pre-impregnated UDs were tested in transverse direction by Leutz (2015) and Margossian et al. (2016), but results were scattered with low repeatability. In the case of HiTape®, the internal structure with no resin support inside plies exhibits little cohesion only brought by a slight entanglement of the UD dry fibers network, whereas TP veil on either side of the tape brings a superficial cohesion just sufficient for handling and processing. This structure makes

both elongational responses in transverse directions almost impossible to characterize. Therefore, for modeling purpose, it can be assumed that transverse stiffnesses correspond to a little fraction of the longitudinal modulus, i.e. $E_T = \alpha E_L$ with $\alpha \ll 1$.

In transverse compression, the response of the fibrous material is driven by the fiber network rearrangement as well as fiber-to-fiber contacts (Xiong et al., 2019). Recalling that we consider double-membrane forming processes where the individual layers of the stack are compacted along thickness direction ($E_3$) before forming, compression behavior understanding appears to be of prime importance. Such non-linear behavior is already known as it was determined using a standard testing machine to compact HiTape® stacks between circular platens in previous works (Celle et al., 2008; Blais, 2016; Blais et al., 2017).

## Shears

Considering the small thickness of HiTape® (approximately 0.2 mm), we can assume that transverse shear–i.e. in planes $(O, E_2, E_3)$ and $(O, E_1, E_3)$ as depicted in **Figure 3B**–is negligible. Conversely, in-plane shear (**Figure 3A**) may develop during forming. This deformation mechanism indeed dominates the formability of composites semi-products such as woven ones (Cao et al., 2008) due to the relative mobility of the fiber tows (Boisse, 2004). Two types of specific characterization methods were hence developed to evaluate the formability of both dry and pre-impregnated semi-products: the bias-extension (Orawattanasrikul, 2006; Charmetant, 2011) and the picture-frame (Kawabata, 1980) tests.

To the knowledge of the authors, literature on the characterization of the in-plane shear behavior of dry UD reinforcements under forming conditions is quite sparse. Related studies may be found for UD-NCF (Trejo et al., 2020) as well as for pre-impregnated UD tapes with various strategies proposed to maintain the specimen integrity during testing: Leutz (2015) tested 4-ply stacks (**Figure 4**), McGuinness and Bradaigh (1998) used membranes to maintain a single ply but struggled to isolate the ply response from the membrane one, and Larberg et al. (2012) considered cross-plied. However, all these tests demonstrated a rapid decohesion following the occurrence of
<figure><img src="image_204.png" /><figcaption>**FIGURE 4** | An example of in-plane shear characterization of pre-impregnated UD carbon/epoxy: picture-frame test on a 4-ply stack of UD of same orientation (left to right: initial configuration, intermediate shear angle of 2°, final shear angle of 40°) (Leutz, 2015).</figcaption></figure>

out-of-plane instabilities, either fiber micro-buckling or ply buckling. Due to the only superficial cohesion of HiTape®, in-plane shear characterization seems out of reach from a technological point of view. More importantly, due to the compaction which limits the fiber transverse displacements during forming, the formability of a stack of HiTape® plies is expected to be driven more by the relative motion of thin HiTape® individual plies than by intra-ply in-plane shear as it is highlighted in the following.

## Out-of-Plane Bending

Many recent studies showed that accounting for ply scale out-of-plane bending behavior of either dry or pre-impregnated semi-products (Cao et al., 2008; Haanappel et al., 2014) in forming simulations allows the prediction of wrinkles occurrence in thick architectures (Boisse et al., 2011), and more generally is a step forward to reliable simulations (Madeo et al., 2015). Moreover, even for thin plies, the bending behavior controls the overall formability to a great extent (Liang et al., 2014). Although the hypothesis of continuity is often made for fibrous media—especially those with high fiber volume fraction—, bending behavior cannot be directly deduced from the elongation stiffness for these media, mainly because of internal relative motion between fibers during deformation (Kang and Yu, 1995). For this reason, reinforcements bending behavior should be characterized.

Due to the low bending stiffness of composite semi-products, classical characterization methods such as three- or four-point bending tests have classically been discarded by the textile industry. Peirce (Peirce, 1930) early proposed a ply scale linear elastic bending model for cloth, extending Euler-Bernoulli's beam theory to large deflections but small deformations using a correction term. This model led to a simple, cost-effective and rapid bending characterization method called fixed angle flexometer, and involving a strip of fabric bending under its own weight. Using this method Leutz (2015) measured the bending stiffnesses $B_i$ about direction $E_i$ ($i = 1, 2$) of pre-impregnated UD's in thermally controled environments and showed that the ratio of bending stiffnesses $B_2/B_1$ is about $10^3$. Indeed, bending about $E_2$ is dominated by the fiber

longitudinal stiffness whereas bending about $E_1$ is almost solely controled by the molten resin as illustrated in **Figure 3B**. In the case of HiTape® this bending ratio is expected to be even higher due to the absence of molten resin. We can therefore infer that $B_2$ ($B$ thereafter) dominates in forming modeling, i.e. the influence of $B_1$ on forming simulation may be negligible. Still, the question of how to measure bending stiffness $B_1$ remains open for future works.

Going a step forward, fibrous semi-products were shown to exhibit a non-linear bending response (Liang et al., 2014). Their tangent behavior can be more generally defined by the differential form $B(\kappa) = dM(\kappa)/d\kappa$ in the corresponding plane, where $M$ is the bending moment per unit width, $B$ the bending stiffness per unit width, and $\kappa$ the curvature (Clapp et al., 1990). Following Peirce (1930), various models were introduced to account for the non-linearity of the moment-curvature relationship: bi-linear (Huang, 1979), non-linear for small curvatures and linear for larger curvatures (Abbott et al., 1971; Clapp et al., 1990), with a threshold (Grosberg, 1966), etc.

In order to characterize this particular behavior, a specific apparatus called Kawabata Evaluation System KES-FB2 (Kawabata, 1980) was designed and successfully used on fabrics (Ngo Ngoc et al., 2002) and non-crimp fabric (NCF) multi-plies (Lomov et al., 2003) for example; however this type of apparatus is rather cumbersome and limited to room temperature. Other characterization methods were further derived from the fixed angle flexometer test to characterize non-linear behaviors: de Bilbao et al. (2008) increased step by step the free length of the strip for NCFs and interlocks, while other approaches were based on recording the deformed shape of the sample and using image processing to relate the bending moment to the sample curvature. In order to characterize a wide range of curvatures, Liang et al. (2014) also added a lumped mass at the end of the strip.

Besides, some methods were developed to take into account the temperature-dependency of the bending behavior of pre-impregnated reinforcements, for UD carbon/PEEK<fnref n="1" /> and UD

<footnote n="1">¹PEEK: PolyEtherEtherKetone</footnote>
<figure><img src="image_205.png" /><figcaption>FIGURE 5 | (A) Principle of the pull-through test and (B) schematic of the inter-ply friction characterization setup.</figcaption></figure>

carbon/PPS<fnref n="2" /> for example (Liang, 2016). In those studies, managing the increase of temperature in the sample before testing was reported to be challenging. A specific apparatus called hot vertical cantilever (Soteropoulos et al., 2011; Alshahrani and Hojjati, 2017) was hence proposed in order to load the sample only when the desired temperature was reached. However, a deformation of the sample was already observed before applying the load due to the temperature rise (Angel and Graef, 2016). In addition, the homogeneity and control of the environment are critical and can significantly influence the results (Angel and Graef, 2016). To overcome this, other tests were proposed on smaller volumes: for instance in a 3-point bending DMA machine (Margossian et al., 2015) or with a rheometer (Sachs et al., 2014).

## Inter-ply Response

The two inter-ply deformation modes are depicted on **Figure 3C**, namely inter-ply opening (decohesion) and sliding. In the case of double-membrane forming processes, it seems wise to assume that decohesion along direction $E_3$ is not likely to occur due to external loading. Conversely, sliding between plies is expected to be a key factor controlling the formability of a stack (Guzman-Maldonado et al., 2019); inter-ply sliding needs therefore to be characterized, and is dominated by the presence of the TP veil.

## Characterization Methods of Inter-Ply Sliding

Standard methods were developed for the characterization of the frictional behavior of sheets of materials-standards related to inter-fiber or inter-tow friction are not discussed here. They involve a fixed plane, either horizontal or inclined, onto which a first sample of tested material is attached, and over which a mobile block holding a second sample is translated. The static (i.e. when relative motion is initiated) and dynamic (i.e. during relative motion) friction coefficients are calculated using

Coulomb's law, which assumes that friction coefficients are independent of both relative sliding velocity and normal pressure:

$$CoF = \frac{F_T}{F_N} \qquad (1)$$

where CoF is the friction coefficient, and $F_T$ and $F_N$ are the applied tangent and normal forces, respectively.

However, these methods have been developed to represent friction conditions that sheet materials may encounter during handling or automated processing, but they are not representative of composite forming conditions. Indeed, not only are they limited to room temperature, but the required force corresponding to forming (i.e. of the order of magnitude of a few hN for standard tests) cannot be applied with such methods. Those are the reasons why other tests were developed for composite semi-products friction characterization. A horizontal plane with a specific system for applying normal pressure was proposed by Hivet et al. (2012), Montero et al. (2017), but the pressure still seems too low for our forming conditions. To achieve normal pressures of about 1 bar, the *pull-out* and the *pull-through* tests (illustrated on **Figure 5A**) are the most common ones; they are described in the benchmark of Sachs et al. (2012). In these tests, a reinforcement sample (referred to as central specimen) is placed between two plates covered with reinforcement samples (referred to as lateral specimens) and a normal force is applied. The central sample is then loaded in tension and forced to slide, while the corresponding force is recorded. The difference between the two apparatuses lies in the contact area, which is constant during the test for the pull-through while it decreases for the pull-out.

## Inter-ply Friction Behavior of Reinforcements

In the literature, Coulomb's law (Eq. 1) for dry friction is generally used to model the friction behavior during inter-ply sliding of solid-state reinforcements, i.e. either dry (Hivet et al., 2012) or containing thermoplastic below its melting temperature (Sachs et al., 2012) as illustrated in **Figure 6B** (top). Indeed, after a transient phase, the coefficient of friction reaches a plateau value often interpreted as the dynamic coefficient of friction. However,

<footnote n="2">²PPS: PolyPhenylene Sulfide</footnote>
<figure><img src="image_206.png" /></figure>

the hypothesis of a single constant friction coefficient for solid-state reinforcements is too coarse; normal pressure in particular has a strong influence on the friction coefficient Hivet et al. (2012). Thus more complex models were proposed; for example Howell's model (Howell and Mazur, 1953) describes a power law between tangent and normal forces such that $F_T = k F_N^m$ where k and m are parameters to be determined experimentally (Das et al., 2005; Cornelissen et al., 2014; Najjar et al., 2014), which leads to a dynamic friction coefficient which depends on the normal force. A review of these models is proposed by Yuksekkay (2009). In the case of dry reinforcements, an inter-ply friction coefficient of 0.3 is frequently assumed (Fetfatsidis et al., 2013) in the absence of experimental data. This is consistent, for example, with the results obtained by Hivet et al. (2012) who identified a friction coefficient between 0.15 and 0.4 for a dry woven reinforcement.

As soon as the reinforcement contains resin in a molten state, this resin may play the role of a lubricant and strongly influences the friction response; this is called lubricated contact. In that case, the friction coefficient depends on resin viscosity η, sliding speed v and normal force applied $F_N$; three parameters assembled in the Hersey number H defined by:

$$H = \frac{\eta v}{F_N} \qquad (2)$$

Based on a phenomenological approach, Stribeck's theory (Stribeck, 1903) expresses the existence of a minimum on the curve of friction coefficient vs. Hersey number. Three different lubrication regimes illustrated on **Figure 6A** can then be highlighted. In the limit regime, the lubricating film is very thin, i.e. its thickness is of the same order of magnitude as the roughness of the two surfaces; in this case the friction behavior depends both on the lubricant and the physico-chemical properties of the surfaces. In the mixed regime, there co-exist some areas where the two surfaces are in direct contact and some areas where they are separated by a lubricating film (Frêne and Zaïdi, 2011). In the hydrodynamic regime, the two surfaces are completely separated by a lubricating film; the roughness of the

surfaces is then negligible. In this case, the relationship between the friction coefficient and the Hersey number is increasing linearly (Hersey, 1914). Stribeck's theory was frequently used for either thermoplastic (Gorczyca et al., 2007; Ten Thije et al., 2011; Vanclouster et al., 2010; Fetfatsidis et al., 2013; Zhu et al., 2011; Wang et al., 2013) or thermoset (Rashidi et al., 2020) pre-impregnated reinforcements friction characterization for forming modeling. In this case of reinforcements containing a resin in the fluid state, the same type of curve with stabilization of the friction coefficient is obtained as illustrated in **Figure 6B** (bottom), and a hydrodynamic regime is generally observed. Affine (Najjar et al., 2014) or power laws (Haanappel et al., 2014) for example can be used to relate the friction coefficient to the Hersey number; in other cases, the hydrodynamic regime is directly assumed (Bel et al., 2012). If the normal pressure is high, the resin at the interface may migrate and limit lubrication may occur: a mixed lubrication model is then more suitable (Larberg and Akermo, 2011; Leutz, 2015). For reinforcements containing molten thermoplastic, the friction coefficient varies over a wider range, for example values between 0.1 and 1.8 were measured by Sachs et al. (2012).

In both cases (presence or not of a resin in a fluid state), during the transient regime prior to stabilization of the tangential force, the latter can reach a maximum value (Figure 6B) which becomes more evident as the velocity increases (Ten Thije et al., 2011). If Fetfatsidis et al. (2013) used this maximum force value to deduce a static coefficient of friction, this interpretation is however not unanimous. Indeed, in the case of the work of Ten Thije et al. (2011), the tangent force peak did not correspond to the initiation of relative motion but to a displacement value already higher than the one matching the maximum elongation of the sample before rupture. This means that the sample had already started to slide, and this maximum coefficient of friction did not correspond to a static coefficient of friction. Instead, the transition from peak to plateau was interpreted as a rearrangement of the reinforcement.
<figure><img src="image_207.png" /><figcaption>FIGURE 7 | (A) Principle and notations for the bending moment calculation with Eq. 7 and (B) schematic of the bending test.</figcaption></figure>

Finally, a dependence of the friction coefficient upon fiber orientation was observed. Testing a carbon/epoxy UD, Leutz (2015) for example obtained a friction coefficient varying from 0.88 for orientations of 0°/90°<fnref n="3" /> to 1.37 for 0°/0°, with an intermediate value for 90°/90°. Larberg and Akermo (2011), who identified a mixed lubrication model, obtained on the same type of material a maximum coefficient of friction for 0°/ 45° and 0°/30°, about 25% higher than that obtained with orientations of 0°/90°and 0°/0°.

Conclusions on HiTape® Characterization
Strategy

The literature review conducted so far highlights the two
deformation mechanisms that primary require characterization in
order to feed forming simulations of HiTape®: bending and inter-ply
friction. In particular for inter-ply friction, the presence of the TP veil
(which differs from a continuous film) raises the question of which
of the mixed or hydrodynamic lubrication regime predominates.
Such characterization works need to be conducted under conditions
as close as possible from double-membrane vacuum hot forming
process, hence at the forming process temperature (which is above
the melting temperature of the TP veil) and for external pressures
between 0.5 and 1.5 bar. In the following, the selected methods for
both non-linear bending behavior and interply friction behavior
characterizations are introduced.

MATERIALS AND METHODS

Materials

The characterized HiTape®<fnref n="4" /> reinforcement was composed of Hexcel
HexTow® fibers with a thermoplastic veil on each side of the tape.
This tape was 12.7 mm wide and had an areal weight of 210 g/m².

Out-of-Plane Bending Characterization
Method

In order to characterize the non-linear bending behavior of HiTape®
at ply scale and at the forming process temperature (T<sub>proc</sub>) which is
above the melting temperature of the TP veil (T<sub>melt</sub>), a modified

flexometer based on the same principle than previous studies (de
Bilbao et al., 2010; Liang et al., 2014) was designed. The selected
method to identify the non-linear bending response consisted in
recording the deformed shape, calculating the bending moment over
the sample length, and deriving a bending moment-curvature
relationship. Such procedure was extended from previous studies
(Clapp et al., 1990; de Bilbao et al., 2008; 2010; Liang et al., 2014)
carried out on both pre-impregnated and dry materials.

Test Setup

The bending test apparatus, represented in **Figure 7B**, was composed of a frame with a specimen holder developed to limit edge effect due to the clamp. Thermocouples were used for temperature monitoring. This apparatus was placed in a thermoregulated oven equipped with a glass window. A CCD camera was positioned outside the oven and a LabVIEW program⁵ recorded in real time both photographs of the sample and temperature of the thermocouples, with a chosen acquisition frequency and duration. A set of masses to be positioned at the free end of the sample was available to increase sample curvature. The grip system to attach the selected mass to the sample tip justified to model it as a point force oriented with gravity.

Test Procedure

The sample (either unit plies or stacks obtained with DFP) was
first pre-consolidated under vacuum at T<sub>proc</sub> in order to be as
close as possible to double-membrane vacuum hot forming
process conditions. The sample was then placed in the oven at
T<sub>proc</sub>. When it no longer deformed, a photograph of the sample
was recorded and post-processed.

Data Post-processing

Using the output photograph of the deformed sample, the
cartesian coordinates (x, z) of the mean-line of the deformed
sample in the orthonormal basis (E₁, E₃) were obtained by
digital image processing. Then, a program ran the following steps
(corresponding notations are illustrated in **Figure 7A»):

1. Derive the relationship $\hat{z} = f(x)$ using a Levenberg-Marquardt algorithm (Marquardt, 1963) to find the optimal couple $(k_1, k_2)$

<footnote n="3">³Orientation of samples in contact, in degrees with respect to loading direction.</footnote>

<footnote n="4">⁴Permission from the product owner was obtained for this study.</footnote>

<footnote n="3">³www.ni.com/fr-fr/shop/labview.html</footnote>
<figure><img src="image_208.png" /><figcaption>FIGURE 8 | Comparison between experimental data (dashed lines) and model (solid lines) for a single ply tested at room temperature: (A) deformed sample mean line with err<sub>%</sub><sup>(z)</sup> = 1.3% and (B) normalized bending moment-curvature relationship with err<sub>%</sub><sup>(M)</sup> = 3.8%.</figcaption></figure>

<figure><img src="image_209.png" /></figure>

minimizing the distance between the experimental data (x, z) and the following function:

$$ \tilde{z} = f(x) = -k_1 \ln(1 + (k_2 x)^2) \quad (3) $$

which is referred to as deformed shape function thereafter. This particular form illustrated for instance in **Figure 8A** is justified in **Section 4.1.2.1**.

2. Compute curvature $κ(x)$ from the established relationship $\tilde{z} = f(x)$:

$$ \kappa(x) = \frac{f''(x)}{(1 + f'(x)^2)^{3/2}} \quad (4) $$

3. Compute curvilinear abscissa $s(x)$ as follows:

$$ s(x) = \int_{t=0}^{x} \sqrt{1 + f'(t)^2} dt. \quad (5) $$

4. Compute bending moment per unit width $M(s)$ for each point $A(s)$ of the deformed sample. To this end, a first analytical expression could be derived for a sample of width $b$ loaded under its own areal weight $A$ (ply areal weight multiplied by number of plies) under gravity $g$ (see also **Figure 7A**):

$$ M(s) = - \int_{\xi=s}^{L} \frac{AB(\xi) \wedge AgE_3}{\xi} d\xi \approx \int_{\xi=s}^{L} Ag(\xi - s)\cos(\alpha(\xi))d\xi \quad (6) $$

where $s$ is the curvilinear abscissa at point $A$ such that $s = 0$ at clamping $O$ and $s = L$ at free end $C$, and the curvilinear abscissa $ξ$ and angle $α$ are the coordinates of the current $B(ξ)$ point from $A(s)$ to $C$ in the associated Frénet's basis. The additional mass $m$ placed at the free end of the sample also contributed to the bending moment which then reads:

$$ M(s) \approx \int_{\xi=s}^{L} Ag(\xi - s)\cos(\alpha(\xi))d\xi + \frac{mg}{b} A(s)C \cdot E_1. \quad (7) $$

5. Derive a relationship between bending moment and curvature $\hat{M} = g(κ)$ using a Levenberg-Marquardt algorithm (Marquardt, 1963) to optimize the triplet $(R_0, R_{inf}, κ_{lim})$ by minimizing the distance between the moment $M$ computed from expression seven and Voce's model (Voce, 1955) which reads:

$$ \hat{M} = g(\kappa) = M_0 + R_0\kappa + R_{inf} (1 - e^{-\kappa/\kappa_{lim}}) \quad \text{with } M_0=0. \quad (8) $$

Voce's model (Pannier, 2006) is illustrated in **Figure 8B** and its justification is given in **Section 4.1.2.1**.

The outputs of the corresponding program were the following: the range of curvatures of the deformed sample, the optimized couple $(k_1, k_2)$ in Eq. 3 describing the shape of the deformed sample, the optimized triplet $(R_0, R_{inf}, κ_{lim})$ in Eq. 8 giving the bending moment vs. curvature relationship, and finally, in any point of the deformed sample: the curvilinear abscissa, the curvature and the bending moment.

The two successive regressions—firstly on coordinates $(x, z)$ in step 3.2.3 and secondly on points $(κ, M)$ in step 3.2.3—induced an uncertainty which could be quantified by the two following errors that describe the mean scatter between experimental data $(z_i, M_i)$ and their homologous from the corresponding model $(\tilde{z}_i(x_i), \hat{M}_i(κ_i))$ at the ith point of the deformed sample:

$$ \text{err}_{\%}^{(z)} = \frac{\sqrt{\sum_{i=1}^{N} (\hat{z}_i(x_i) - z_i)^2}}{|z|_{\text{max}} - |z|_{\text{min}}}; \quad \text{err}_{\%}^{(M)} = \frac{\sqrt{\sum_{i=1}^{N} (\hat{M}_i(κ_i) - M_i)^2}}{|M|_{\text{max}} - |M|_{\text{min}}}. \quad (9) $$

For commodity reasons, these errors were normalized by the range of values obtained experimentally. We could then define the two following validity criteria for any tested sample: $err_{\%}^{(z)} < 5\%$ and $err_{\%}^{(M)} < 5\%$. In some cases an angle locally formed on the samples (due to the clamping condition for
example); such phenomenon was related to the tendency of the quasi-inextensible fibers to out-of-plane buckle in the inner radius of the curvature, especially for thick stacks. In this case, the samples exhibited an angular mean-line with a poor quality of regression that the validity criteria allowed to reject. However it was not always sufficient and a complementary visual check of the quality of the regression was realized. Eventually, the coefficient of variation $\sigma_{\%}^{(a)}$ was of interest to measure the mean scatter of a set of N data $a_i$ around its mean $\bar{a}$:

$$ \sigma_{\%}^{(a)} = \frac{\sqrt{\sum_{i=1}^{N} (a_i - \bar{a})^2}}{\bar{a}}. \qquad (10) $$

## Inter-ply Friction Characterization

Based on previous studies (Gorczyca et al., 2007; Vancooster et al., 2010; Ten Thije et al., 2011; Zhu et al., 2011; Fetfatsidis et al., 2013; Wang et al., 2013), the pull-through test principle was selected to characterize HiTape® inter-ply lubricated friction at ply scale. Let us recall the principle of this test illustrated on Figure 5A: a central sample is placed in between two fixed lateral samples and a controled normal pressure is applied across the thickness of the reinforcements, the central sample is then pulled through the setup. The tensile traction force $F_T$ is recorded as a function of the central sample displacement. This test principle is selected because it enables to apply a normal pressure that is both homogeneous on the surface and constant during the test.

The test setup was designed such that conditions representative of double-membrane vacuum hot forming were covered. As forming pressure is about 1 bar-but may vary locally -, normal stresses in the range 0.5-1.5 bar were considered. The relative speed between two adjacent plies was tricky to estimate; in the case of vacuum forming, which is a slow forming process, a relative speed of the order of magnitude of a few millimeters in some seconds was considered, i.e., 1-10 mm min⁻¹. Finally, the forming process temperature $T_{proc}$ such that $T_{proc} > T_{melt}$ was prescribed. In this study, temperature was kept constant as it is considered that the stack is heated up to $T_{proc}$ before starting forming; in future works the influence of a variation of temperature around $T_{proc}$ may be studied.

## Test Setup

A fixture was designed both to hold the lateral specimens in place and to apply the desired normal pressure and temperature conditions. Figure 5B provides a schematic of the final bench. The closing pressure was prescribed by an actuator. The lateral samples were placed in contact with heating pads and were fixed by means of a specific system. This apparatus was designed to be mounted on a standard traction machine, and the central sample attachment was fixed to the crosshead of the machine. The test was controled using the traction machine software, and an acquisition station enabled the force, displacement, and time data of the traction machine to be recorded in real time, as well as the temperature of the heating pads.

Temperature homogeneity inside the samples was checked by placing thermocouples in the central sample just upstream of the test area and carrying out the test: it was then verified that when the thermocouple came into the test area, i.e. in contact with the heating pads, the setpoint temperature $T_{proc}$ was reached almost immediately. Indeed, the isothermal temperature regime was quasi-instantaneously reached since carbon is a good heat conductor, but also because the pressure applied compacted the sample hence increased fiber-to-fiber contacts, the reinforcement was thin, and the low speed resulted in little contact changes over time.

## Test Procedure

The test procedure consisted first in fixing both lateral specimens to the heating pads and clamping the central specimen in the upper attachment of the traction machine crosshead. The actuator was then activated to apply the closing pressure, and the temperature of the pads increased up to the setpoint temperature $T_{proc}$. When $T_{proc}$ was reached homogeneously both on the pads and in the central sample, the test started. It began with a preload of 10 N at a crosshead displacement rate of 10 mm min⁻¹, after which the crosshead kept pulling the central specimen at the preselected displacement rate. The input parameters of the test were: the displacement rate of the crosshead representing the relative speed between tow adjacent plies, and the normal stress applied. The outputs were: displacement, tensile force, normal force, and temperature as a function of time.

## Data Post-processing

The friction coefficient $CoF$ and modified Hersey number $H^*$ were calculated as follows:

$$ CoF = \frac{F_T}{2F_N} ; \quad H^* = \frac{H}{\eta} = \frac{\gamma}{F_N} \qquad (11) $$

where $F_T$ and $F_N$ are, respectively, the tangential force and the normal force, $\eta$ is the TP veil viscosity and $\gamma$ the sliding velocity. Let us notice that the factor 2 in the expression of the friction coefficient appears (compared to Eq. 1) since it was assumed that the problem is perfectly symmetrical and consequently the contact area is twice the tested surface. We decided to consider $H^*$, i.e., the Hersey number without the influence of viscosity, since in this study isothermal conditions were fulfilled and hence viscosity was assumed to be constant in the considered normal pressure range.

# RESULTS AND DISCUSSION

## Out-of-Plane Bending

The out-of-plane bending characterization was achieved in three steps: first some tests on single plies were performed in order to validate the proposed methodology, then single plies were characterized, and finally the method was extended to UD stacks to get closer to industrial needs. For confidentiality reasons, all the moment values presented are normalized.
<figure><img src="image_210.png" /></figure>

<figure><img src="image_211.png" /><figcaption>**FIGURE 9** | Results of bending tests on single plies at room temperature using three different lumped masses $m_4$, $m_6$ and $m_7$ such that $m_7 > m_6 > m_4$: (A) experimental data of the deformed sample mean line and (B) normalized bending moment-curvature relationship using Voce's model. Each curve corresponds to one single test.</figcaption></figure>

<figure><img src="image_212.png" /></figure>

<figure><img src="image_213.png" /><figcaption>**FIGURE 10** | Results of bending tests on single plies at $T_{proc}$ (experimental data in dashed line, optimized model in solid line): (A) deformed sample mean line with $err_{95\%} < 2.5\%$, and (B) normalized bending moment-curvature relationship using Voce's model with $err_{95\%} < 3\%$. Curves 01 to 04 correspond to single tests; curve moy is the average model.</figcaption></figure>

# Results

## Preliminary Tests on Single Plies at Room Temperature

Preliminary tests were first conducted in order to validate the proposed methodology; these tests were carried out at room temperature to validate the approach in lighter conditions. In order to assess the test repeatability, 10 samples were first tested at room temperature with the same lumped mass $m_7$ and same operator. The deformed shapes are plotted in **Figure 9A** (orange solid lines) and Voce's models corresponding to moment-curvature relationships are plotted in **Figure 9B** (orange solid lines). The scattering in deformed shape led to 8% relative standard deviation on maximum bending moments, computed using Eq. 7 at clamping for each sample. Note that these values are not visible on **Figure 9B** and correspond to different

curvature values depending on the sample. This scatter was mainly due to variability on the positioning of the lumped mass at the tip of the sample. However, although rather low at room temperature, it was verified that this scatter was even reduced in temperature conditions.

This test was repeated with two other lumped masses $m_4$ and $m_6$ such that $m_7 > m_6 > m_4$ with a difference of 50% between $m_4$ and $m_7$, in order to modify the maximum loading (for confidentiality reasons, values of lumped masses are not given). The deformed shapes obtained for the three samples of the test are plotted in **Figure 9A** (purple and green dashed lines); as the mass increases, the maximum curvature of the sample (i.e. close to clamping) increased. Voce's models corresponding to moment-curvature relationships are plotted in **Figure 9B** (purple
<figure><img src="image_214.png" /><figcaption>**FIGURE 11** | Normalized bending response at $T_{proc}$ for a single ply (1p) and for 2- (2p) and 4-ply (4p) UD stacks, and theoretical response for a ply of double thickness (2p Kirch). Each experimental (exp) curve corresponds to one single test; solid lines (moy) show corresponding average optimized models.</figcaption></figure>

and green dashed lines); one can observe that the moment-curvature relationships obtained are within the envelope defined by the results with $m_7$. Therefore, for different loads, different ranges of curvatures could be described but the characterized bending response remained the same. It was thus the intrinsic response of the material that was characterized by this test. In a previous work, de Bilbao et al. (2010) illustrated the same test characteristic by changing the length of the sample.

## Results for Single Plies at $T_{proc}$

In order to identify the out-of-plane bending behavior of single plies under forming conditions, tests were conducted at $T_{proc}$ for four samples. **Figure 10A** presents the experimental deformed shapes of the four samples along with the corresponding optimized model obtained in step 3.2.3 with **Eq. 3**. The error $err_{%}^{(z)}$ was lower than 2.5% for all four samples which justified the choice of the deformed shape function (**Eq. 3**) to describe the non-linear deformed shape of the samples at $T_{proc}$ (**Figure 10A**).

In **Figure 10B**, curvatures and normalized bending moments computed using **Eqs 4** and **7**, respectively, are plotted for the four samples, as well as the corresponding optimized Voce's models (**Eq. 8**). The average model is also plotted, its parameters $R_0$, $R_{inf}$ and $\kappa_{lim}$ were taken as the average of the parameters of the four samples. This **Figure 10B** illustrates that the selected Voce's model described in a satisfactory manner the moment-curvature relationship at $T_{proc}$, with an error $err_{%}^{(M)}$ lower than 3% for every sample.

## Results on Stacks at $T_{proc}$

The bending test was extended to UD stacks of two, four, six and eight plies, with at least three samples for each. Tests on six- and 8-ply stacks could not be processed since the validity criteria were not satisfied. **Figure 11A** presents the averaged bending moment

$v_s$. curvature for various numbers of plies including a single ply, obtained by averaging Voce's model triplets of parameters. One can notice, again, that the Voce's model enabled to describe properly the moment-curvature evolution for either single plies, two- or 4-ply stacks.

## Discussion

### Comments on the Choice of the Models

Let us first observe in **Figure 8** the deformed shape and the moment-curvature relationship obtained for one sample tested at room temperature, meeting the validity criterion. The results presented in this figure exhibit a very low error between the deformed function and the mean-line obtained using image processing. This was also demonstrated with further results plotted in **Figure 10A** where $err_{%}^{(z)} < 2.5\%$ for the four samples. This justified the choice of the deformed function selected (**Eq. 3**). Indeed, such function exhibited the following characteristics: the origin of the graph was positioned at the clamping ($z(0) = 0$), the slope vanished at this same point ($z'(0) = 0$), and the function tended to a linear behavior at the free end of the sample-three features observed on experimental data.

Besides, **Figure 8B** shows that the selected Voce's model (**Eq. 8**) satisfactorily described the moment-curvature relationship obtained experimentally, with little error. Again, this tendency was more general as illustrated in **Figure 10B** with a very low error on the moment ($err_{%}^{(M)} < 3\%$ for every sample). Voce's model provided a good description of the observed behavior by the mean of a single class function $C^2$ which is non-linear for small curvatures and quasi-linear beyond a certain curvature threshold $\kappa_{lim}$. In **Eq. 8**, $R_0$ corresponds to the slope of the quasi-linear part for curvature higher than this threshold $\kappa_{lim}$ and is referred to as asymptotic tangent bending stiffness (per unit of width) thereafter; it is the minimum value of the bending stiffness. Furthermore, $R_0 + R_{inf}/\kappa_{lim}$ corresponds to the initial slope (i.e. for curvatures close to 0) and is referred to as initial tangent bending stiffness (per unit of width); it is the maximum value of the bending stiffness.

Overall, the proposed models yielded a good fit of the experimental data. In **Figures 8B,10B**, it can be observed that the major differences between the raw results and the Voce's model (Voce, 1955) are located either at the foot of the curve (low moments and low curvatures) corresponding to the free end of the sample where the moment vanishes, or at the end of the curve (high moment and high curvatures) corresponding to the clamp where the moment is maximum. Firstly, with regard to the foot of the curve, it should be noted that while the computed bending moment *exactly* equals zero at the free end of the sample, the curvature tends toward zero but only *approximately* because of the deformed shape function used (**Eq. 3**). This is not physical as it would describe a deformation in a stress-free state. Therefore we chose to impose $M_0 = 0$ in Voce's model; this was not the case in other works (de Bilbao et al., 2010; Liang et al., 2014). Secondly, the large curvatures-large moments extremity of the experimental response seems to reveal a sort of vertical tangent which would correspond to a capping of the curvature. This was actually due to the clamping which prescribed a zero slope to the sample,
whereas Voce's model intrinsically does not account for this effect on the curvature. However, this effect was concentrated and did not affect the overall response identified; at the most it could increase artificially the error err<sub>%</sub><sup>(M)</sup> calculated (Eq. 9).

Using intermediate models to build the moment-curvature relationship required to estimate the error introduced by the two successive regressions. The uncertainty on bending stiffness B due to the uncertainty on curvature κ is given by:

$$ \Delta B = \frac{R_{\text{inf}}}{\kappa_{\text{lim}}^2} e^{-\kappa/\kappa_{\text{lim}}} \Delta \kappa. \qquad (12) $$

This shows that the initial tangent stiffness (i.e. for κ = 0) is influenced by k₁ and k₂ through the curvature fitting. Conversely, the asymptotic tangent stiffness (for larger curvatures) is not affected, since ΔB → 0 when κ → +∞.

## Comments on the Curvature Range Covered

In these tests, the curvatures described by the deformed functions were between 0.0576 and 0.0010 mm⁻¹, which correspond to curvature radii from 17 to 1,000 mm approximately. Our aim was to describe curvatures representative of industrial geometries, for which curvature radii can reach 5 mm (corresponding to a curvature of 0.2 mm⁻¹). Such high curvature could not be reached experimentally with this method without the occurrence of textit{b}reaks in the deformed shape, as explained earlier. However, since Voce's model is linear beyond the threshold curvature which is much lower than 0.2 mm⁻¹ (κ<sub>lim</sub> = 0.0088 mm⁻¹ here), the moment-curvature relationship could then be extrapolated to larger curvatures. Hence, the model obtained here can be extended to curvature range of industrial geometries.

## Comments on the Identified Responses

First, a comparison between Figures 9B,10B highlights that with increasing temperature the bending moment is reduced for a given range of curvature. The softening effect of temperature is observed with a stiffnesses ratio between the two temperatures of the order of 10⁰ for the asymptotic bending stiffness and 10¹ for the initial bending stiffness.

The results obtained at T<sub>proc</sub> with the proposed procedure were cross-checked with preliminary tests using a fixed angle flexometer (Peirce, 1930) as introduced earlier. It was verified that the bending stiffness value given by Peirce model (Peirce, 1930) was within the range of bending stiffnesses (i.e. bounded by the asymptotic and initial tangent stiffnesses) obtained by the proposed method. More precisely, it was observed that Peirce stiffness is closer to the initial tangent bending stiffness, which is consistent since the fixed angle flexometer generates only small curvatures as the sample only bends under its own weight.

One limitation of this bending test was the thickness of the specimens, as it limited the number of plies of the tested stacks. Indeed, a too large thickness (or number of plies) led to the local buckling phenomenon explained earlier: an opening (fish-eye) induced by the through-thickness decohesion of the sample appeared. This phenomenon was due both to the inextensibility of individual fibers and to local singularities. In

that case, the hypothesis of continuity of the material, and consequently the interpretation of the results according to the above method, were no longer valid. It was therefore decided to minimize these effects, which consisted in reducing the sample curvature by decreasing the additional mass. However, lower curvatures were obtained and therefore the quasi-linear response could not be appropriately described as the value of the regression parameter R₀ and consequently the asymptotic tangent bending stiffness could be overestimated. The choice of the mass was thus a compromise between maintaining the cohesion of the sample, and testing over a wide range of curvatures.

With this method we characterized the bending responses of single plies as well as two- and 4-ply stacks at T<sub>proc</sub>, and an average model was derived for each configuration; they are recalled in Figure 11. This figure clearly illustrates the need for identification of plies and stacks bending response for two reasons. Firstly, predicting individual ply non-linear response is a real challenge since a simple calculation of the theoretical linear bending stiffness of a continuous (Kirchhoff) beam with longitudinal stiffness and dimensions corresponding to a HiTape® ply shows two orders of magnitude difference with the identified experimental response. Secondly, a supplementary curve may be observed on this figure: it corresponds to the theoretical bending response of a homogeneous continuous medium according to Kirchhoff theory, with the same experimental behavior than one single ply at T<sub>proc</sub> but of double-thickness. As it can be verified, this latter theoretical behavior of a double-thickness single ply appears to be four times stiffer than the experimental behavior identified for a 2-ply stack.

These differences can certainly be partly attributed to the invalidity of the continuous medium theory in this context, as there may have been some fibrous rearrangements during testing, but we also assume that they are accentuated by the presence of the thermoplastic veil at the inter-ply. Indeed, during the bending of a stack, the fiber inextensibility induces transverse shear which can only concentrate in the soft TP veil, otherwise out-of-plane buckling of plies in the inner curvature would occur (Drapier et al., 1996, 2001). In other words, this means that the overall forming capability of the stack originates from inter-ply sliding. In conclusion, either for single plies or stacks, the continuous media (Cauchy) theory can clearly not describe the response of single HiTape® plies, and the major role of the TP veil is demonstrated in the bending of a stack.

## Inter-Ply Sliding Results

Let us recall the pressure and speed ranges representative of hot vacuum forming conditions that we introduced in Section 3.3: pressures around 1 bar and relative speeds of two adjacent plies of the order of magnitude of a few millimeters in some seconds, i.e. 1–10 mm min⁻¹. In order to assess the relative influence of velocity and normal pressure in those ranges, five configurations were tested with four to five samples each. Keeping speed at 2.0 mm min⁻¹, tests were carried out at 1.35, 1.00 and 0.50 bar, respectively. Tests at bounding pressures of 1.35 and 0.50 bar were also performed at 0.5 mm min⁻¹. As this work was the first
<figure><img src="image_215.png" /><figcaption>**FIGURE 12** | Results of inter-ply sliding test: maximum and static friction coefficients vs. pressure.</figcaption></figure>

<figure><img src="image_216.png" /><figcaption>**FIGURE 13** | Results of inter-ply sliding test: maximum and static friction coefficients vs. adjacent plies relative speed.</figcaption></figure>

inter-ply sliding characterization of HiTape®, we decided to focus on 0°/0° sliding, i.e. both lateral and central specimen were aligned with loading direction.

The resulting tangent force vs. displacement curves (not presented for confidentiality reasons) exhibited the same shape than wet friction results of friction coefficient vs. displacement from literature, for example from Sachs et al. (2012) (**Figure 6B** bottom). Such comparison makes sense as friction coefficient is related to tangent force only by normal force which is kept constant during the test. After tensionning of the sample, displacement initiated (at a tangent force value $F_{stat}$). Then, as the prescribed displacement increased, the tangent force also increased to a maximum value $F_{max}$ and then decreased. $F_{stat}$ was used to compute the static friction coefficient $CoF_{stat}$ using Eq. 11

<figure><img src="image_217.png" /><figcaption>**FIGURE 14** | Stribeck curve for HiTape® over the pressure and speed ranges representative of double-membrane hot vacuum forming conditions.</figcaption></figure>

with $F_T = F_{stat}$. Similarly $F_{max}$ was also used to compute a so-called maximum friction coefficient $CoF_{max}$ using Eq. 11 with $F_T = F_{max}$. Let us highlight that as the normal force $F_N$ was kept constant during the test, the friction coefficient vs. displacement curve also exhibited the same shape than **Figure 6B** (bottom).

The results of static and maximum coefficients of friction are presented with their standard deviations as a function of pressure in **Figure 12** and as a function of the relative speed in **Figure 13**. A strong influence of pressure and velocity was observed in the ranges of pressure and velocity studied. Note that the standard deviations on every configuration were small (less than 10%) compared to the differences observed from one configuration to another.

## Discussion

The shape of the tangent force-displacement response obtained was consistent with literature results presented for composites pre-impregnated with a TP resin tested in a molten state, as illustrated in **Figure 6B** taken from Sachs et al. (2012). However, unlike in Ten Thije et al. (2011) and Fetzatsidis et al. (2013), the tangent force measured in our experiment did not clearly stabilize. Consequently, we could not extract from the presented tests a dynamic friction coefficient. This phenomenon is attributed to a fibrous rearrangement during the test and to possible penetration of the veil inside the UD fiber bed.

When forming a stack, the expected relative displacement between two adjacent plies is of the range of a few millimeters. Moreover, purely geometrical considerations show that forming two 0.2 mm thick inextensible layers on a quadrant implies a relative displacement between of 0.3 mm the two layers. This highlights that only the very first millimeters of sliding of the test are of interest for us. Using bounding values of $F_{stat}$ and $F_{max}$ for the tangent force is therefore justified, and we can state that friction coefficient during forming is bounded by $CoF_{stat}$ and $CoF_{max}$.
In order to conclude on the type of lubricated regime that characterizes the inter-ply response during hot forming, both coefficients of friction obtained CoF<sub>stat</sub> and CoF<sub>max</sub> are plotted as functions of the modified Hersey number (Eq. 11) in Figure 14. In both cases, an increasing linear relationship between the friction coefficients and the modified Hersey number *H*<sup>*</sup> may be observed. Recalling the Stribeck curve (Figure 6), we can therefore conclude that a hydrodynamic lubricated regime prevailed, which corresponds to the presence of a lubricating film of thickness larger than the roughness of the reinforcement. In terms of industrial forming, this means that inter-ply mobility is promoted at low Hersey numbers, i.e. at high forming pressures and low velocities, and may not depend of the relative orientations of the plies.

## CONCLUSION

In this work, the deformation mechanisms of HiTape® reinforcements were inventoried: intra-ply mechanisms control the response of a single ply, and inter-ply mechanisms (especially inter-ply sliding) are involved when considering the forming of stacks. More precisely, HiTape® mainly deforms in bending (and forms wrinkles) which is controled by the fiber response as well as the presence of the TP veil that acts between plies as a lubricating layer providing shear compliance especially at the forming process temperature (which is above the veil melting temperature). This study was therefore dedicated to the characterization of two deformation modes encountered by HiTape® reinforcements during double-membrane hot vacuum forming process and for which the behavior was unknown, namely out-of-plane bending and inter-ply friction. Dedicated apparatuses were designed on purpose, in order to meet the requirements of handling, heating, and finely controling the characterization of non-linear bending response and inter-ply sliding in conditions representative of the industrial double-membrane hot vacuum forming processes.

For bending, a modified Peirce's flexometer was designed and validated. It enabled to assess the complete isothermal moment-curvature non-linear relationship in a single test, taking advantage of the curvature range covered by the deformation of single plies and stacks under a terminal loading that could be adjusted to reach the curvatures in scope. The bending responses

identified correspond to fibrous media bending for which the internal mobility of the fibers leads to a range of stiffnesses much lower that what can be expected in continuous Cauchy's media with the same elongational modulus. Also, the bending stiffnesses identified for stacks clearly showed that the thermoplastic veil between plies has a strong influence and controls to a high extent the stack bending overall response.

Such results justified the development of a second apparatus designed for the characterization of inter-ply friction: a pull-through test with controled temperature and pressure. A hydrodynamic lubricated friction regime was observed through a linearly increasing relationship between the friction coefficients–both static and maximum–and the modified Hersey number (divided by the assumed constant viscosity of the thermoplastic veil). From an industrial point of view, this result showed that in order to limit the occurrence of wrinkle-type defects during forming, inter-ply sliding should be promoted with forming conditions inducing high pressure and low velocity. This guideline should be verified in further work with dedicated forming experiments.

## DATA AVAILABILITY STATEMENT

The datasets presented in this article are not readily available because they are confidential. Requests to access the datasets should be directed to laure.bouquerel@hexcel.com.

## AUTHOR CONTRIBUTIONS

LB: main co-author, the material presented herein comes out from her PhD work. NM: LB's co-supervisor and co-contributor to this article. SD: LB's main supervisor and co-contributor to this article.

## FUNDING

This work produced in LB's PhD was supported both by ANRT, the French National Agency for Research and Technology, under the grant 2016/0094, and by the company Hexcel Reinforcements SASU.

## REFERENCES

* Abbott, G. M., Grosberg, P., and Leaf, G. A. V. (1971). The mechanical properties of woven fabrics Part VII : the hysteresis during bending of woven fabrics. *Textil. Res. J.* 41, 345–358. doi:10.1177/004051757104100411

* Alshahrani, H., and Hojjati, M. (2017). A new test method for the characterization of the bending behavior of textile prepregs. *Composites Part A: Applied Science and Manufacturing*. 97, 128–140. doi:10.1016/j.compositesa.2017.02.027

* Angel, B., and Graef, J. (2016). Study of the bending stiffness of fibre-reinforced thermoplastics at forming temperature. *JEC Composites Magazine.*, 103, 76–79.

* Bel, S., Hamila, N., Boisse, P., and Dumont, F. (2012). Finite element model for NCF composite reinforcement preforming: composites Part A: applied Science and manufacturing: importance of inter-ply sliding. *Compos. Part A: Appl. Sci. Manuf.* 43 (12), 2269–2277. doi:10.1016/j.compositesa.2012.08.005

* Blais, M. (2016). Modélisation et suivi du procédé par infusion de résine sur une nouvelle génération de renforts structuraux pour l'aéronautique. PhD. thesis, Saint-Étienne (France): Ecole Nationale Supérieure des Mines de Saint-Étienne.

* Blais, M., Moulin, N., Liotier, P. J., and Drapier, S. (2017). Resin infusion-based processes simulation : coupled Stokes-Darcy flows in orthotropic preforms undergoing finite strain. *Int. J. Material Form.* 10 (1), 43–54. doi:10.1007/s12289-015-1259-2

* Boisse, P., Hamila, N., Vidal-Sallé, E., and Dumont, F. (2011). Simulation of wrinkling during textile composite reinforcement forming. Influence of tensile, in-plane shear and bending stiffnesses. *Compos. Sci. Technol.* 71, 683–692. doi:10.1016/j.compscitech.2011.01.011
Boisse, P. (2004). *Mise en forme des renforts fibreux de composites*. Techniques de l'ingénieur Plasturgie: procédés spécifiques aux composites Ref: TIB474DUO. Article ref: am3734

Bouquerel, L., Moulin, N., Drapier, S., Boisse, P., and Beraud, J.-M. (2017). Modelling and simulating the forming of new dry automated lay-up reinforcements for primary structures. *AIP conference proceedings*. 1896 (1), 030008. doi:10.1063/1.5007995

Cao, J., Akkerman, R., Boisse, P., Chen, J., Cheng, H. S., de Graaf, E. F., et al. (2008). Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results. *Compos. Appl. Sci. Manuf.* 39, 1037–1053. doi:10.1016/j.compositesa.2008.02.016

Celle, P., Drapier, S., and Bergheau, J. M. (2008). Numerical aspects of fluid infusion inside a compressible porous medium undergoing large strains. *European Journal of Computational Mechanics*. 18, 819–827. doi:10.3166/ remm.17.819-827

Charmetant, A. (2011). Approches hyperélastiques pour la modélisation du comportement mécanique de préformes tissées de composites. PhD thesis, Villeurbanne (France): Institut National des Sciences Appliquées de Lyon

Clapp, T. G., Peng, H., Ghosh, T. K., and Eischen, J. W. (1990). Indirect measurement of the moment-curvature relationship for fabrics. *Textil. Res. J.* 60, 525–533. doi:10.1177/004051759006000906

Cornelissen, B., Sachs, U., Rietman, B., and Akkerman, R. (2014). Dry friction characterisation of carbon fibre tow and satin weave fabric for composite applications. *Composites Part A: Applied Science and Manufacturing*. 56, 127–135. doi:10.1016/j.compositesa.2013.10.006

Creech, G., and Pickett, A. K. (2006). Meso-modelling of non-crimp fabric composites for coupled drape and failure analysis. *J. Mater. Sci.* 41 (20), 6725–6736. doi:10.1007/s10853-006-0213-6

Das, A., Kothari, V. K., and Vandana, N. (2005). A study on frictional characteristics of woven fabrics. *Auten Res. J.* 5 (3), 133–140

de Bilbao, E., Soulat, D., Hivet, G., and Gasser, A. (2010). Experimental study of bending behaviour of reinforcements. *Exp. Mech.* 50, 333–351. doi:10.1007/ s11340-009-9234-9

de Bilbao, E., Soulat, D., Hivet, G., Launay, J., and Gasser, A. (2008). Bending test of composite reinforcements. *Int. J. Material Form.* 1, 835–838. doi:10.1007/ s12289-008-0265-z

Dörr, D., Brymerski, W., Ropers, S., Leutz, D. M., Joppich, T., Kärger, L., et al. (2017). A benchmark study of finite element codes for forming simulation of thermoplastic UD-Tapes *Procedia CIRP*. 66. 101–106. doi:10.1016/j.procir. 2017.03.223

Drapier, S., Grandidier, J.-C., Gardin, C., and Potier-Ferry, M. (1996). Structure effect and microbuckling. *Composites Science and Technology*. 56 (7), 861–867. doi:10.1016/0266-3538(96)00033-4

Drapier, S., Grandidier, J. C., and Potier-Ferry, M. (2001). A structural approach of plastic microbuckling in long fibre composites: comparison with theoretical and experimental results. *Int. J. Solid Struct.* 38, 3877–3904. doi:10.1016/S0020- 7683(00)00247-X

Drapier, S., and Wisnom, M. R. (1999). Finite-element investigation of the compressive strength of non-crimp-fabric-based composites. *Compos. Sci. Technol.* 59, 1287–1297. doi:10.1016/S0266-3538(98)00165-1

Dufort, L., Drapier, S., and Grédiac, M. (2001). The cross section warping in short beams under three point bending: an analytical study. *Compos. Struct.* 52, 1233–1246.

Fagiano, C. (2010). Computational modeling of tow-placed composite laminates with fabrication features. PhD thesis. Delft (Netherlands): Technische Universiteit Delf.

Fetfatsidis, K. A., Jauffrès, D., Sherwood, J. A., and Chen, J. (2013). Characterization of the tool/fabric and fabric/fabric friction for woven-fabric composites during the thermostamping process. *Int. J. Material Form.* 6, 209–221. doi:10.1007/s12289-011-1072-5

Frène, J., and Zaidi, H. (2011). *Introduction à la tribologie*. Techniques de l'ingénieur : frottement et usure Ref: TIB464DUO. Article ref: TRI 100-1

Gorczyca, J. L., Sherwood, J. A., and Chen, J. (2007). A friction model for thermostamping commingled glass-polypropylene woven fabrics. *Composites Part A: Applied Science and Manufacturing*. 38 (2), 393–406. doi:10.1016/j. compositesa.2006.03.006

Grosberg, P. (1966). The mechanical properties of woven fabrics Part II: the bending of woven fabric. *Textil. Res. J.* 36, 205–211. doi:10.1177/ 004051756603600301

Guzman-Maldonado, E., Wang, P., Hamila, N., and Boisse, P. (2019). Experimental and numerical analysis of wrinkling during forming of multi-layered textile composites. *Compos. Struct.* 208, 213–223. doi:10.1016/j.compstruct.2018.10.018

Haanappel, S. P., ten Thije, R. H. W., Sachs, U., Rietman, B., and Akkerman, R. (2014). Formability analyses of uni-directional and textile reinforced thermoplastics. *Composites Part A: Applied Science and Manufacturing*. 56, 80–92. doi:10.1016/j.compositesa.2013.09.009

Hersey, M. D. (1914). The laws of lubrication of horizontal journal bearings. *J. Wash. Acad. Sci.* 4, 542–552. doi:10.2307/24520857

Hexcel (2015a). Direct processes technology manual. *Tech. rep.* Available at: https://www.hexcel.com/Resources/Technology-Manuals.

Hexcel (2015b). HiTape® A new efficient composite technology for primary aircraft structures. *SAMPE J.* 51, 7–15.

Hivet, G., Allaoui, S., Cam, B. T., Ouagne, P., and Soulat, D. (2012). Design and potentiality of an apparatus for measuring yarn/yarn and fabric/fabric friction. *Exp. Mech.* 52, 1123–1136. doi:10.1007/s11340-011-9566-0

Howell, H. G., and Mazur, J. (1953). Amontons' law and fibre friction. *Journal of the Textile Institute Transactions*. 44, T59–T69. doi:10.1080/ 19447025308659728

Huang, N. C. (1979). Finite biaxial extension of completely set plain woven fabrics. *J. Appl. Mech.* 46, 651. doi:10.1115/1.3424621

Kang, T. J., and Yu, W. R. (1995). Drape simulation of woven fabric by using the Finite-element method. *J. Textil. Inst.* 86, 635–648. doi:10.1080/ 00405009508659040

Kawabata, S. (1980). *The standardization and analysis of hand evaluation*. 2nd edn. (Osaka, Japan: Hand Evaluation Standardization Committee).

Larberg, Y. R., Akermo, M., and Norrby, M. (2012). On the in-plane deformability of cross-plied unidirectional prepreg. *J. Compos. Mater.* 46, 929–939. doi:10. 1177/0021998311412988

Larberg, Y. R., and Akermo, M. (2011). On the interply friction of different generations of carbon/epoxy prepreg systems. *Composites Part A: Applied Science and Manufacturing*. 42, 1067–1074. doi:10.1016/j.compositesa.2011. 04.010

Leutz, D. M. (2015). *Forming simulation of AFP material layups: material characterization, simulation and validation*. Munich, Germany: Technische Universität München.

Liang, B. (2016). Experimental and numerical study of the bending behaviour of textile reinforcements and thermoplastic prepregs. PhD thesis, Villeurbanne, France: Institut National des Sciences Appliquées de Lyon.

Liang, B., Hamila, N., Peillon, M., and Boisse, P. (2014). Analysis of thermoplastic prepreg bending stiffness during manufacturing and of its influence on wrinkling simulations. *Compos. Part A Appl. Sci. Manuf.* 67, 111–122. doi:10.1016/j.compositesa.2014.08.020

Lightfoot, J. S., Wisnom, M. R., and Potter, K. (2013). A new mechanism for the formation of ply wrinkles due to shear between plies. *Compos. Part A Appl. Sci. Manuf.* 49, 139–147. doi:10.1016/j.compositesa.2013.03.002

Lomov, S. V., Verpoest, I., Barburski, M., and Laperre, J. (2003). Carbon composites based on multiaxial multiply stitched prepforms. Part 2. KES-F characterisation of the deformability of the preforms at low loads. *Compos. Part A Appl. Sci. Manuf.* 34, 359–370. doi:10.1016/S1359-835X(03)00025-3

Madeo, A., Ferretti, M., Dell'Isola, F., and Boisse, P. (2015). Thick fibrous composite reinforcements behave as special second-gradient materials: three-point bending of 3D interlocks. *Z. Angew. Math. Phys.* 66, 2041–2060. doi:10.1007/s00033-015-0496-z

Margossian, A., Bel, S., and Hinterhoelzl, R. (2015). Bending characterisation of a molten unidirectional carbon fibre reinforced thermoplastic composite using a Dynamic Mechanical Analysis system. *Compos. Part A Appl. Sci. Manuf.* 77, 154–163. doi:10.1016/j.compositesa.2015.06.015

Margossian, A., Bel, S., and Hinterhoelzl, R. (2016). On the characterisation of transverse tensile properties of molten unidirectional thermoplastic composite tapes for thermoforming simulations. *Compos. Part A Appl. Sci. Manuf.* 88, 48–58. doi:10.1016/j.compositesa.2016.05.019

Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. *J. Soc. Ind. Appl. Math.* 11, 431–441. doi:10.1137/0111030
McGuinness, G. B., and Bradaigh, C. O. (1998). Characterisation of thermoplastic composite melts in rhombus-shear : the picture-frame experiment. *Compos. Part A Appl. Sci. Manuf.* 29, 115-132. doi:10.1016/S1359-835X(97)00061-4

Montero, L., Allaoui, S., and Hivet, G. (2017). Characterisation of the mesoscopic and macroscopic friction behaviours of glass plain weave reinforcement. *Compos. Appl. Sci. Manuf.* 95, 257-266. doi:10.1016/j.compositesa.2017.01.022

Najjar, W., Pupin, C., Legrand, X., Boude, S., Soulat, D., and Dal Santo, P. (2014). Analysis of frictional behaviour of carbon dry woven reinforcement. *J. Reinforc. Plast. Compos.* 33, 1037-1047. doi:10.1177/0731684414521670

Ngo Ngoc, C., Bruniaux, P., and Castelain, J. M. (2002). "Modeling friction for yarn/fabric simulation Application to bending hysteresis," in *14th European Simulation Symposium*. Editor V. Krug (Dresden, Germany: SCS Europe BVBA).

Orawattanasrikul, S. (2006). *Experimentelle Analyse des Scherdeformation biaxial verstärkter Mehlagengesticke*. PhD thesis. Dresden, Germany: Technische Universität Dresden.

Pannier, Y. (2006). Identification de paramètres élastoplastiques par des essais statiquement indéterminés : mise en oeuvre expérimentale et validation de la méthode des champs virtuels. PhD thesis, Paris (France): Arts et Métiers ParisTech.

Peirce, F. T. (1930). The "handle" of cloth as a measurable quantity. *Journal of the Textile Institute Transactions*. 21, T377-T416. doi:10.1080/19447023008661529

Pickett, A. K. (2018). *Process and mechanical modelling of engineering composites*. Stuttgart, Germany: IFB University of Stuttgart.

Rashidi, A., Montazerian, H., Yesilcimen, K., and Milani, A. (2020). Experimental characterization of the inter-ply shear behavior of dry and prepreg woven fabrics: significance of mixed lubrication mode during thermoset composites processing. *Compos. Appl. Sci. Manuf.* 129, 105725. doi:10.1016/j.compositesa.2019.105725

Rosen, B. (1964). Mechanics of composite strengthening. in *Seminar of the American society for metals*. Cleveland, OH: Metal Parks, 37-75.

Sachs, U., Akkerman, R., and Haanappel, S. P. (2014). Bending characterization of UD composites. *Key Eng. Mater.* 611-612, 399-406. doi:10.4028/www.scientific.net/KEM.611-612.399

Sachs, U., Fetfatsidis, K. A., Schumacher, J., Ziegmann, G., Allaoui, S., Hivet, G., et al. (2012). A friction-test benchmark with Twintex PP. *Key Eng. Mater.* 504-506, 307-312. doi:10.4028/www.scientific.net/KEM.504-506.307

Senner, T., Kreissl, S., Merklein, M., Meinhardt, M., and Lipp, A. (2014). Bending of unidirectional non-crimp-fabrics: experimental characterization, constitutive modeling and application in finite element simulation. *J. Inst. Eng. Prod.* 9, 1-10. doi:10.1007/s11740-014-0568-5

Soteropoulos, D., Fetfatsidis, K. A., Sherwood, J. A., and Langworthy, J. (2011). Digital method of analyzing the bending stiffness of non-crimp fabrics. AIP Conference Proceedings, 1353, 913-917. doi:10.1063/1.3589632

Stribeck, R. (1903). *Die wesentlichen Eigenschaften der Gleit- und Rollenlager*. Berlin, Germany: Springer.

Ten Thije, R. H. W., Akkerman, R., Ubbink, M., and Van Der Meer, L. (2011). A lubrication approach to friction in thermoplastic composites forming processes. *Composites Part A: Applied Science and Manufacturing*. 42, 950-960. doi:10.1016/j.compositesa.2011.03.023

Trejo, E. A., Ghazimoradi, M., Butcher, C., and Montesano, J. (2020). Assessing strain fields in unbalanced unidirectional non-crimp fabrics. *Compos. Appl. Sci. Manuf.* 130, 105758. doi:10.1016/j.compositesa.2019.105758

Vanclooster, K., Lomov, S. V., and Verpoest, I. (2010). Simulation of multi-layered composites forming. *Int. J. Material Form.* 3, 695-698. doi:10.1007/s12289-010-0865-2

Voce, E. (1955). A practical strain-hardening function. *Metallurgica*. 51, 219-226.

Wang, P., Hamila, N., and Boisse, P. (2013). Thermoforming simulation of multilayer composites with continuous fibres and thermoplastic matrix. *Compos. B Eng.* 52, 127-136. doi:10.1016/j.compositesb.2013.03.045

Xiong, H., Hamila, N., and Boisse, P. (2019). Consolidation modeling during thermoforming of thermoplastic composite prepregs. *Materials*. 12, 2853. doi:10.3390/ma12182853

Yuksekkay, M. E. (2009). More about fibre friction and its measurements. *Textil. Prog.* 41 (3), 141-193. doi:10.1080/00405160903178591

Zhu, B., Yu, T. X., Zhang, H., and Tao, X. M. (2011). Experimental investigation of formability of commingled woven composite preform in stamping operation. *Compos. B Eng.* 42, 289-295. doi:10.1016/j.compositesb.2010.05.006

**Conflict of Interest:** LB was employed by the company Hexcel Reinforcements SASU.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 Bouquerel, Moulin and Drapier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_218.png" /></figure>

Check for
updates

# An Incremental Cohesive Law for Delamination Under a Mixed Mode Loading

Man Zhu<sup>1,2</sup>, Larissa Gorbatikh<sup>2*</sup> and Stepan V. Lomov<sup>2,3</sup>

<sup>1</sup>State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing, China, <sup>2</sup>Department of Materials Engineering, KU Leuven, Leuven, Belgium, <sup>3</sup>Center for Design, Manufacturing and Materials, Skolkovo Institute of Science and Technology, Moscow, Russia

Cohesive zone models rely on the formulation of a cohesive constitutive law. The latter describes the relation between displacement and traction in a cohesive element at an integration point. Cohesive constitutive laws in the presence of opening and shearing modes are less studied in comparison with those formulated for a single mode, particularly when the mode mixity changes. The mode mixity at an integration point is determined by the load history at the point. In this study, a formulation of the cohesive constitutive law is proposed for a mixed mode loading condition with the ability to deal with the variation in mode mixity. The proposed law is constructed incrementally and takes into account the load history. The validation is performed by simulating delamination in carbon fiber/epoxy composites in the mixed-mode bending test that is commonly used to characterize the inter-laminar fracture toughness. Although the mode mixity is fixed in this test at the specimen level, it varies locally at the element level. Cohesive constitutive laws proposed in the literature predict macroscopic delamination behavior that is dependent on the strength of the interface, while, according to the analysis of linear elastic fracture mechanics, the dependence is expected to be only on the fracture toughness. Predictions with the current formulation, where the cohesive law is updated incrementally, show low sensitivity to the interface strength. The structural response simulated with it had a good agreement with the analytical solution of linear elastic fracture mechanics.

**Keywords:** cohesive law, delamination, mixed mode, composite, finite element method

## OPEN ACCESS

**Edited by:**

Patricia Krawczak,
IMT Lille Douai, France

**Reviewed by:**

Frédéric Lachuad,
Institut Supérieur de l'Aéronautique et
de l'Espace (ISAE-SUPAERO), France
Paolo S. Valvo,
University of Pisa, Italy

**\*Correspondence:**

Larissa Gorbatikh
larissa.gorbatikh@kuleuven.be

**Specialty section:**

This article was submitted to
Mechanics of Materials,
a section of the journal
Frontiers in Materials

**Received:** 15 June 2020

**Accepted:** 02 October 2020

**Published:** 12 November 2020

**Citation:**

Zhu M, Gorbatikh L and Lomov SV (2020) An Incremental Cohesive Law for Delamination Under a Mixed Mode Loading. Front. Mater. 7:572995.

doi: 10.3389/fmats.2020.572995

## INTRODUCTION

Laminated materials can experience de-cohesion or delamination under different loading conditions, from simple single mode to more complex mixed mode conditions. The mixed mode conditions, where both normal and shear tractions act on the inter-laminar interface, are, for example, present in structures that are curved, like flanges of the ribs in aircraft wings (Gözlüklü and Coker, 2012; Gozluklu et al., 2015). The mixed mode can also develop in plane structures subjected to simple loading conditions after the onset of delamination, like in the case of laminated plates subjected to compression after impact (González et al., 2012). Modeling tools facilitate design of reliable structures.

Cohesive zone models (CZMs) have successfully proven their efficiency in modeling interface behavior in laminated materials including fibre-reinforced composites (Alfano and Crisfield, 2001; Borg et al., 2002; Turon et al., 2006; Jokić et al., 2016; Jensen et al., 2019; Liu et al., 2019). The concept
of CZMs was originally proposed by Barenblatt (Barenblatt, 1959; Barenblatt, 1962) and Dugdale (Dugdale, 1960) as an approach to depicting crack formation and development. In CZMs, cohesive elements are placed in the potential crack path, and the propagation of a crack is described by the behavior of the integration points in cohesive elements. In CZMs, the damage status of the integration points develops in a progressive way controlled by a cohesive constitutive law. This cohesive constitutive law is formulated as a relationship between the traction and the displacement jump caused by normal and shear tractions acting at the interface. An accurate cohesive constitutive law is essential for simulating delamination behavior correctly. Cohesive laws for single modes (i.e., pure mode I and pure mode II) have been studied extensively with many formulations available. Modeling of delamination under a mixed-mode condition is still an active research subject. In some studies it is proposed to treat the mixed mode as a special case of single modes and to formulate their cohesive laws in the way it is done for single modes. An additional precondition is often assumed that the mode mixity should be fixed during the whole loading history. However the fixed mode mixity is hardly ensured in practical situations. Even if the mode mixity is considered to be fixed at the macro-scale, it may vary locally, like in the mixed-mode bending (MMB) test which is designed to characterize inter-laminar fracture toughness of the composite at a fixed mode mixity at macro-scale (ASTM D6671, 2013). Moreover, whether the constant mode mixity at the macro-scale can be ensured in MMB test has been questioned in the literature (Bennati et al., 2013a; Bennati et al., 2013b). Detailed simulations of the MMB test have also shown that the mode mixity varies at different integration points of cohesive elements and that even at the same integration point the mode mixity changes as the test progresses (Turon et al., 2010; Harper et al., 2012; De Moura et al., 2016). Therefore, the question rises “How to include variations in mode mixity into the cohesive constitutive law?”.

In the state-of-the-art simulations (Turon et al., 2006; Turon et al., 2010; Turon et al., 2018; Jensen et al., 2019), a popular way to consider the variable mode mixity is based on a set of pre-defined values of mode mixities. Firstly, different cohesive constitutive laws are defined for these fixed mode mixities. Then, during the loading history each integration point is

checked for its mode mixity and once the latter is changed the cohesive law acting there is revised to account for the new mode mixity. By using this methodology, a strong assumption is made, namely that the damage evolution at the integration point under variable mode mixity is the same as the one with a fixed mode mixity. This means the effect of the loading history is not considered.

Another issue with the mixed-mode simulation by using existing cohesive constitutive laws is that the force-displacement relation in the MMB test at macro-scale depends on the inter-laminar strength of cohesive elements. This conflicts with the linear elastic fracture mechanics (LEFM) analysis according to which delamination is only controlled by the fracture toughness of the interface. Indeed, let us take a look at an analytical expression of the force-displacement relation in the MMB test (as shown in Figure 1) (ASTM D6671, 2013):

$$U_{\text{lever}} = \frac{P_{\text{lever}}}{8bE_{11}h^3L^2} \left[ 4(3c - L)^2 (a + h\chi)^3 + (c + L)^2(2L^3 + 3(a + 0.42h\chi)^3) \right] \quad (1)$$

$$P_{\text{lever}} = \sqrt{\frac{16G_c b^2 E_{11} h^3 L^2}{12(3c - L)^2 (a + h\chi)^2 + 9(c + L)^2 (a + 0.42h\chi)^2}} \quad (2)$$

$$\chi = \sqrt{\frac{E_{11}}{11G_{13}} \left\{ 3 - 2 \left( \frac{\Gamma}{1 + \Gamma} \right)^2 \right\}}, \quad \Gamma = 1.18 \frac{\sqrt{E_{11}E_{22}}}{G_{13}} \quad (3)$$

where $U_{\text{lever}}$ and $P_{\text{lever}}$ are the applied displacement on MMB specimen and corresponding reaction force; $c$ and $L$ are the geometry parameters of the test set-up; $b$ and $h$ are the width and thickness of the specimen; $a$ is the crack length; $E_{11}$, $E_{22}$, and $G_{13}$ are the longitudinal Young's modulus, transverse Young's modulus and shearing modulus of the composite laminates, respectively; $G_c$ is the fracture toughness of the inter-laminar interface. From Eqs 1–3, it can be seen that once the test set-up and specimen geometry are fixed, the structural response of the specimen is only determined by the inter-laminar fracture toughness and has no dependency on the inter-laminar strength. However, this is not the case when the MMB test was simulated with the existing mixed-mode cohesive constitutive law (Turon et al., 2010; Harper et al., 2012). More specifically, different load-displacement relations were obtained

<figure><img src="image_219.png" /><figcaption>FIGURE 1 | Apparatus of the MMB test.</figcaption></figure>
for different values of the inter-laminar strength in cohesive elements. The application of this model relies on the selection of a specific value of the inter-laminar strength.

To get rid of this inter-laminar strength dependency, Turon et al. revised their original cohesive constitutive law twice (Turon et al., 2010; Turon et al., 2018), which is one of the mostly used CZMs for simulating delamination in composite structures. In their first update (Turon et al., 2010), a relationship between the pure mode I and mode II inter-laminar strength was given:

$$ \frac{\sigma_{20}}{\sigma_{10}} = \sqrt{\frac{G_{c\_II}}{G_{c\_I}}} \qquad (4) $$

where $\sigma_{10}$ and $\sigma_{20}$ are pure mode I and II inter-laminar strength, respectively; $G_{c\_I}$ and $G_{c\_II}$ are inter-laminar fracture toughness under pure mode I and II loading conditions, respectively. Only when this relationship is satisfied, the dependence of MMB test's structural response on the inter-laminar strength can be eliminated. However, the mandatory requirement of the relationship between pure mode I and II inter-laminar fracture toughness and strength values makes no sense from a physical point of view. This relationship was then abandoned in the second update by Turon et al. (2018). Instead, a mode-dependent stiffness was adopted as a function of mode mixity, namely

$$ K_m = (1 - B)K_I + BK_{II} \qquad (5) $$

where $K_I$, $K_{II}$ and $K_m$ are the stiffness used in the cohesive law for pure mode I, mode II and mixed-mode with a mode mixity of B, respectively (see details of the cohesive law definition in the next section). With this modified formulation of the cohesive constitutive law, the dependence of the macro-scale response on the inter-laminar strength can be eliminated but only during the delamination propagation after the peak load. The peak load is still affected by the inter-laminar strength, which is against the conclusion of the study of Lu et al. (2019). Lu et al. systematically analyzed the effect of inter-laminar strength on the delamination behavior, and it was found that in the presence of a stress concentration, such as a pre-delamination, the structural response is strength independent including both peak load and the response in the propagation part.

In this paper, a new formulation of a cohesive constitutive law is proposed for simulating delamination with variable mode mixity. The constitutive law is constructed incrementally based on the deformation, traction and mode mixity at the integration point obtained from the previous simulation step. In this way the effect of loading history is included in the cohesive law. For simplicity, the formulation of the proposed cohesive constitutive law is firstly explained for a pure mode condition in *Bi-Linear Cohesive Constitutive Law in the Case of a Single Mode*, and then extended to the variable mode mixity cases in *Cohesive Constitutive Law for Variable Mode Mixity*. The developed cohesive constitutive law is validated by simulating MMB test, and the results are shown in *Results of FEM Validation and Discussion*. Discussions about the mixed-mode constitutive law is also presented in the same section. Finally, *Conclusions and Outlook* summarizes conclusions and highlights of this study.

## BI-LINEAR COHESIVE CONSTITUTIVE LAW IN THE CASE OF A SINGLE MODE

A classical bilinear cohesive law is adopted in this study. It has been utilized many times in the literature both for the simulation of static and fatigue delamination due to its simplicity in implementation (Alfano and Crisfield, 2001; Turon et al., 2006; Kawashita and Hallett, 2012; Nojavan et al., 2016). The bilinear cohesive law developed by Turon (Turon et al., 2006) is widely used for modeling delamination in composite structures, and serves as basis in this study.

In the three-dimensional formulation, the displacement jump $\delta_i$ and traction $\sigma_i$ of the integration point in cohesive element have three components, one indicating opening deformation and the other two describing shearing deformation. However, in the simulation of a pure mode, only one component of the displacement jump is non-zero and the cohesive law can be simplified to the relation between two scalars: displacement jump and traction for the relevant mode. The bilinear cohesive relation between the jump and the traction consists of two sections: an initial linearly elastic response, which is followed by linear softening. **Figure 2A** gives a sketch of the bilinear cohesive law. The integration point behaves linearly-elastically with the initial stiffness, $K_0$, until the initial displacement $\delta_0$ is

<figure><img src="image_220.png" /><figcaption>FIGURE 2 | Sketch of the bilinear cohesive law: (A) damage variable $\delta_0$; (B) incremental type damage variable $c^{(t+\Delta t)}$.</figcaption></figure>
reached, after which the interface traction $\sigma$ decreases as the displacement jump increases. The initial displacement, $\delta_0$, is the displacement jump corresponding to the inter-laminar strength of the integration points in cohesive element $\sigma_0$: $\delta_0 = \sigma_0/K_0$. A damage variable, $d_0$, is defined as the degradation of the initial stiffness, $K = (1 - d_0) K_0$. The integration point becomes traction free, when displacement jump reaches the critical value $\delta_f$.

The traction at the integration point $\sigma_{coh}$ can be determined as a function of the displacement jump, $\delta_{coh}$, and damage variable, $d_0$:

$$ \sigma_{coh} = (1 - d_0)K_0\delta_{coh} \quad (6) $$

The damage variable, $d_0$, is equal to zero in the elastic region OA, and increases from zero at point A to unity at point C in the softening region AC:

$$ d_0 = \begin{cases} 0, & 0 \le \delta_{coh} < \delta_0, \\ \frac{\delta_f (\delta_{coh} - \delta_0)}{\delta_{coh} (\delta_f - \delta_0)}, & \delta_0 \le \delta_{coh} < \delta_f, \\ 1, & \delta_f \le \delta_{coh}. \end{cases} \quad (7) $$

The area under the traction-displacement jump curve (i.e., area OAC in Figure 2A) equals $G_c$, which is the critical energy release rate of the integration point. $G_c$ is also interpreted as the inter-laminar fracture toughness. The displacement at which the element becomes traction free (sometimes referred as failure displacement) can be calculated as follows:

$$ \delta_f = \frac{2G_c}{K_0\delta_0} \quad (8) $$

The strain energy release rate (SERR), which has been supplied to the cohesive element up to the moment when the displacement jump reaches the value of $\delta_{coh}$, noted here as $G_{sup}$, is equal to the area under the cohesive law, as shown by area OABD in Figure 2A:

$$ G_{sup} = \int_{0}^{\delta_{coh}} \sigma(\delta) d\delta \quad (9) $$

The energy dissipated due to accumulated damage, $G_{dis}$, is represented by area OAB in Figure 2A. It is the amount of energy that cannot be recovered after unloading: $G_{dis} = G_{sup} - (1/2)K\delta_{coh}^2$, and a new traction free surface (in other words, a crack) will be created once it is equal to the fracture toughness, $G_c$.

Based on the traditional bilinear constitutive law, we propose an incremental type of damage variable $d^{(t+\Delta t)}$, which describes reduction of stiffness $K^{(t)}$ to $K^{(t+\Delta t)}$ (Figure 2B). In the following analysis, the superscripts (t) and (t+\Delta t) indicate the variables that are obtained at time t and t+\Delta t, respectively. The definition of the incremental type of damage variable allows the cohesive constitutive law to be updated in an incremental way which will be beneficial to the simulation of the case with a variable mode mixity. At simulation time t, the residual stiffness of the

cohesive element is $K^{(t)}$, as simulation time (deformation) progresses the integration point behaves following the constitutive curve O-B-D, as shown in Figure 2B. No reduction of $K^{(t)}$ happens until the displacement reaches $\delta_0^{(t+\Delta t)}$ at time $t+\Delta t$, followed by which the incremental damage variable, $d^{(t+\Delta t)}$, increases from zero to one at $\delta_f^{(t+\Delta t)}$. The area under the constitutive curve O-B-D, is the residual SERR to be dissipated after time t, and its value is the difference between the fracture toughness and dissipated SERR: $G_{res}^{(t+\Delta t)} = G_c - G_{dis}^{(t)}$. It is the total energy to be dissipated from time t till the integration point becomes traction free.

Here we define a parameter, “no-damage displacement” at time t, $\delta_{ND}^{(t)}$, and its value is obtained at the last simulation step from time $t - \Delta t$ to t:

$$ \delta_{ND}^{(t)} = \begin{cases} \delta_{coh}^{(t)}, & d^{(t)} > 0, \\ \delta_0^{(t)}, & d^{(t)} = 0, \end{cases} \quad (10) $$

where $d^{(t)}$ is the damage accumulated between time $t - \Delta t$ and $t$; $\delta_{coh}^{(t)}$ and $\delta_0^{(t)}$ are the displacement jump of the cohesive element and initial displacement at time $t$, respectively. The “no-damage displacement,” $\delta_{ND}^{(t)}$, means that if the element is unloaded and reloaded from time $t$, no further reduction of residual stiffness will happen before the displacement jump reaches this value. It is used for determining the initial displacement for the next simulation step, i.e., from time $t$ to $t + \Delta t$: $\delta_0^{(t+\Delta t)} = \delta_{ND}^{(t)}$. The traction level corresponding to this displacement jump is then defined as “no-damage traction”:

$$ \sigma_{ND}^{(t)} = K^{(t)} \delta_{ND}^{(t)} \quad (11) $$

The initial value of initial displacement at $t = 0$, $\delta_0^{(t=0)}$ is equal to $\delta_0$ that is corresponding to the inter-laminar strength $\sigma_0$. Based on the constitutive law at time $t + \Delta t$, i.e., curve OBD in Figure 2B, the displacement at which the traction becomes zero, $\delta_f^{(t+\Delta t)}$, can be determined by the residual SERR $G_{res}^{(t+\Delta t)}$, stiffness $K^{(t)}$, and initial displacement $\delta_0^{(t+\Delta t)}$:

$$ \delta_f^{(t+\Delta t)} = \frac{2G_{res}^{(t+\Delta t)}}{K^{(t)} \delta_0^{(t+\Delta t)}} \quad (12) $$

The incremental damage variable, $d$, stiffness, $K$, and traction, $\sigma$, at time $t + \Delta t$ can finally be updated as:

$$ d^{(t+\Delta t)} = \begin{cases} 0, & 0 \le \delta_{coh}^{(t+\Delta t)} < \delta_0^{(t+\Delta t)}, \\ \frac{\delta_f^{(t+\Delta t)} [\delta_{coh}^{(t+\Delta t)} - \delta_0^{(t+\Delta t)}]}{\delta_{coh}^{(t+\Delta t)} [\delta_f^{(t+\Delta t)} - \delta_0^{(t+\Delta t)}]}, & \delta_0^{(t+\Delta t)} \le \delta_{coh}^{(t+\Delta t)} < \delta_f^{(t+\Delta t)}, \\ 1, & \delta_f^{(t+\Delta t)} \le \delta_{coh}^{(t+\Delta t)}. \end{cases} \quad (13a) $$

$$ K^{(t+\Delta t)} = (1 - d^{(t+\Delta t)}) K^{(t)} \quad (13b) $$

$$ \sigma_{coh}^{(t+\Delta t)} = K^{(t+\Delta t)} \delta_{coh}^{(t+\Delta t)} \quad (13c) $$

With the definition of the incremental type damage variable, the update of constitutive curve is conducted in an incremental
way: it is only determined by the constitutive curve at the last simulation step, and the effect of history loading is included.

## COHESIVE CONSTITUTIVE LAW FOR VARIABLE MODE MIXITY

In the simulation of mixed-mode delamination, the displacement jump and traction cannot be treated as scalars, and therefore the displacement jump and traction in the cohesive constitutive equation are considered as vectors with three components. With the definition of incremental damage variable, the constitutive equation for a three-dimensional cohesive element can be expressed as follows:

$$ \sigma_i^{(t+\Delta t)} = K_{ij}^{(t+\Delta t)} \delta_j^{(t+\Delta t)}, \quad i = 1, 2, 3; \quad j = 1, 2, 3 \quad (14a) $$

$$ K_{ij}^{(t+\Delta t)} = \left[ 1 - d^{(t+\Delta t)} \left( 1 + \bar{\delta}_{1j} \frac{\langle -\delta_j^{(t+\Delta t)} \rangle}{\delta_j^{(t+\Delta t)}} \right) \right] \bar{\delta}_{ij} K^{(t)} \quad (14b) \\ i = 1, 2, 3; \quad j = 1, 2, 3 $$

Where subscript i and j denote the mode of deformation: i, j = 1-opening mode (Mode I), i, j = 2-shearing mode (Mode II) and i, j = 3-shearing mode (Mode III); $\delta_j$ (j = 1, 2, 3) is the displacement jump for the corresponding modes; $\sigma_i$ (i = 1, 2, 3) is the corresponding traction in each single mode; $K_{ij}^{(t+\Delta t)}$ is the stiffness tensor describing the relation between displacement jump vector and traction vector; $\bar{\delta}_{ij}$ is the Kronecker delta; $\langle \cdot \rangle$ is the Macaulay bracket defined as $\langle x \rangle = 1/2(|x| + x)$ to consider no damage under compression. The traction vector, $\sigma_i^{(t+\Delta t)}$, is determined by the displacement jump vector, $\delta_j^{(t+\Delta t)}$, the value of residual stiffness at time t, $K^{(t)}$, and incremental damage variable $d^{(t+\Delta t)}$ which is determined by Eq. 13a.

With the three-dimensional constitutive law, the mixed-mode deformation can be considered, as a mixture of deformations under different (pure) modes. The isotropic properties are assumed for shear modes (i.e., the plane of Mode II and III), and a displacement jump in general shear mode is defined as the resultant of the two shearing displacement jump:

$$ \delta_s = \sqrt{\delta_2^2 + \delta_3^2} \quad (15) $$

The corresponding resultant shear traction is calculated as:

$$ \sigma_s = \sqrt{\sigma_2^2 + \sigma_3^2} \quad (16) $$

An equivalent mixed-mode displacement jump is then defined as the resultant of opening and shearing displacement jump:

$$ \delta_m = \sqrt{\langle \delta_1 \rangle^2 + \delta_s^2} \quad (17) $$

The mode mixity B is defined as the ratio of supplied SERR in shear mode to the total supplied SERR:

$$ B = \frac{G_{s\_sup}}{G_{sup}} \quad (18) $$

where $G_{sup}$ is the total supplied SERR under mixed-mode loading, equal to the sum of supplied SERR in mode I, $G_{1\_sup}$, and shear mode $G_{s\_sup}$.

Different from the one-dimensional constitutive law for pure mode simulation, under the mixed-mode load condition, the constitutive law is three-dimensional as shown in Figure 3A by the black curve. It can be decomposed into two independent constitutive curves, representing the cohesive behavior under mode I and the resultant behavior under shear modes, respectively, as shown in Figure 3 by the red and green curves. The subscripts 1, s, and m in Figure 3 represent mode I, general shear mode and mixed-mode, respectively.

The development of the cohesive law in the mixed-mode space relies on the cohesive behavior in mode I and the resultant shear mode behavior. Figures 3B,C depicts the decomposed constitutive laws. The areas of the red regions in Figures 3B,C are the dissipated SERR at time t in mode I, $G_{1\_dis}^{(t)}$, and shear mode, $G_{s\_dis}^{(t)}$, which are the SERR dissipated in each mode caused by the historical loading. The supplied SERR at time t, $G_{sup}^{(t)}$, is the sum of the dissipated SERR and recoverable SERR:

$$ G_{i\_sup}^{(t)} = G_{i\_dis}^{(t)} + \frac{1}{2} K^{(t)} [\delta_i^{(t)}]^2, \quad i = 1, s \quad (19) $$

At simulation time $t + \Delta t$ the displacement jump and traction are updated based on the cohesive law, and the new dissipated SERR is calculated as the green area in Figures 3B,C, and therefore the dissipated SERR at $t + \Delta t$ is updated as:

$$ G_{i\_dis}^{(t+\Delta t)} = \begin{cases} G_{i\_dis}^{(t)} + \frac{1}{2} [K^{(t)} - K^{(t+\Delta t)}] \delta_{i0}^{(t+\Delta t)} \delta_i^{(t+\Delta t)}, & d^{(t+\Delta t)} > 0 \\ G_{i\_dis}^{(t)}, & d^{(t+\Delta t)} = 0 \end{cases}, \quad i = 1, s \quad (20) $$

The supplied SERR at time $t + \Delta t$ is equal to the dissipated SERR at time $t + \Delta t$ plus the area of blue region in Figures 3B,C:

$$ G_{i\_sup}^{(t+\Delta t)} = G_{i\_dis}^{(t+\Delta t)} + \frac{1}{2} K^{(t+\Delta t)} [\delta_i^{(t+\Delta t)}]^2, \quad i = 1, s \quad (21) $$

The mode mixity at time $t + \Delta t$ is then determined by Eq. 18. For the mixed-mode cases with only mode I and mode II deformation, the critical SERR or fracture toughness under such a mode mixity can be obtained by applying B-K criterion (Benzeggagh and Kenane, 1996):

$$ G_{c\_m}^{(t+\Delta t)} = G_{c\_I} + (G_{c\_I} - G_{c\_II}) [B^{(t+\Delta t)}]^\eta \quad (22) $$

where $G_{c\_I}$ and $G_{c\_II}$ are the fracture toughness for pure mode I and mode II, and $\eta$ is a fitting parameter fitted from the experimental results. The validation case in this study is a MMB test, where only the mixture of mode I and II is considered, and therefore the use of B-K criterion is feasible. For those interested to use the proposed mixed-mode cohesive constitutive law to the mixture of all three modes, a different criterion for calculating the fracture toughness may be needed or validation of the B-K criterion first.
<figure><img src="image_221.png" /><figcaption>FIGURE 3 | Mixed-mode constitutive curve: (A) constitutive curves in the mixed-mode space; (B) decomposition of mixed-mode constitutive curve on the mode I plane; (C) decomposition of the mixed-mode constitutive curve on the shear mode plane.</figcaption></figure>

The cohesive constitutive law for each single mode should be updated based on the latest mode mixity, i.e., $B^{(t+\Delta t)}$. However, the determination of $B^{(t+\Delta t)}$ depends on the updated cohesive law, as described by Eqs 18–21. To avoid the complicated iteration, the mode mixity at time $t$, $B^{(t)}$, is used as the value of $B^{(t+\Delta t)}$ for updating the constitutive law at time $t + \Delta t$, with an assumption that the mode mixity $B$ doesn't change significantly if the time increment $\Delta t$ is small. The updated value of $B^{(t+\Delta t)}$ is then used to obtain the cohesive law for the next simulation step from time $t + \Delta t$ to time $t + \Delta t + \Delta t'$. The feasibility of this assumption will be discussed in the next section.

The cohesive law at time $t + \Delta t$ depends on three variables: the initial stiffness, $K^{(t)}$, initial displacement, $\delta_{i0}^{(t+\Delta t)}$ ($i = 1, s$), and displacement at which traction becomes zero, $\delta_{if}^{(t+\Delta t)}$ ($i = 1, s$). The initial stiffness, $K^{(t)}$, is obtained during the previous simulation step. Similar to the one-dimensional cohesive law, the initial displacement is equal to the no-damage displacement in each mode: $\delta_{i0}^{(t+\Delta t)} = \delta_{IND}^{(t)}$ ($i = 1, s$). Their initial values at $t = 0$ are correspond to the interfacial strengths, $\sigma_{i0}$ ($i = 1, s$):

$$ \delta_{i0}^{(t=0)} = \delta_{IND}^{(t=0)} = \frac{\sigma_{i0}}{K_0} \qquad (23) $$

The displacement for the traction free condition is determined by the residual SERR $G_{res}^{(t+\Delta t)}$. At the time of $t+\Delta t$, with mode mixity $B^{(t+\Delta t)}$, if the integration point is loaded till traction is zero, the total SERR dissipated is equal to the fracture toughness $G_{c-m}^{(t+\Delta t)}$ obtained from Eq. 22. On the other hand, the fracture toughness, $G_{c-m}^{(t+\Delta t)}$, is the sum of critical SERR on mode I and shear mode under mode mixity of $B^{(t+\Delta t)}$:

$$ \begin{aligned} G_{m-c}^{(t+\Delta t)} &= G_{1-c}^{(t+\Delta t)} + G_{s-c}^{(t+\Delta t)}, \\ G_{1-c}^{(t+\Delta t)} &= G_{1-dis}^{(t)} + \frac{1}{2} K^{(t)} \delta_{10}^{(t+\Delta t)} \delta_{1f}^{(t+\Delta t)}, \\ G_{s-c}^{(t+\Delta t)} &= G_{s-dis}^{(t)} + \frac{1}{2} K^{(t)} \delta_{s0}^{(t+\Delta t)} \delta_{sf}^{(t+\Delta t)}, \end{aligned} \qquad (24) $$

The ratio of the critical SERR for the shear modes to the (mixed-mode) fracture toughness also satisfies the definition of
mode mixity, which should be identical to the value calculated in Eq. 18:

$$B^{(t+\Delta t)} = \frac{G_{s\_c}^{(t+\Delta t)}}{G_{m\_c}^{(t+\Delta t)}} \quad (25)$$

Substituting Eq. 25 into Eq. 24, the displacements at zero tractions for the opening and shear modes can be obtained as follows:

$$\delta_{1f}^{(t+\Delta t)} = \frac{2[B^{(t+\Delta t)}G_{m\_c}^{(t+\Delta t)} - G_{1\_dis}^{(t)}]}{K\delta_{10}^{(t+\Delta t)}}, \quad (26)$$

$$\delta_{sf}^{(t+\Delta t)} = \frac{2[(1 - B^{(t+\Delta t)})G_{m\_c}^{(t+\Delta t)} - G_{s\_dis}^{(t)}]}{K\delta_{s0}^{(t+\Delta t)}}.$$

Once the cohesive law on a single plane is determined, the cohesive law in the mixed-mode space (as shown in Figure 3A) can be updated. The mixed-mode initial displacement and failure displacement are calculated as the resultant of the ones under a single mode, namely $\delta_{m0}^{(t+\Delta t)} = \sqrt{(\delta_{10}^{(t+\Delta t)})^2 + (\delta_{s0}^{(t+\Delta t)})^2}$ and $\delta_{mf}^{(t+\Delta t)} = \sqrt{(\delta_{1f}^{(t+\Delta t)})^2 + (\delta_{sf}^{(t+\Delta t)})^2}$. Substituting these into Eqs 13a and 13b, the incremental damage variable d and stiffness K at time t + Δt can be obtained as in the one-dimensional cohesive law. The three-dimensional cohesive law is finally determined by Eqs. 14a and 14b.

Besides the stiffness tensor, $K_{ij}$, defined by cohesive constitutive law in Eq. 14, a tangent tensor, $\dot{K}_{ij}$, describing the relation between the displacement increment and traction increment is needed in finite element simulation for predicting the next displacement increment. It is obtained by differentiating the cohesive constitutive law in Eq. (14):

$$\dot{\sigma}_i^{(t+\Delta t)} = [1 - d^{(t+\Delta t)}]\bar{\delta}_{ij}K^{(t)}\dot{\delta}_j^{(t+\Delta t)} - [\delta_j^{(t+\Delta t)} + \bar{\delta}_{1j}\langle -\delta_1^{(t+\Delta t)} \rangle]\dot{\delta}_i^{(t+\Delta t)}, \quad i, j = 1, 2, 3 \quad (27)$$

where $\dot{\delta}_j^{(t+\Delta t)}$, $\dot{\sigma}_i^{(t+\Delta t)}$ and $\dot{d}^{(t+\Delta t)}$ are derivatives of the displacement jump vector, traction vector and damage variable. The damage variable is a function of the displacement jump as shown by Eq. 13a, and therefore its derivative should be expressed as:

$$\dot{d}^{(t+\Delta t)} = \begin{cases} \frac{\delta_{mf}^{(t+\Delta t)} \delta_{m0}^{(t+\Delta t)} [\delta_i^{(t+\Delta t)} + \bar{\delta}_{1i} \langle -\delta_1^{(t+\Delta t)} \rangle]}{[\delta_{mf}^{(t+\Delta t)} - \delta_{m0}^{(t+\Delta t)}] (\delta_m^{(t+\Delta t)})^3} \dot{\delta}_i^{(t+\Delta t)}, & \delta_{m0}^{(t+\Delta t)} < \delta_m^{(t+\Delta t)} < \delta_{mf}^{(t+\Delta t)} \\ 0, & \delta_m^{(t+\Delta t)} < \delta_{m0}^{(t+\Delta t)} \text{ or } \delta_m^{(t+\Delta t)} > \delta_{mf}^{(t+\Delta t)} \end{cases} \quad (28)$$

Substituting Eq. 28 into Eq. 27, the tangent tensor is calculated as:

$$\dot{K}_{ij} = \begin{cases} [1 - d^{(t+\Delta t)}]\bar{\delta}_{ij}K^{(t)} - \bar{\delta}_{ij}K^{(t)}D_{ij}^{(t+\Delta t)}, & \delta_{m0}^{(t+\Delta t)} < \delta_m^{(t+\Delta t)} < \delta_{mf}^{(t+\Delta t)} \\ [1 - d^{(t+\Delta t)}]\bar{\delta}_{ij}K^{(t)}, & \delta_m^{(t+\Delta t)} < \delta_{m0}^{(t+\Delta t)} \text{ or } \delta_m^{(t+\Delta t)} > \delta_{mf}^{(t+\Delta t)}, \end{cases} \\ D_{ij}^{(t+\Delta t)} = \frac{[\delta_j^{(t+\Delta t)} + \bar{\delta}_{1j} \langle -\delta_1^{(t+\Delta t)} \rangle] [\delta_i^{(t+\Delta t)} + \bar{\delta}_{1i} \langle -\delta_1^{(t+\Delta t)} \rangle] \delta_{mf}^{(t+\Delta t)} \delta_{m0}^{(t+\Delta t)}}{[\delta_{mf}^{(t+\Delta t)} - \delta_{m0}^{(t+\Delta t)}] (\delta_m^{(t+\Delta t)})^3} \quad (29)$$

## RESULTS OF FEM VALIDATION AND DISCUSSION

The developed formulation of the mixed-mode cohesive constitutive law is tested by simulating mixed-mode delamination in fibre-reinforced composites. The incremental mixed-mode cohesive law was implemented into the cohesive element as a user-defined material (UMAT) subroutine in FE software Abaqus Implicit. In the simulation, 8-node solid element C3D8I was chosen for modeling composite plies and 8-node cohesive element COH3D8 was applied to the inter-laminar interface. The cohesive elements have a zero thickness.

### FEM Model Description

The MMB test defined in ASTM D6671 (ASTM D6671, 2013) is frequently used to characterize the mixed-mode delamination behavior, as shown in Figure 1. The MMB test set-up is designed such as to combine Mode I and II loadings. The supplied SERR at the crack tip can be calculated as (Reeder and REWS, 1990):

$$\begin{align} G_I &= \frac{3P_{\text{lever}}^2 (3c - L)^2}{4b^2 h^3 L^2 E_{1f}} (a + \chi h)^2, \notag \\ G_{II} &= \frac{9P_{\text{lever}}^2 (c + L)^2}{16b^2 h^3 L^2 E_{1f}} (a + 0.42\chi h)^2, \tag{30} \\ G_m &= G_I + G_{II}, \notag \end{align}$$

where $G_I$ and $G_{II}$ are the Mode I and II components of SERR; $G_m$ is the total mixed-mode SERR; $P_{\text{lever}}$ is the applied load at the end of the load lever; a is the crack length; b and h are the specimen width and half thickness, L is the half span length, and c is the lever length, χ is the crack length correction parameter calculated by Eq. 3, and $E_{1f}$ is the bending modulus of the composite:

$$E_{1f} = \frac{8(a + \chi h)^3 + [6(a + 0.42\chi h)^3 + 4L^3](c + L)^2}{16L^2 b h^3 C} \quad (31a)$$

where C is the compliance obtained from the load displacement curve. The mode mixity B at the specimen level is then calculated as $B = G_{II}/G_m$.

The displacement applied at the end of the specimen $U_{\text{end}}$, in the middle point $U_{\text{mid}}$ and at the end of the load lever $U_{\text{lever}}$ satisfy the following relation (Camanho et al., 2003):

$$U_{\text{lever}} = \left(\frac{3c - L}{4L}\right) U_{\text{end}} + \left(\frac{c + L}{L}\right) \left(U_{\text{mid}} + \frac{1}{4}U_{\text{end}}\right) \quad (32)$$
<figure><img src="image_222.png" /><figcaption>TABLE 1 | Material properties of AS4/PEEK UD laminates (Turon et al., 2010).</figcaption></figure>

<table><thead><tr><th></th><th>E<sub>11</sub><br>(GPa)</th><th>E<sub>22</sub> = E<sub>33</sub><br>(GPa)</th><th>v<sub>12</sub> = v<sub>13</sub><br>(-)</th><th>v<sub>23</sub><br>(-)</th><th>G<sub>12</sub> = G<sub>13</sub><br>(GPa)</th><th>G<sub>23</sub><br>(GPa)</th></tr></thead><tbody><tr><td>Elastic<br>properties</td><td>122.7</td><td>10.1</td><td>0.25</td><td>0.45</td><td>5.5</td><td>3.7</td></tr><tr><td>Interlaminar<br>properties</td><td>G<sub>c, I</sub><br>(N/mm)</td><td>G<sub>c, II</sub><br>(N/mm)</td><td>η</td><td>σ<sub>10</sub><br>(MPa)</td><td>σ<sub>20</sub><br>(MPa)</td><td></td></tr><tr><td></td><td>0.969</td><td>1.719</td><td>2.0</td><td>80</td><td>100</td><td></td></tr></tbody></table>

The model to simulate the MMB test is shown in **Figure 4**. To simplify the simulation, the displacement is applied at three positions respectively, instead of simulating the load lever explicitly. A constraint equation is defined to meet the boundary condition defined by Eq. 32, which can be viewed as the application of a virtual load lever.

The same geometry of MMB specimen simulated by Turon et al. (2018) is used in this study. It is 102.0 mm long and 3.12 mm thick with a precrack placed at the middle plane. The material properties of the laminate are listed in **Table 1**.

The details of the geometry model and applied boundary conditions are illustrated in **Figure 4**. A layer of cohesive elements with damage variable $d = 1$ are placed at the position of the pre-crack to simulate the traction free surface, which cannot hold tension or shear but compression to prevent the penetration between the upper and lower arms. The potential path for the delamination propagation is simulated as one layer of cohesive elements with the incremental cohesive constitutive law proposed in *Cohesive Constitutive Law for Variable Mode Mixity* and the input parameters of the law are reported in **Table 1**. The upper and lower arms of MMB specimen were simulated by solid elements. To accelerate the simulation, only one element is placed along the width direction of the specimen, i.e., *z* direction in **Figure 4**, and the width of the model was equal to the element width 0.1 mm. A state of generalized plane strain was assumed in the width direction: the front and back edges of the model were kept parallel during the deformation. This assumption means that the displacement $δ_z$ is constant throughout the *x-y* plane of the model. A fine mesh along the *x* direction, was set at the end of the pre-crack and beginning of the potential delamination path, depicted as gray zone in

<table><caption>TABLE 2 | Geometry parameters of MMB tests.</caption><thead><tr><th>Mode mixity at specimen level<br>B = G<sub>II</sub>/(G<sub>I</sub> + G<sub>II</sub>)</th><th>20%</th><th>50%</th><th>80%</th></tr></thead><tbody><tr><td>Pre-crack a<sub>0</sub> (mm)</td><td>33.70</td><td>34.10</td><td>31.40</td></tr><tr><td>Lever length c (mm)</td><td>109.89</td><td>44.60</td><td>28.47</td></tr></tbody></table>

<table><caption>TABLE 3 | Four sets of values for the interlaminar strength used in the simulation.</caption><thead><tr><th>Set no</th><th>Set1</th><th>Set2</th><th>Set3</th><th>Set4</th></tr></thead><tbody><tr><td>σ<sup>static</sup><sub>10</sub> (MPa)</td><td>80</td><td>80</td><td>80</td><td>30</td></tr><tr><td>σ<sup>static</sup><sub>50</sub> (MPa)</td><td>100</td><td>106.55</td><td>120</td><td>100</td></tr></tbody></table>

**Figure 4.** The preliminary mesh sensitivity study shows that the simulation result is not sensitive to the mesh size when it is smaller than 0.05 mm. Therefore, a value of 0.05 was chosen as the refined mesh size. Coarse mesh was used in the other parts, where the element length was 0.2 mm. Four layers of solid elements were placed along the thickness direction, i.e., *y* direction, on the upper and lower arms respectively.

The MMB test with a mode mixity *B* at the specimen level equal to 20%, 50% and 80% are simulated. According to **30** and **31a**, the corresponding pre-crack length $a_0$ and lever length *c* are listed in **Table 2**. The displacement is applied at the end of the load lever $U_{lever}$, and the reaction force obtained at the same position is recorded.

## Effect of the Inter-Laminar Strength

The analytical solution of the reaction force and displacement can be obtained by Eqs 1–3 based on LEFM (ASTM D6671, 2013). Four sets of values for the inter-laminar strength were chosen in this study for testing its effect on the delamination growth at a specimen level (Table 3).

The inter-laminar strength in Set1 is used by Turon. et al. (2006) for carbon fiber/PEEK composite AS4/PEEK; the strengths in Set2 is determined to satisfy the condition defined by Eq. 4 to eliminate the strength dependency in Turon's model (Turon et al., 2010); the values in Set3 and Set4 are chosen arbitrarily. The resulting
<figure><img src="image_223.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_224.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_225.png" /><figcaption>**FIGURE 5** | The dependence of the load-displacement curve on the inter-laminar normal and shear strength with: (A) mode mixity of 20%; (B) mode mixity of 50%; (C) mode mixity of 80%.</figcaption></figure>

reaction force-displacement curves with these different sets of the inter-laminar strength are plotted in **Figure 5**, where the reaction force at the lever end is expressed as the linear force along the width direction (i.e., *z* direction). No matter which set of inter-laminar strength is used, the simulation results meet with each other and all have a good agreement with the analytical solution under all the mode mixity load conditions. In the simulations of mode mixity *B* equal to 50% and 80%, there is no dependence of structural response on the selected inter-laminar strength. In the case of *B* = 20%, the selection of inter-laminar strength only has a minor effect on the maximum force, and the largest difference on the maximum force is 2.22%. The propagation behavior i.e., structural response after the maximum load, has no dependence on the strength value in this case. With the developed cohesive constitutive law, there is no need to select a special strength value for a correct mixed-mode delamination behavior.

In this study, the mode mixity is calculated using the supplied SERR, which is as the ratio between the supplied SERR for shearing to the total supplied SERR expressed by Eq. 18. This is identical to the definition of the mode mixity described in the MMB test standard (ASTM D6671, 2013). The definition of the mode mixity used by Turon et al. (2006) is based on the displacement jump: $B = \delta_s^2/\delta_m^2$,

which gives very different predictions of the mode mixity from the supplied SERR based on Eq. 18. The difference becomes even more significant when the mode mixity at the integration point undergoes significant changes during the loading history. Only when the cohesive law for each single mode mixity has a similar shape, the two definitions of the mode mixity produce the same results. The difference in definition of mode mixity explains why Turon's model can only predict the correct structural response of MMB test when a specific relation between the mode I and shear mode inter-laminar strength is satisfied.

## Variable Mode Mixity at the Integration Point

MMB test is designed for characterizing the mixed-mode delamination behavior under fixed mode mixity at specimen level. The mode mixity is determined by specimen geometry, the lever length *c* and material elastic properties, as depicted by Eq. 30. In the current simulations, all these parameters are fixed to obtain a mode mixity of 20%, 50% and 80% at specimen level. However, the mode mixity at the integration point does not stay constant during the loading process. This phenomenon has also been observed by Turon et al. (2010).

<figure><img src="image_226.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_227.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_228.png" /><figcaption>**FIGURE 6** | Change of the mode mixity at the integration point which is 5 mm away from the pre-crack tip with: (A) mode mixity of 20%; (B) mode mixity of 50%; (C) mode mixity of 80%.</figcaption></figure>
<figure><img src="image_229.png" /></figure>

To illustrate how the mode mixity changes, it was evaluated at the integration point 5 mm away from a pre-crack tip during the whole loading history. It is shown in **Figure 6**, where the x axis is the normalized displacement jump (that includes both normal and shear displacements) at the integration point. The displacement jump is normalized by the value of the displacement jump when the traction at the integration point is zero. As shown in **Figure 6**, similarly to the observation of Turon et al. (2010), the mode mixity is not fixed as the beam theory at specimen level predicts. Instead it decreases monotonically approaching the level close to the one defined at the specimen level when the integration point becomes traction free, i.e., 20%, 50% and 80% for **Figures 6A-C**, respectively. Because the predicted structural behavior from the proposed cohesive law has little sensitivity to the inter-laminar strength, the discrepancy of mode mixity at the specimen level and at the integration point level cannot be attributed to the dependence of delamination behavior on the inter-laminar strength, which was used to explain the dependence of predicted load-displacement on the selection of the inter-laminar strength in Turon's study (Turon et al., 2010).

## Mixed-Mode Constitutive Law at Integration Point

Due to the variable mode mixity at the integration point, the constitutive behavior obtained at the integration point also differs from the one for a fixed mode mixity. **Figure 7** shows the mixed mode constitutive law at the integration point which is 5 mm to a pre-crack tip. The mixed mode constitutive law at the integration point deviates from the typical triangle shape for

the case of a fixed mode mixity. The total SERR dissipated when the integration point becomes traction free is output, which can be considered as the fracture toughness of the interface under such a load condition. The values of fracture toughness used in LEFM analysis for the fixed mode mixity defined at specimen level are also listed. It can be found that these two values agrees with each other for all the mode mixities and sets of inter-laminar strengths, as listed in **Table 4**. The similar values of fracture toughness output from FEM and used in LEFM analysis for fixed mode mixity at specimen level explains why the predicted reaction force-displacement curve at specimen level has a good agreement with the LEFM analysis even if the mode mixity at integration point deviates from the one defined at specimen level for all the mode mixities and selections of inter-laminar strengths.

## Feasibility of the Assumption of Small Mode Mixity Change

In *Cohesive Constitutive Law for Variable Mode Mixity*, an assumption was made for simplifying the development of three-dimensional cohesive constitutive law: the change of mode mixity is small, and therefore the mode mixity obtained from previous simulation step is directly used as the mode mixity for the current simulation increment. The error between the actual mode mixity in the simulation increment and the value used for developing the cohesive law for this simulation increment is calculated, as shown in **Figure 8**. The results are from the simulation with the inter-laminar strengths from Set1. It can be found that the error of the mode mixity is always below 7% for all the mode mixities.

## CONCLUSION AND OUTLOOK

A new formulation of a cohesive constitutive law was developed in this study, which is capable of simulating delamination with variable mode mixity and has little sensitivity to the inter-laminar strength. An incremental type of damage variable is proposed, a factor accounting for the reduction of the residual stiffness in

<table><caption>**TABLE 4** | Fracture toughness of interface output from different simulation cases.</caption><thead><tr><th>Mode-mixity at specimen level</th><th>LEFM<br>(N/mm)</th><th>Set1<br>(N/mm)</th><th>Set2<br>(N/mm)</th><th>Set3<br>(N/mm)</th><th>Set4<br>(N/mm)</th></tr></thead><tbody><tr><td>20%</td><td>1.00</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.99</td></tr><tr><td>50%</td><td>1.16</td><td>1.15</td><td>1.15</td><td>1.15</td><td>1.15</td></tr><tr><td>80%</td><td>1.45</td><td>1.44</td><td>1.43</td><td>1.45</td><td>1.44</td></tr></tbody></table>
<figure><img src="image_230.png" /><figcaption>**FIGURE 8** | Error of mode mixity from the assumption of small change of mode mixity.</figcaption></figure>

the traction-displacement cohesive law. It allows updating the constitutive law step by step (incrementally) and includes the effect of loading history. The mode mixity at the integration point is updated during the loading, and calculated in the formulation of the supplied SERR.

The model was validated by simulating MMB test with three different mode mixities. Simulations were run for four different values of the inter-laminar strength, and for these simulation cases, the reaction force-displacement for the loaded MMB specimen has a good agreement with the analytical solution from LEFM analysis. There is only a minor dependence on the inter-laminar strength. In the case of *B* = 20% the largest difference in the maximum load is 2.22%, and this dependence disappears in the propagation part after the maximum load. The dependence on the strength is eliminated in the simulations with mode mixity of 50% and 80%. The mode mixity at the integration point is found to be different from the constant mode mixity applied macroscopically at the specimen level. However, the total SERR dissipated (up to the displacement when the integration point becomes traction free) is almost equal to the fracture toughness input to LFEM analysis for the fixed mode mixity

## REFERENCES

Alfano, G., and Crisfield, M. A. (2001). Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. *Int. J. Numer. Methods Eng.* 50, 1701–1736. doi:10.1002/nme.93

ASTM D6671/D6671M-13 (2013). *Standard test method for mixed mode I-mode II interlaminar fracture toughness of unidirectional fiber reinforced polymer matrix composites*. West Conshohocken, PA: ASTM International.

Barenblatt, G. I. (1959). The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks. *J. Appl. Math. Mech.* 23, 622–636. doi:10.1016/0021-8928(59)90157-1

Barenblatt, G. I. (1962). The mathematical theory of equilibrium cracks in brittle fracture. *Adv. Appl. Mech.* 7, 55–129. doi:10.1016/s0065-2156(08)70121-2

at specimen level and is independent of the inter-laminar strength. The agreement of the output fracture toughness with the fracture toughness used in LEFM analysis for fixed mode mixity at specimen level explains the correctly predicted reaction force-displacement curve for all the studied cases.

For those who are interested to perform characterization under true constant mode-mixity conditions Sørensen et al. (2006) proposed a test method where a double cantilever beam specimen is loaded with uneven bending moments (DCB-UBM). An FEM study of the DCB-UBM test with the proposed incremental cohesive law would be interesting to confirm that the mode mixities are indeed the same at the macro- and micro-scale in this test.

The validation case in this study is an MMB test, for which experimental data are available. A case study with all three modes of delamination supported by experimental data is needed to complete the validation of the proposed mixed-mode cohesive law.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article can be made available upon request by the authors, without undue reservation.

## AUTHOR CONTRIBUTIONS

MZ: model development, simulation, analysis, article writing; LG: data interpretation, supervision, article editing; SL: data interpretation, supervision, article editing.

## FUNDING

The work has been funded by the SBO project "M3Strength," which fits in the MacroModelMat (M3) research program, coordinated by Siemens (Siemens Digital Industries Software, Belgium) and funded by SIM (Strategic Initiative Materials in Flanders) and VLAIO (Flanders Innovation & Entrepreneurship). S.V. Lomov holds Toray Chair for Composite Materials at KU Leuven, support from which is gratefully acknowledged.

Bennati, S., Fisicaro, P., and Valvo, P. S. (2013a). An enhanced beam-theory model of the mixed-mode bending (MMB) test-part I: literature review and mechanical model. *Meccanica* 48, 443–462. doi:10.1007/s11012-012-9686-3

Bennati, S., Fisicaro, P., and Valvo, P. S. (2013b). An enhanced beam-theory model of the mixed-mode bending (MMB) test-part II: applications and results. *Meccanica* 48, 465–484. doi:10.1007/s11012-012-9682-7

Benzeggagh, M. L., and Kenane, M. (1996). Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. *Compos. Sci. Technol.* 56, 439–449. doi:10.1016/0266-3538(96)00005-x

Borg, R., Nilsson, L., and Simonsson, K. (2002). Modeling of delamination using a discretized cohesive zone and damage formulation. *Compos. Sci. Technol.* 62, 1299–1314. doi:10.1016/s0266-3538(02)00070-2
Camanho, P. P., Davila, C. G., and De Moura, M. F. (2003). Numerical simulation of mixed-mode progressive delamination in composite materials. *J. Compos. Mater.* 37, 1415–1438.

De Moura, M. F. S. F., Gonçalves, J. P. M., and Silva, F. G. A. (2016). A new energy based mixed-mode cohesive zone model. *Int. J. Solid Struct.* 102-103, 112–119. doi:10.1016/j.ijsolstr.2016.10.012

Dugdale, D. S. (1960). Yielding of steel sheets containing slits. *J. Mech. Phys. Solid.* 8, 100–104. doi:10.1016/0022-5096(60)90013-2

González, E. V., Maimí, P., Camanho, P. P., Turon, A., and Mayugo, J. A. (2012). Simulation of drop-weight impact and compression after impact tests on composite laminates. *Compos. Struct.* 94, 3364–3378. doi:10.1016/j.compstruct.2012.05.015

Gözlüklü, B., and Coker, D. (2012). Modeling of the dynamic delamination of L-shaped unidirectional laminated composites. *Compos. Struct.* 94, 1430–1442. doi:10.1016/j.compstruct.2011.11.015

Gozluklu, B., Uyar, I., and Coker, D. (2015). Intersonic delamination in curved thick composite laminates under quasi-static loading. *Mech. Mater.* 80, 163–182. doi:10.1016/j.mechmat.2014.07.013

Harper, P. W., Sun, L., and Hallett, S. R. (2012). A study on the influence of cohesive zone interface element strength parameters on mixed mode behaviour. *Compos. Appl. Sci. Manuf.* 43, 722–734. doi:10.1016/j.compositesa.2011.12.016

Jensen, S. M., Martos, M. J., Bak, B. L. V., and Lindgaard, E. (2019). Formulation of a mixed-mode multilinear cohesive zone law in an interface finite element for modelling delamination with R-curve effects. *Compos. Struct.* 216, 477–486. doi:10.1016/j.compstruct.2019.02.029

Joki, R. K., Grytten, F., Hayman, B., and Sørensen, B. F. (2016). Determination of a cohesive law for delamination modelling-accounting for variation in crack opening and stress state across the test specimen width. *Compos. Sci. Technol.* 128, 49–57. doi:10.1016/j.compscitech.2016.01.026

Kawashita, L. F., and Hallett, S. R. (2012). A crack tip tracking algorithm for cohesive interface element analysis of fatigue delamination propagation in composite materials. *Int. J. Solid Struct.* 49, 2898–2913. doi:10.1016/j.ijsolstr.2012.03.034

Liu, Q., Gorbatikh, L., and Lomov, S. V. (2019). A combined use of embedded and cohesive elements to model damage development in fibrous composites. *Compos. Struct.* 223, 110921. doi:10.1016/j.compstruct.2019.110921

Lu, X., Ridha, M., Chen, B. Y., Tan, V. B. C., and Tay, T. E. (2019). On cohesive element parameters and delamination modelling. *Eng. Fract. Mech.* 206, 278–296. doi:10.1016/j.engfracmech.2018.12.009

Nojavan, S., Schesser, D., and Yang, Q. D. (2016). A two-dimensional *in situ* fatigue cohesive zone model for crack propagation in composites under cyclic loading. *Int. J. Fatig.* 82, 449–461. doi:10.1016/j.ijfatigue.2015.08.029

Reeder, J. R., and Crews, J. H. (1990). Mixed-mode bending method for delamination testing. *AIAA J.* 28, 1270–1276. doi:10.2514/3.25204

Sørensen, B. F., Jørgensen, K., Jacobsen, T. K., and Østergaard, R. C. (2006). DCB-specimen loaded with uneven bending moments. *Int. J. Fract.* 141, 163–176. doi:10.1007/s10704-006-0071-x

Turon, A., Camanho, P. P., Costa, J., and Dávila, C. G. (2006). A damage model for the simulation of delamination in advanced composites under variable-mode loading. *Mech. Mater.* 38, 1072–1089. doi:10.1016/j.mechmat.2005.10.003

Turon, A., Camanho, P. P., Costa, J., and Renart, J. (2010). Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: definition of interlaminar strengths and elastic stiffness. *Compos. Struct.* 92, 1857–1864. doi:10.1016/j.compstruct.2010.01.012

Turon, A., González, E. V., Sarrado, C., Guillamet, G., and Maimí, P. (2018). Accurate simulation of delamination under mixed-mode loading using a cohesive model with a mode-dependent penalty stiffness. *Compos. Struct.* 184, 506–511. doi:10.1016/j.compstruct.2017.10.017

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Zhu, Gorbatikh and Lomov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_231.png" /></figure>

Check for
updates

# Characterization of Surface Fatigue Crack Nucleation and Microstructurally Small Crack Growth in High Strength Aluminum Alloys

Robert Fleishel<fnref n="1" />, Cole Cauthen<fnref n="2" />, Steven Daniewicz³, Andrew Baker⁴, J. Brian Jordon<fnref n="2" /> and Stephanie TerMaath<fnref n="1" />*

It is well established that fatigue crack nucleation and small crack growth in high strength aluminum alloys are highly influenced by the surrounding microstructure including grain boundaries, texture, inclusion barriers, among other factors. As such, specific and targeted experimental and computational methods are necessary to accurately capture and predict the discrete behavior of microstructurally small fatigue cracks. In this study, surface fatigue crack nucleation and microstructurally small crack growth in high strength aluminum alloys, commonly used in aerospace applications, are evaluated through a holistic approach encompassing fatigue testing, crack measurement, and computational prediction of crack growth rates. During fatigue testing, crack shapes and growth are quantified using a novel surface replication technique that is applied to investigate crack nucleation, as well as to collect validation data that includes an accurate description of crack shape during crack propagation, a challenging and essential component in predicting crack growth. Computational simulation of fatigue crack growth in non-straight, complex surface crack arrays typically requires high fidelity analysis using computationally expensive methods to account for the mathematical and geometrical complexities inherent in the solution. A dislocation distribution based technique has been previously demonstrated to rapidly and accurately predict the stress intensity factors for through cracks of complex shape. This method was expanded and investigated as an approach for rapidly predicting the crack growth rate of kinked and tortuous surface crack arrays, using the crack configuration and bulk material properties as inputs. To investigate the accuracy and effectiveness of this characterization approach, surface crack growth in AA7075-T7351 was experimentally analyzed and modeled under high cycle and low cycle fatigue conditions. This comprehensive approach was determined to be an expedient and applicable method for characterizing and evaluating the nucleation and crack growth rate of non-planar microstructurally small and short crack configurations.

**Keywords:** crack growth rate analysis, fatigue-cracks, aluminum alloy, computational modeling, surface crack behavior

OPEN ACCESS

**Edited by:**

Patricia Krawczak,
IMT Lille Douai, France

**Reviewed by:**

Frédéric Charles Lebon,
Aix-Marseille Université, France
Andrea Spagnoli,
University of Parma, Italy

*Correspondence:

Stephanie TerMaath
stermaat@utk.edu

**Specialty section:**

This article was submitted to
Mechanics of Materials,
a section of the journal
Frontiers in Materials

**Received:** 02 August 2020

**Accepted:** 10 December 2020

**Published:** 11 February 2021

**Citation:**

Fleishel R, Cauthen C, Daniewicz S, Baker A, Jordon J B and TerMaath S (2021) Characterization of Surface Fatigue Crack Nucleation and Microstructurally Small Crack Growth in High Strength Aluminum Alloys. *Front. Mater.* 7:590747
doi: 10.3389/fmats.2020.590747

<footnote n="1">¹Department of Mechanical, Aerospace, and Biomedical Engineering, The University of Tennessee, Knoxville, TN, United States,</footnote>

<footnote n="2">²Department of Mechanical Engineering, The University of Alabama, Tuscaloosa, AL, United States, ³Department of Metallurgical</footnote>

<footnote>and Materials Engineering, The University of Alabama, Tuscaloosa, AL, United States, ⁴Boeing Research and Technology, The</footnote>

<footnote>Boeing Company, Berkeley, MO, United States</footnote>
# INTRODUCTION

Predicting the crack growth rate and direction in the early stages of crack growth is challenging due to several phenomena. The rate and direction of crack growth immediately following initiation is heavily dependent on microstructure (Hussain 1997; Stephens et al., 2001), and both can change rapidly as a crack grows across grain boundaries, through imperfections, and around inclusions. This behavior can cause microstructurally small cracks to kink or for several smaller cracks to form, later coalescing into a single dominant crack, in the early stages of crack growth. To characterize crack nucleation and to predict the growth of small cracks requires accurate correlation of crack shape with crack growth rate and a computational approach that accounts for complex crack shapes due to direction changes and crack branching. The objective of this research is to evaluate the accuracy and efficiency of a comprehensive approach to characterizing crack nucleation and predicting the crack growth rate of small cracks using a novel surface replication technique to obtain high resolution crack shapes and a computational predictive capability for the propagation of such non-straight cracks. The research was performed through experimental and computational components.

Recent analytical modeling of microstructurally small surface cracks has been successfully performed (Cauthen et al., 2017), in which a strip yield model (Daniewicz and Aveline, 2000) and CTOD (crack tip opening displacement) method were applied for a straight surface fatigue crack. An analytical multistage fatigue model, also assuming a straight crack, has been used to predict microstructurally small crack growth of aluminum alloys (McDowell et al., 2003; Xue et al., 2007; Bernard et al., 2013; McCullough et al., 2019). Research has also been performed to explore the non-straight nature of microcracks, including the effects of grain and crystal orientations (Herbig et al., 2011; Castelluccio and McDowell, 2012; Castelluccio and McDowell, 2016; Wilson et al., 2018). While many methods are available for predicting crack growth that account for both microstructure and crack shape using high fidelity modeling, including finite element methods (Wilson and Dunne, 2019), boundary element methods (Sun and Ma, 1997), and peridynamics (Silling and Abe, 2014), these methods require explicit modeling of the microstructure and significant computational time for a single analysis. The need to know the exact microstructure present around the crack tip makes applying such methods to existing crack configurations, especially in time critical situations such as during an inspection to determine remaining life, difficult.

Experimental approaches for quantification of nucleation and fatigue crack growth mechanisms in metals are well documented. Some methods employ micronotches where the initiation event can be readily located and observed. However this method does limit the observation of naturally occurring cracks where the initiation mechanisms which come from microstructural features or defects, such as intermetallic particles and persistent slip bands (PSBs) can be quantified (Gall et al., 1999; Xue et al., 2007; Ma et al., 2019; Zerbst et al., 2019; Dai et al., 2020). As such, without knowing the location of where the fatigue crack initiated, taking in-situ measurements without the use of a micronotch is

improbable. An alternate approach that allows for naturally occurring fatigue cracks to initiate from microstructural features and/or defects is the surface replica approach. The main idea of the surface replica technique is to take intermittent surface replicas of the fatigue specimen over the lifetime of the test. Once the test is complete, analysis of the surface replicas can be examined in reverse order since the fatigue crack(s) that lead to ultimate fracture are the largest and easy to identify under low magnification and tracking them is straightforward. A relatively new technique for taking surface replicas uses a fast drying two-part silicon rubber compound that does not experience shrinkage compared to the more traditional acetate tape method (Newman et al., 1994). The procedure to make the surface replicas using the fast drying two-part silicon rubber compound, commercially known as Repliset, was first demonstrated by Newman et al. (2009). To demonstrate the use of fast drying two-part silicon rubber compound on quantifying fatigue mechanisms, previous work (Jordon et al., 2012) showed that the method could be used to quantify the fatigue crack nucleation stage and small crack growth rates in a high strength aluminum alloy. While the fast drying two-part silicon rubber replica method is not ideal for directly measuring subsurface crack growth, it has been used to further understand the mechanisms of surface dominated crack nucleation and growth from naturally occurring cracks. Recently, other metals have been examined using the fast drying two-part silicon rubber replica method including titanium (Hu et al., 2019; Zhu et al., 2019), magnesium alloys (Bernard et al., 2013; Cauthen et al., 2017), stainless steel (Deng et al., 2016), and nickel based super alloys (Deng et al., 2015; Qin et al., 2015; Deng et al., 2016; Zhu et al., 2016; Ye et al., 2017; Zhu et al., 2018), where these studies explored the effects of variation in grain size, maximum stress, and stress ratios (R values) on fatigue crack growth behavior. In addition, the fast drying two-part silicon rubber replica method has also been employed on a range of specimen geometries from standard notched specimen types (Bernard et al., 2013; Cauthen et al., 2017; Zhu et al., 2018), to laser created micro notches (Qin et al., 2015; Zhu et al., 2019), and aerospace components (Newman et al., 2009; Altıntaş and Güngör, 2013). The present study focuses on using the novel fast drying two-part silicon rubber replica method to quantify the fatigue crack nucleation and microstructurally small crack growth mechanisms of AA7075-T7351. In particular, the use of the surface replica approach provides direct and high resolution observation of both the surface crack initiation event and the subsequent propagation of small cracks where the crack kinking can be directly quantified. One of the advantages of this experimental method is the opportunity to directly correlate crack growth to sources of crack initiation that are also confirmed under high resolution scanning electron microscopy. As such, to the best of the authors' knowledge, this is the first investigation to quantify surface crack nucleation that is directly linked to microstructurally small fatigue crack growth in AA7075-T7351 using the two-part silicon rubber replica method.

Using these measurements, the fatigue crack growth rates (FCGR) are numerically determined and compared to the experimental results to validate a fatigue crack growth model
using a dislocation distribution based technique (DDT). DDT has been previously demonstrated to rapidly and accurately predict the stress intensity factors for a variety of problems (Lardner 1974; Barber 1992; Hills et al., 1996), including long, through crack arrays of complex shape (TerMaath et al., 2006; Yavuz et al., 2006). Dislocation distributions have also been used to model fatigue crack growth behavior (Yavuz et al., 2020). This method was expanded to predict the fatigue crack growth rate (FCGR) of short, surface cracks. To develop a rapid prediction capability for predicting the fatigue crack growth rate (FCGR) in kinked and tortuous surface crack arrays, a DDT was expanded to include short, surface crack analysis using the crack configuration and bulk material properties as inputs. While a goal of this research was to evaluate the DDT using experimental data that was correlated to the replica data, this method was chosen because it has been demonstrated for large arrays of cracks of complex shape and is robust as cracks grow near each other and coalesce, an important feature for future development and applicability to real world structures. FCGR prediction is achieved by combining a method for calculating the stress intensity factor (SIF, designated as *K* for linear elastic behavior) of complex shaped cracks (TerMaath, 2001) with empirical SIF based crack growth models (Paris, 1961) and geometric factors to account for surface crack depth (Raju and Newman, 1979). The use of a geometric factor vs. a strip-yield model to account for surface cracking reduces the complexity of the analysis and is less computationally expensive, an important consideration for future development for large arrays of surface cracks. The use of bulk material crack growth data enables this model to be applicable without specialized data collection or knowledge of the exact microstructure surrounding the crack due to the fact that *K*-based experimental crack growth data is more commonly available than CTOD experimental data. This model is explored for its applicability to determine the FCGR of small cracks in non-straight configurations, with the change in crack driving force, the stress intensity factor, accounting for the non-straight nature of the crack.

## MATERIALS AND METHODS

To investigate the applicability of the integrated experimental and computational approach for characterizing the fatigue crack growth behavior of kinked surface cracks, microstructurally small crack growth in AA7075-T7351 was experimentally analyzed and computationally modeled under high cycle and low cycle fatigue conditions. This legacy material was chosen to demonstrate the characterization approach due to its widespread use and well-established microstructural variability. The crack growth was measured on the surface of a shallow notched flat fatigue specimen using a novel, fast drying two-part silicon rubber surface replication technique. Fractography and Electron Backscatter Diffraction (EBSD) mapping of specimens were performed post-fracture to determine the source of crack initiation and distribution of intermetallic particles within the material, as well as to characterize the microstructure at the location of crack

initiation. The crack configurations captured via surface replication were then input into the computational model to calculate the crack growth rate at varying stages of crack propagation and the computational results were compared to the experimental data.

## Experimental Testing

The material used for testing in this study was a 63.5 mm thick rolled plate of AA7075-T7351. In this study, the fatigue specimens were cut out of the plate in the longitudinal orientation at the quarter thickness (t/4) and half thickness (t/2) in order to assess the effect of location dependent mechanisms in the rolled plate. This specimen design used in this study was a flat specimen with a reduced gauge, where this reduced gauge contains a shallow notch to facilitate fatigue crack initiation (Figure 1). While most of the specimen geometry was machined using wire electrical discharge machining (EDM), the shallow notch was prepared using a traditional milling tool. In addition, the shallow notch was machined into one side of the specimen geometry to facilitate natural crack initiation within a general location. Following machining, all surfaces were lightly polished using progressively finer polishing paper until a smooth surface was achieved. This specimen geometry was first used by (Tokaji et al., 2004) and has been successfully implemented to measure microstructurally small surface fatigue crack growth for magnesium alloys (Bernard et al., 2013; Cauthen et al., 2017; Cauthen et al., 2020).

Fatigue crack growth was evaluated for the test specimens in both the low cycle fatigue (LCF) and high cycle fatigue (HCF) regime at a load ratio of $R = 0.1$ and a test frequency below 10 Hz. Load levels for these regimes were chosen such that cycles to failure for LCF testing was approximately $10^4$ cycles, and cycles to failure for HCF testing was approximately $10^5$ cycles. Preliminary testing was conducted to find these desired load levels, and a stress range of 420 and 260 MPa was chosen for LCF and HCF

<figure><img src="image_232.png" /><figcaption>**FIGURE 1** | (A) Flat tensile specimen with shallow notch used for fatigue testing. (B) Specimen dimensions (mm).</figcaption></figure>
<figure><img src="image_233.png" /><figcaption>FIGURE 2 | Examples of idealized kinked crack configurations overlaid on surface replication images. (A) a crack with two segments and a single kink, (B) a crack with three segments and two kinks.</figcaption></figure>

testing, respectively. Eight samples, four in the HCF and four in the LCF regimes, were tested. The fatigue crack growth tests were conducted in ambient laboratory air and relative humidity and were terminated when full separation of the specimen was achieved. The application of the replica method for quantifying fatigue crack growth required the regular interruption of the fatigue tests at a mean cyclic load. Once stopped, the fatigue sample was dusted with particle free compressed air to remove any surface contaminates. Following the procedure detailed by Jordon et al. (2012), the Repliset® material, which is a two-part silicon rubber replica, was dispensed on to the shallow notch of the fatigue specimen using a mixing nozzle and dispensing gun. Tape was placed over the Repliset® material to hold it in place during the curing phase, where the curing time for the Repliset® material was approximately 5–10 min. Once cured, the tape was carefully pulled away from the specimen revealing an inverted replication of the surface of the shallow notch. This procedure was repeated such that 20–40 replicas were made over the life of each fatigue specimen. In the post-test crack evaluation procedure described by Jordon et al. (2012) each replica was coated with gold-palladium to improve the conductivity of the surface for scanning electron microscopy (SEM). However, they noted that the downside to using SEM to measure cracks was the tendency to melt the micron-sized features of the replicas under high magnification. As an alternate approach to using an SEM to measure the cracks, Bernard et al. (2013) demonstrated the efficacy of using an optical microscopy z-axis stitching method to measure cracks instead of the SEM approach. Similar, to Bernard et al. (2013), the surface replications in this present study were viewed and measured using a Keyence vhx-1000 optical microscope equipped with a Z axis imaging stacking feature that allowed for the entire surface of the replica to remain in focus. Crack measurements were made for one LCF and one HCF specimen at both considered thicknesses, for a total of four samples. Note that due to the nature of preparing each replica by hand, some of the

replicas were applied with insufficient pressure during the curing process resulting in poor impressions, thus rendering some images of the crack not usable. For surface replications with low quality images, the data points were removed from the analysis. Between 15–30 data points were available for each of these four specimens. The replica images were used to create simplified crack geometries for input into the SIF prediction code. Depending on the crack shape, the surface cracks were dimensioned and modeled as a single, straight crack consisting of two crack segments, a two segment crack with a single kink, or a three segment crack with two kinks. Examples of the idealized kinked crack geometries overlaid on surface replication images are shown in Figure 2. A high and low value were measured for each crack segment length to capture uncertainty due to image quality and measurement error. The total crack length was the sum of the lengths of all crack segments.

Experimentally observed crack growth rates were calculated using the crack length measurements obtained from each surface replication. A two-point secant method was used to calculate the fatigue crack growth rate (Eq. 1).

$$ \left(\frac{dc}{dN}\right)_n = \frac{c_{n+1} - c_{n-1}}{N_{n+1} - N_{n-1}} \quad [1] $$

where c is the total crack length, N is the cycle count, and n is the surface replication number. Note that because c is the total crack length, the computed growth rate is that of the entire crack, not an individual crack tip.

To characterize the microstructure and crack initiation, the fracture surface of each specimen was imaged using a scanning electron microscope (SEM). The crack initiation area of each specimen was examined, and the source of crack initiation was determined. Lastly, focus ion milling (FIB) was performed on the surface of the two LCF fractured specimens near the crack initiation site. Misorientation of grains was quantified at this location using EBSD. The observed fatigue crack’s position, crack length, and crack depth, d, were measured.

## Fatigue Crack Growth Rate Prediction Methodology

The computational approach developed and evaluated in this study consists of several steps to predict FCGR as described in this section. The DDT approach to calculate SIFs at the crack tips assumes a 2D infinite plate with a through crack and linear elastic behavior under remote stress. Therefore, the surface crack location in the bar and finite bar geometry must be measured and accounted for in another step in the analysis. Once $\Delta K_I$ and $\Delta K_{II}$ (representing mode I and II crack behavior respectively) are determined for each crack tip, a single equivalent stress intensity factor range ($\Delta K_{eq}$) is calculated according to Tanaka (Tanaka, 1974). A geometry factor is then applied to account for the surface crack depth. Using this corrected $\Delta K_{eq}$, the crack growth rate for the total crack growth is calculated via Paris law (Paris, 1961).
## Kinked Crack SIF Prediction

The computational model used to predict SIFs for microstructurally small crack growth is based on dislocation distribution theory that was originally developed for static loading (TerMaath et al., 2006; Yavuz et al., 2006) and adapted to fatigue loading (Yavuz et al., 2020). This model was chosen because it can solve for SIFs at the crack tips of kinked and branched cracks as well as combinations of the two. It also calculates the SIF at kinks by including the effects of wedge singularities so that crack growth at these locations can be evaluated. Importantly, it has been proven to be robust for cracks growing close to each other (TerMaath et al., 2006) and for large arrays of cracks of complex shapes (Yavuz et al., 2006), such that future development can capture this important scenario. Integrals are solved analytically (without numerical integration) eliminating this source of error. While not needed for this study, a significant benefit of this method is the capability to solve the stress and displacement field in a plate containing a large array of cracks of complex shapes, which can be highly interacting. This method is based on the isotropic bulk materials of the plate and assumes linear elastic behavior in an infinite body containing through cracks. A summary of this SIF prediction method is included in this section and full mathematical derivations are provided in (TerMaath, 2001).

The method is built on the superposition of two elastic problems, the trivial solution of the uncracked plate and the auxiliary problem that contains the cracks with prescribed tractions applied to the boundaries such the crack faces are traction free when the two problems are added (Eq. 2). Cracks are modeled as a series of crack segments of sufficient resolution to capture the crack configuration (as demonstrated in Figure 2). Thus, the foundation of the analysis is to solve for the crack opening displacement profiles for each crack segment that satisfy the global traction free condition. Once the opening displacement profile of each crack segment is determined, SIFs at crack tips and kinks, as well as the full stress and displacement fields in the plate, can be determined.

$$ \sigma_{xy}^{\infty} n_y + \sigma_{xx}^{\infty} n_x = -n_y \sum_{i=1}^{ncs} s_{xy}(Z_{1i}, Z_{2i}) - n_x \sum_{i=1}^{ncs} s_{xx}(Z_{1i}, Z_{2i}) \quad [2a] $$

$$ \sigma_{yy}^{\infty} n_y + \sigma_{xy}^{\infty} n_x = -n_y \sum_{i=1}^{ncs} s_{yy}(Z_{1i}, Z_{2i}) - n_x \sum_{i=1}^{ncs} s_{xy}(Z_{1i}, Z_{2i}) \quad [2b] $$

The variable $\sigma^{\infty}$ represents the far field loading conditions, $n_y$ and $n_x$ are the normals of the crack segment, $ncs$ is the total number of crack segments in the analysis, and $s$ denotes the tractions induced in a specific crack segment from the dislocation distributions. $Z$ is the Cauchy singular integral used to solve for local stresses (Eq. 3), and one and two denote the tangential and normal direction of the crack segment, respectively.

$$ Z_{1i} = \int_{0}^{a_i} \frac{\mu_i(t) dt}{z - t} \quad [3a] $$

$$ Z_{2i} = \int_{0}^{a_i} \frac{\mu_i(t) dt}{(z - t)^2} \quad [3b] $$

The variable $\mu$ denotes the dislocation distribution at a point due to an opening displacement profile, $i$ represents a specific crack segment, $z$ is the location of the traction free point in complex notation, and $t$ is the distance along the crack from the kink.

The traction free conditions are enforced at a user specified number of points along the crack segments and the point allocation is biased resulting in a higher density of points near the crack tips and kinks (TerMaath, 2001). The dislocation distributions are derived from opening displacement profiles that are constructed by summing terms from three types of critical crack behaviors (Eq. 4) to provide an overall crack shape with sufficient flexibility to capture the actual crack shape that meets the traction free boundary condition within a specified error. Eqs 4A–C are referred to as the Wedge, Polynomial, and Tip series, respectively, and each term in a series represents a degree of freedom. The total number of terms used for each series is designated as $nw$, $np$, and $nt$. Each series is used independently in the tangential and normal modes of deformation. Wedge series account for the stress singularity at a kink, polynomial series provide the opening and sliding crack displacements, and tip terms capture the singularity at a crack tip.

$$ W(t) = \sum_{j=0}^{nw-2} A_{pj} \left( \left( \frac{t}{a} \right)^{\rho+j} - (nw-j) \left( \frac{t}{a} \right)^{\rho+nw-1} + (nw-j-1) \left( \frac{t}{a} \right)^{\rho+nw} \right) \quad [4a] $$

$$ P(t) = \sum_{j=0}^{np-2} A_{0j} \left[ \left( \frac{t}{a} \right)^{j} - (np-j) \left( \frac{t}{a} \right)^{np-1} + (np-j-1) \left( \frac{t}{a} \right)^{np} \right] \quad [4b] $$

$$ T(t) = \sum_{j=0}^{nt-2} A_{\frac{1}{2}j} \left[ \left( \frac{a-t}{a} \right)^{\frac{2j+1}{2}} - (nt-j) \left( \frac{a-t}{a} \right)^{\frac{2(nf-1)+1}{2}} + (nt-j-1) \left( \frac{a-t}{a} \right)^{\frac{2nt+1}{2}} \right] \quad [4c] $$

The variable $A$ is the weighting coefficient of each term in the profile, and $\rho$ is a measure of the wedge distortion per kink. After the opening profiles are constructed, they are converted to dislocation distributions and are then integrated according to Eq. 3. The opening displacements are turned on and off at crack segment ends via the use of Dirac delta and Heaviside step functions. Non-trivial physical constraints are applied at adjoining crack segments to ensure material continuity.

The weighting coefficients from the opening displacement profiles are solved using a least squares fit to determine the contribution of each dislocation distribution to the overall problem solution. The weighting coefficients of the $j = 0$ tip series directly represent the singularity of the crack tip, and thus the mode I and II stress intensity factor are solved according to Eq. 5, where $G$ is the shear modulus of the material and $\kappa$ is
<figure><img src="image_234.png" /><figcaption>FIGURE 3 | Crack configuration for two crack segments and a single kink.</figcaption></figure>

Kosolov's constant which allows for plane strain or plane stress solutions.

$$ \begin{bmatrix} K_I \\ K_{II} \end{bmatrix} = \frac{2G}{1 + \kappa} \begin{bmatrix} A_{\frac{1}{2}0}^{[1]} \\ A_{\frac{1}{2}0}^{[2]} \end{bmatrix} \quad [5] $$

## Convergence and Uncertainty Quantification to Inform Kinked Crack Modeling

It is essential for the analytical model to provide converged and accurate results with minimized computation time. Therefore, a convergence study on the model form parameters and uncertainty quantification (UQ) on the parameters representing measurement error were performed using the SIF prediction method to inform code optimization. The crack configuration investigated is shown in Figure 3 for the case of a V shaped crack with a single kink. This configuration was chosen as it is representative of the crack growth patterns observed from the replications of the experimental specimens, and literature values for this crack configuration are available for validation. For this case, the crack is loaded in normal unit tension only. The baseline length of crack segment one is 0.6 with an angle $\theta_1$ equal to 60° and the length of crack segment two is 1 with an angle $\theta_2$ of 0. It should be noted that all parameters investigated in the current effort are considered as independent parameters.

A convergence study on the four model form parameters that are user-specified and defined as integer values consisted of the number of points assigned to each crack segment (the traction free condition is enforced at these points) and the number of wedge, tip and polynomial terms in approximating the opening displacement profile of each crack segment. Based on the results, the model form parameters were held at the following converged values for all analyses: 50 points were allocated along all crack

segments and six terms were designated for all opening displacement profiles.

Measurement error is introduced when obtaining crack lengths and angles from images and is dependent on both image quality and human error. This experimental error was defined by 4 parameters: the measured length and angle of each crack segment. Measurement parameters were defined as normal distributions using the mean and standard deviations provided in Table 1 and were sampled using a Latin hypercube sampling function. Latin hypercube sampling stratifies the individual variables of the input vector, sampling from these strata, and randomly pairing them to build samples of the full input vector (Helton and Davis, 2003). A 10% standard deviation was chosen for the measurement error of the lengths and the 1.75° for the angles based on the resolution of the images obtained from the specimen replications.

For automated sampling and generation of input files, we utilized Distribution-based Input for Computational Evaluations (DICE). DICE is a general-purpose tool used for the design of experimental parametric studies (Nahshon and Reynolds 2016). DICE receives a parameter file and replaces the parameter values in text files based on user-defined parameters for the study being performed. DICE outputs a series of ready-to-execute analysis input files. The four input parameters were defined in a parameter definition file and sampled as probabilistic parameters as previously described.

UQ results are provided in Table 1. According to the standard deviation values, some parameters contributed higher uncertainty than others to the results. However, all standard deviations are on the order of two decimal places in the accuracy of the SIFs. This error is well within current capability to measure the SIF of a material and it was deemed that the uncertainty due to measurement error was acceptable for this study.

## Straight Crack SIF Predictions

To generate straight crack results, the mode I and II SIFs were calculated according to Eq. 6 using the direct distance and incline between the crack tips and then converted to the effective SIF for use in the computational model.

$$ K_I = \sigma_y^\infty \sqrt{\pi a} (\cos^2\theta + \alpha \cos^2\theta) \quad [6a] $$

$$ K_{II} = \sigma_y^\infty \sqrt{\pi a} (1 - \alpha) (\cos \theta) (\sin \theta) \quad [6b] $$

where $\theta$ is the angle to the x-axis and $\sigma_x^\infty = \alpha \sigma_y^\infty$.

## Validation of SIF Predictions

While the computational code was rigorously validated in previous publications on this method (TerMaath, 2001), a few additional cases representing the microcracks under investigation were analyzed using the SIF prediction method and compared with known values. A single straight crack at an incline in an infinite plate crack was divided into two equal crack segments for varying crack lengths such that $c = 2a$. The crack was loaded in unit biaxial tension and compared to the theoretical solution given in Eq. 6 (Sih et al., 1962). For unit biaxial tension loading, $\sigma_y^\infty = 1$ and the loading variable $\alpha$ is equal to one so that the tension loading in both directions is equal to the unit value. The results of this comparison are provided in Table 2 and it is demonstrated the SIF predictions agree with the calculated values in all cases.
<table><caption>**TABLE 1** | UQ results for each input parameter relative to the output parameters.</caption><thead><tr><th rowspan="3">Varied input</th><th colspan="8">Output variation</th></tr><tr><th colspan="2">K<sub>I,1</sub></th><th colspan="2">K<sub>II,1</sub></th><th colspan="2">K<sub>I,2</sub></th><th colspan="2">K<sub>II,2</sub></th></tr><tr><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td>cs 1 length mean: 0.6 std: 0.06</td><td>0.8440</td><td>0.0079</td><td>0.9387</td><td>0.0115</td><td>1.4796</td><td>0.0143</td><td>0.0200</td><td>0.0079</td></tr><tr><td>cs 1 angle mean: 135° std: 1.75°</td><td>0.8442</td><td>0.0459</td><td>0.9379</td><td>0.0126</td><td>1.4796</td><td>0.0076</td><td>0.0200</td><td>0.0459</td></tr><tr><td>cs 2 length mean: 1.0 std: 0.1</td><td>0.8440</td><td>0.0302</td><td>0.9378</td><td>0.0303</td><td>1.4789</td><td>0.0527</td><td>0.0199</td><td>0.0302</td></tr><tr><td>cs 2 angle mean: 0.0° std: 1.75°</td><td>0.8436</td><td>0.0040</td><td>0.9390</td><td>0.0013</td><td>1.4783</td><td>0.0022</td><td>0.0199</td><td>0.0040</td></tr></tbody></table>

<table><caption>**TABLE 2** | Validation of SIF predicted values for an inclined straight crack SIF and a single kinked crack (specifically the crack tip of segment 1) with known values.</caption><thead><tr><th rowspan="2">Crack Configuration</th><th colspan="2">K<sub>I</sub></th><th colspan="2">K<sub>II</sub></th></tr><tr><th>Predicted</th><th>Known</th><th>Predicted</th><th>Known</th></tr></thead><tbody><tr><td>Inclined straight crack (Sih et al., 1962)</td><td></td><td></td><td></td><td></td></tr><tr><td>c = 2.010, θ = 0.49°</td><td>1.7769</td><td>1.7769</td><td>0</td><td>0</td></tr><tr><td>c = 2.107, θ = 4.79°</td><td>1.8194</td><td>1.8194</td><td>0</td><td>0</td></tr><tr><td>c = 2.227, θ = 8.95°</td><td>1.8704</td><td>1.8704</td><td>0</td><td>0</td></tr><tr><td>cs 1 length variation of a single kinked crack (Ishida and Nishino, 1982)</td><td></td><td></td><td></td><td></td></tr><tr><td>a = 0.02</td><td>1.1678</td><td>1.1747</td><td>0.6478</td><td>0.6537</td></tr><tr><td>a = 0.2</td><td>1.2857</td><td>1.2850</td><td>0.6459</td><td>0.6464</td></tr><tr><td>a = 0.4</td><td>1.4028</td><td>1.4034</td><td>0.6358</td><td>0.6324</td></tr></tbody></table>

<table><caption>**TABLE 3** | Comparison of literature values, DDT results, and FE outputs for serval crack configurations.</caption><thead><tr><th>case</th><th>DDT error (%)</th><th>DDT run time (seconds)</th><th>FE error (%)</th><th>FE run time (seconds)</th></tr></thead><tbody><tr><td>Straight crack</td><td>0.02</td><td>0.18</td><td>2.51</td><td>7.50</td></tr><tr><td>Symmetric crack with one kink</td><td>0.02</td><td>0.46</td><td>2.09</td><td>4.00</td></tr><tr><td>Symmetric crack with two kinks</td><td>0.05</td><td>0.81</td><td>3.49</td><td>4.10</td></tr></tbody></table>

A second case for unit biaxial loading was performed using the crack configuration in **Figure 3** for varying lengths of crack segment one which is held at a constant angle, $\theta_1 = 60°$ and for constant values of crack segment two of $a_2 = 2$ and $\theta_2 = 0°$. The comparison of the predicted SIFs with a known solution are provided in **Table 2** (Ishida and Nishino, 1982). Again, the values are a close match.

To compare the prediction code with another computational method, two dimensional finite element (FE) models of through cracks were created in Abaqus (2014) for several different crack configurations. These were used to validate the SIF prediction method for cases where literature values are not available and to provide a comparison in computational cost differences relative to FE. The FE models evaluated mode I and II SIFs via the contour integral method, with quarter tip quadratic elements surrounding the crack tips. To obtain relevant comparisons to the SIF predictive code, which assumes a crack in an infinite plate, cracks were created in the center of a square with large dimensions relative to the crack, with the edge length being at least 10x larger than the total crack length for all cases studied. The converged number of elements required for the FE models ranged from 50,058 to 95,728. **Table 3** shows the results for serval test cases, comparing SIF outputs and run times. In all cases studied, the predictive code was more accurate and required less

computation time. For the problems studied in this paper, the time reduction is not of significant importance, but for cases with hundreds or thousands of crack segments, the reduction in computation time will be considerable. When considering meshing a large crack array in FE or another high resolution method, the modeling and computational times can become prohibitive and a lower fidelity method is essential.

## Mixed Mode Calculation to Determine the Effective SIF

The mode I and mode II SIF values at each crack tip were combined into a single effective SIF, $K_{eq}$, to be input into a fatigue crack growth rate model. A wide range of methods to calculate $K_{eq}$ are available in literature. The method by Tanaka (1974), Eq. 7, was used due to its simplicity and applicability to the Paris model for crack growth (Qian and Fatemi 1996).

$$K_{eq} = \sqrt[4]{K_{I}^{4} + 8K_{II}^{4}} \quad [7]$$

## Edge and Surface Effects

The SIF predictive method used to compute K assumes through cracks in an infinite plate, whereas the experimental testing produced surface cracks in specimens with a finite bar configuration. Therefore, both surface and edge effects were
<figure><img src="image_235.png" /><figcaption>**FIGURE 4** | (A) Comparison of SIF results for varying crack location and plate size. Note that the solutions for the eccentric crack are piecewise due to the tabular look up format for some of the data. (B) Diagram of surface fatigue crack.</figcaption></figure>

evaluated to determine the effects of these two assumptions on the model and to account for each.

To determine the range of validity of the infinite plate assumption and to account for edge effects, literature values for straight, off-center cracks in a finite bar were compared to solutions for a crack of the same dimensions in an infinite plate. For short cracks located away from the edge, edge effects are minimal and the results for the infinite plate solution are in reasonable agreement with those for a finite bar. As the crack grows closer to the edge and the crack length is on the order of the dimensions of the bar, the differences become significant and the limit of validity of the SIF predictions must be determined. To determine this limit, the center of each observed crack was measured relative to the specimen edge on the fracture surface and input into the solution from a SIF handbook (Murakami, 1987). An error in the range of 5–10% between finite bar and infinite plate results was considered acceptable. Once the crack grew to a length that caused an error over the specified limit, the SIF results were considered invalid. A representative example of this calculation using specimen HCF T/4 is shown **Figure 4**.

SIF solutions for straight surface cracks are available in a variety of handbooks and publications (Murakami, 1987; Dowling, 2012). A simple solution to a surface crack problem is to apply a geometry factor, F, which is a function of the aspect ratio to the through crack solution to obtain the SIF at either the surface of the crack or the deepest point (Eq. 8). ΔK<sub>eff</sub> will be used to an effective stress intensity factor range after this surface crack correction factor has been applied. For example, a crack with an aspect ratio of 1 (the

ratio of crack depth, *d*, to the surface crack length, *c*) has a geometry factor of *F* = 0.73 that is applied to find the SIF at deepest point of the crack (Dowling 2012).

$$K_{\text{eff}} = F * K_{\text{eq}} \quad [8]$$

To determine the aspect ratio of the experimentally observed surface cracks, both literature estimates and aspect ratios on the final fracture surfaces were collected. Prior research suggests that crack aspect ratio remains consistent regardless of the initiating flaw geometry (Toribio et al., 2016) and is governed by the Paris law crack growth exponent, *m*, and loading configuration. The aspect ratios of the fatigue cracks on the final fracture surfaces were also measured, which ranged from 1.08 to 0.645 for the specimens tested. An average aspect ratio of 2d/c = 0.9 was applied to all experimentally investigated crack configurations. This value was chosen as representative to the majority of aspect ratios seen on the fracture surface, as well as its agreement to values seen in literature. A table with surface crack geometry factors for different crack aspect ratios is shown in **Table 4**. Using this average crack aspect ratio, a surface crack geometry factor of *F* = 0.705 was computed using the results presented in (Murakami, 1987), which accounts for crack depth and aspect ratio in a finite bar. The effects of using this average aspect ratio vs. the actual aspect ratio measured from the replications will be presented and reviewed in the discussion section.

## Paris Law FCGR Calculations

Once an effective SIF range was calculated, the values were input into a Paris Law predictive model (Paris, 1961) to calculate FCGR. The general form of this model relates ΔK<sub>eff</sub> to crack growth via an exponential relationship (Eq. 9). This model was chosen because material data to calculate the constants is commonly available for AA7075-T7351 as well as many other material types.

$$\frac{dc}{dN} = C(\Delta K_{\text{eff}})^m \quad [9]$$

Where dc/dN is the crack growth per loading cycle, ΔK is the stress intensity factor range, and C, *m* are constants fitted to experimental crack growth data. The C and *m* constants for the Paris law crack growth model were calculated using datasets from (Newman et al., 2005) and the Metallic Materials Properties Development and Standardization, MMPDS, (FAA, 2015) for 7075-T7351 tested at a load ratio of *R* = 0.1. A combination of data sets was used to capture the wide range of SIFs produced in

<table><caption>**TABLE 4** | Surface crack correction factors for various surface crack aspect ratios.</caption><thead><tr><td>2d/c</td><td>Correction factor</td></tr></thead><tbody><tr><td>0.2</td><td>0.2623</td></tr><tr><td>0.4</td><td>0.422</td></tr><tr><td>0.5</td><td>0.483</td></tr><tr><td>0.6</td><td>0.550</td></tr><tr><td>0.7</td><td>0.609</td></tr><tr><td>0.8</td><td>0.660</td></tr><tr><td>0.9</td><td>0.705</td></tr><tr><td>1</td><td>0.745</td></tr><tr><td>2</td><td>0.959</td></tr></tbody></table>