Florian Malard

Helo' Ise Verdier

Lara Konecny-Dupre1

David Eme

Christophe J Douady

Tristan Lefebure

ABSTRACT: Since the first discoveries in the 1860s of subterranean organisms in France, there has been a sustained effort in sampling and identifying subterranean species. Yet, the species inventory is still largely incomplete. Since the nineties, breakthrough in molecular methods has revolutionized concepts and practices guiding species inventory. Here, we provide a brief overview of how molecular methods have changed ways of delimiting, identifying and sampling species. Within the framework of integrative taxonomy, species taxa are now routinely delimited both on morphological and molecular grounds. Molecular delimitation methods often provide many more taxa with smaller range size than morphology. Species molecular inventory has simultaneously shifted from the sampling of specimens to that of environmental DNA (eDNA) and from the identification of a single organism (i.e. barcoding) to the simultaneous identification of a whole community (i.e. metabarcoding). Yet, there are still a number of challenges to be addressed for eDNA and metabarcoding to be applied routinely to the inventory of subterranean biodiversity. KEYWORDS: biodiversity, subterranean fauna, integrative taxonomy, DNA taxonomy, DNA barcoding, metabarcoding, environmental DNA.

RESUME: Depuis les premieres decouvertes dans les annees 1860 d'organismes souterrains en France, !'effort d'echantillonnage et d'identification des especes souterraines est reste soutenu. Pourtant, l'inventaire des especes est encore largement incomplet. Depuis les annees 1990, la percee des methodes moleculaires a revolutionne les concepts et les pratiques guidant l'inventaire des especes. lei, nous expliquons comment les methodes moleculaires ont change la fa�on de delimiter, d'identifier et d'echantillonner les especes. Dans le cadre de la taxinomie integrative, les especes sont desormais delimitees a la fois sur des bases morphologiques et moleculaires. Les methodes de delimitation moleculaire fournissent souvent beaucoup plus de taxons avec des aires de distribution plus restreintes que la morphologie. L'inventaire moleculaire des especes est simultanement passe de l'echantillonnage des specimens a celui de I' ADN environnemental (AD Ne) et de !'identification d'un seul organisme (barcoding) a !'identification en masse de toute une communaute (metabarcoding). Toutefois, ii existe encore un certain nombre de defis a relever pour que l'ADN environnemental et le metabarcoding soient appliques en routine a l'inventaire de la biodiversite souterraine.

Species inventory: a long tradition

In France, the first discoveries of subterranean organisms date back to the 1860s [START_REF] Ferreira | Niches within a niche: Ecological differentiation of subterranean amphipods across Europe's interstitial waters[END_REF]bibliographic synthesis in Ferreira, 2005]. Since then, sustained efforts have been made to sample and identify subterranean species, both aquatic [START_REF] Gibert | Environmental DNA in subterranean biology: range extension and taxonomic implications for Proteus[END_REF][START_REF] Zagmajster | Geographic variation in range size and beta diversity of groundwater crustaceans: insights from habitats with low thermal seasonality[END_REF] and terrestrial Uuberthie and Ginet, 1994;Faille, 2006;Mammola et al., 2019a]. Yet, the species inventory is still largely incomplete. Among the groundwater fauna, crustaceans have received much research attention because of their high diversity. The number of groundwater crustacean species in Europe exceeds that of surface-dwelling crustacean species [START_REF] Stoch | Stygobiotlc crustacean species richness: a question of numbers, a matter of scale[END_REF] (figure 1). Since the nineties, breakthrough in molecular methods has revolutionized concepts and practices guiding species inventory.

I. DNA taxonomy: species concept and delimitation of species taxa

We refer to DNA taxonomy as the delimitation and/or description of species using molecular methods. To unders tand the contribution of molecular methods to taxonomic thought, let's start by distinguishing between two separate meanings of the word "species": a concept and a taxon [de Queiroz, 2007]. The species concept accepted by most biologists treats species as evolutionary entities whichare temporally delimited by speciation events. Species taxa-for example the amphipod Niphargus rhenorhodanensis -are hypotheses of the evolutionary entities emerging during speciation. They are used by biologists as operational taxonomic units (OTU) to quantify biodiversity. To make these hypotheses, biologists use many different criteria -among which morphological and molecular criteriawhich evolve differently during speciation (figure 2). In literature, species hypotheses based on molecular criteria are often called OTU (or MOTU for molecular operational taxonomic units). However, note, that species hypotheses based on morphological criteria are also OTUs. Much of the debate around the "species problem" arose from the fact that numerous criteria were unnecessarily used to for mulate multiple definitions or concepts of the species [de Queiroz, 1998;[START_REF] Hey | On the failure of modern species concepts[END_REF]. Not surprisingly, h yp otheses may differ depending on the criteria used: these hypo theses are iteratively tested in a process called "integrative taxonomy" [START_REF] Yeates | -Integrative taxonomy, or iterative taxonomy?[END_REF]. Thatprocess has provided newspecies databases in which individual specimens can be assigned to different species taxa according to different criteria (figure 2). A key feature of such databases is that each morphologically distinguishable species taxon can be associated to one or several DNA b (l rcocles -a section of a DNA sequence from a specific gene. The DNA barcodes can repartition traduit en partie un biais d'echantillonnage lie a la localisation geographique des laboratoires de recherche (figure 1). Depuis les annees 1990, la percee des methodes moleculaires a revolutionne Jes concepts et !es pratiques guidant l'inventaire des especes.

I. Taxinomie moleculaire : concept d'espece et delimitation des taxons

Nous consicl erons ici que la taxinomie moleculaire correspond a la delimitation et/ ou a la description d'es peces a l'aide de methodes moleculaires. Pour comprendre la contribution des methodes moleculaires a la taxinomie, commern ;ons par cl istinguer cl eux sens distincts du mot « espece » : un concept et un taxon [de Queiroz, 2007]. le concept d'espece retenu par la grande majorite des biologistes traite Jes especes comme des entites evolutives qui sont temporellement cl elimitees par des evenements de speciation. Les taxons -par exemple l'amphipode Niplwrgus rhe11orl10da11e11sis -sont des hypotheses cl'entites evolutives emergeant au cours de la speciation. Ils sont utilises par !es biologistes comme des unites taxonomiques operation nelles (OTU) pour quantifier la biodiversite. Pour realiser ces hypotheses, !es biologistes utilisent de nombreux criteres -parmi lesquels des criteres morphologiques et moleculaires -qui evoluent cl ifferemment au cours de la speciation (figure 2). Dans la litterature, Jes hypotheses d'especes realisees sur la base de criteres moleculaires sont sou vent appelees OTU (ou MOTU pour molecular operational taxonomic units). Soulignons toutefois que Jes hypotheses d'especes realisees sur la base de criteres morphologiques sont aussi des OTUs. Les cl ebats sur ce qui a longtemps ete nomme « le probleme de l'espece » sont nes du fait que ces nombreux criteres ont eux-memes ete utilises, sans reelle necessite, pour formuler de multiples definitions ou concepts de l'espece [de Queiroz, 1998 ;[START_REF] Hey | On the failure of modern species concepts[END_REF]. Ainsi, les hypotheses peuvent differer selon Jes criteres utilises : elles sont testees de maniere iterative selon une methode appelee « taxinomie integrative» [START_REF] Yeates | -Integrative taxonomy, or iterative taxonomy?[END_REF]. Cette methode fournit des bases de donnees cl 'un nouveau genre cl ans lesquelles !es individus peuvent etre attribues a cl ifferents taxons en fonction des criteres utilises (figure 2). Une caracteristique de de ces bases de donnees est que chaque taxon morphologiquement distinct peut etre associe a un ou plusieurs code(s) -barre(s) ADN -une partie cl 'une sequence d'ADN cl'un gene specifique. Ces codes-barres peuvent ensuite etre utilises pour identifier ce taxon s'il est a nouveau collecte clans l'environnement ( voir section III ci-clessous).

Les methocles de delimitation moleculaire sont main tenant couramment utilisees pour delimiter Jes especes souterraines aquatiques et terrestres [Faille et al., 2015 ;[START_REF] Zagmajster | -Subterranean biodiver sity patterns from global to regional scales[END_REF][START_REF] Parimuchova | Hexapoda, Collembola) of the Western Carpathians Troglomorphy at the northern distributional limit in Europe[END_REF]. Erne et al.

[2018] ont utilise 2205 sequences du gene cytochrome be used, afterwards, to identify thetaxon if it is collected once again (see section Ill below).

Molecular delimitation methods are now commonly used to delimit aquatic and terrestrial subterranean spe cies [Faille et al., 2015;Erne et al., 2018;[START_REF] Parimuchova | Hexapoda, Collembola) of the Western Carpathians Troglomorphy at the northern distributional limit in Europe[END_REF]. [START_REF] Zagmajster | -Subterranean biodiver sity patterns from global to regional scales[END_REF] used 2205 sequences of the cytochrome c oxidase subunit I (COI) gene to delimit species of obligate groundwater Aselloidea (Jsopoda) and Niphargidae (Amphipoda) from Europe. These sequences came from 1075 sampling sites and collectively represented 263 morphologically distinguishable species (figure 3A). As is often the case, the authors found that molecular deli mitation provided many more species hypotheses (OTUs) in their study, 2.5 times as manythan morphological delimitation. Both delimitations revealed a ridge of high species richness at latitudes ranging from ca 4 2 ° to 46 °N, although the ridge was more pronounced when using molecular delimitation (figure 3 B&C). Moreover, the increase in range size of species at higher latitudes was stronger using COI-based species (figure 3 D&E).

II. An extending panel of molecular approaches and tools

DNA barcodingdefined herein as the molecular identification of individuals of already known species -has transformed the way species are identified and invento rized [Hebert et al., 2003] (figure 2). In subterranean environments, barcoding provided important information on many species for which morphological analysis was uninformative or ambiguous [START_REF] Zagmajster | Availability of DNA barcodes in subterranean amphipods in Europe[END_REF]. Since 2010, the emergence of new molecular technologies has extended barcoding in two ways: 1) a shift from specimens sampling to environmental samples (i.e. environmental DNA) and 2) a shift from the identification of a single oxydase I (COi) pour delimiter Jes especes souterraines d'Aselloidea (Jsopoda) et de Niphargidae (Amphipoda) en Europe. Ces sequences provenaient de 1075 stations d'echantillonnage et representaient collectivement 263 especes morphologiquement distinctes (figure 3A). Comme souvent, !es auteurs ont constate que la delimitation mole culaire fournissait plus d'hypotheses d'especes (OTU) -dans leur etude 2,5 fois plusque la delimitation mor phologique. Les deux delimitations ont revele une crete de plus forte richesse en especes a des latitudes s'etendant de 42 ° a 46 ° N. Toutefois, cette crete est plus clairement revelee par la delimitation moleculaire (figure 3 B&C). De plus, !'augmentation de la taille de l'aire de repartition des especes avec la latitude apparaH aussi plus clairement avec !es methodes de delimitation moleculaire (figure 3 D&E).

II. Un panel etendu d'approches et d'outils moleculaires

le barcodingdefini ici comme !'identification mole culaire d'individus d'especes deja connuesa transforme la fa <; on dont Jes especes sont identifiees et inventoriees [Hebert et al., 2003] (figure 2). En milieu souterrain, le barcoding a fourni des informations importantes sur de nombreuses especes pour lesquelles !'analyse morpholo gique etait non informative ou ambigue [START_REF] Zagmajster | Availability of DNA barcodes in subterranean amphipods in Europe[END_REF]. Depuis 2010, !'emergence de nouvelles technologies moleculaires a etendu le barcoding de deux manieres: 1) un passage de l'echantillonnage de specimens a des echantil lons environnementaux (l'ADN environnemental) et 2) un passage de !'identification d'un seul organisme (barcoding) a !'identification de toute une communaute (metabarco ding) (figure 2). l'ADN environnemental (ADNe)-defini comme le materiel genetique obtenu directement a partir d'echantillons environnementaux (par exemple de l'eau, du organism (i.e. barcoding) to the simultaneous identification of a whole community (i.e. metabarcoding) (figure 2). Environmental DNA (i.e eDNA)defined as genetic mate rial obtained directly from environmental samples (e.g. water, soil, sediment, etc.)has proved to be a sensitive approach to detect and identify rare and hard-to-catch species [START_REF] Thomsen | Environmental DNA-An emerging tool In conservation for monitoring past and present biodiversity[END_REF]. Several methods are used to collect eDNA from water samples, among which the ethanol precipitation method and the water filtration method. Water filtration, has been successfully applied to detect subterranean species such as Proteus anguinus [Goricki et al., 2017;Voros et al., 2017] and Stygobromus hayi [START_REF] Niemiller | Evaluation of eDNA for groundwater invertebrate detection and monitoring: a case study with endangered Stygobromus (Amphipoda: Crangonyctidae)[END_REF]. While eDNA is now used to inventorize soil fauna [START_REF] Kirse | Unearthing the poten tial of soil eDNA metabarcoding -Towards best practice advice for invertebrate biodi versity assessment[END_REF], it has not been applied to terrestrial cave fauna yet. The second extension from barcoding to metabarcodingemerged thanks to the capacity of parallel sequencing in second generation high-throughput sequencers. Instead of individually pro cessing and sequencing the DNA of each organism, this technology allows the simultaneous identification of mul tiple species from a large collection of specimens or even directly from environmental DNA [Critescu 2014]. eDNA metabarcoding has recently been applied to subterranean aquatic communities [START_REF] West | Under the karst: detecting hidden subterranean assemblages using eDNA metabarcoding in the caves of Christmas Island, Australia[END_REF].

sol, des sediments) -s'est avere etre une approche sensible pour detecter et identifier Jes especes rares et difficiles a capturer [START_REF] Thomsen | Environmental DNA-An emerging tool In conservation for monitoring past and present biodiversity[END_REF]. Plusieurs methodes sont utilisees pour recuperer l'ADN contenu dans l'eau des milieux, notamment la precipitation de cet ADN et la filtration d'un volume d'eau. Cette derniere methode a ete utilisee avec succes pour detecter des especes souterraines telles que Proteus anguinus [ Goricki et al., 2017;Voros et al., 2017] el Stygobromus hayi [Nierniller et al., 2018]. Alors que l'ADNe est desormais utilise pour inventorier la faune du sol [START_REF] Kirse | Unearthing the poten tial of soil eDNA metabarcoding -Towards best practice advice for invertebrate biodi versity assessment[END_REF], ii n'a pas encore ete applique a la faune terrestre cavernicole. la deuxieme extensionest apparue grace a la capacite de sequern ;age en parallele des sequenceurs a haut debit de deuxieme generation. Au lieu de traiter et de sequencer individuellement l'ADN de chaque organisme, cette technologie permet !'identification simultanee de plusieurs especes a partir d'une grande col lection de specimens ou meme directement a partir d'ADN environnernental [Critescu, 2014). le metabarcoding de l'ADNe a recemment ete applique aux communautes aquatiques souterraines [START_REF] West | Under the karst: detecting hidden subterranean assemblages using eDNA metabarcoding in the caves of Christmas Island, Australia[END_REF]. Bien que prometteuses pour l'inventaire de la diver site souterraine, !es approches moleculaires soulevent de nouveaux defis. Premierement en raison de la faible Although promising for the inventory of subter ranean diversity, these new molecular approaches raise new challenges. First, most approaches developed for surface habitat may turn out to be unsuitable to sample subterranean eDNA. For example, due to the low density of organisms in subterranean habitats, eDNA is likely to be in very low concentration, requiring much larger volumes of environmental samples, whose collection may happen tobeunfeasible in many situations. Second, one crucial element of a molecular identification analysis is the comprehensiveness of the reference database (figure 2). With many species or, even worse, entire groups not referenced, the identification will fail. Another pernicious effect of the lack of molecular reference data for many subterranean groups is that molecular methods are likely to be biased when applied to poorly known groups. Once collected, eDNA are extracted and then barcodes are amplified before being sequenced in parallel. Ateach step of the process, but particularly at the amplification one, protocols are by definition not optimized for poorly known groups. Therefore, the discovery of the latter may decrease and their mischaracterization be emphasized.

Ill. Beyond species inventory

Molecular approaches also foster the understanding of processes causing geographic variation in species richness. The number of species in a region ultimately depends on three processes: speciation, extinction and dispersal. Phylogenetic inferences made from DNA sequences on the evolutionary relationships among species can, for example, help to disentangle the relative influence of the three pro cesses in shaping the European ridge of high subterranean species richness (see figure 3) [START_REF] Zagmajster | -Subterranean biodiver sity patterns from global to regional scales[END_REF]. The occurrence within that ridge of many short-branch closely-related species may reflect a process of recent and rapid in-situ speciation. Conversely, the presence of species showing long independent evolutionary histories may reflect a low extinction rate which promotes the preservation of diversity over time. Finally, a high species richness may be due to species from neighboring regions dispersing into the ridge. The role of dispersal can be explored using methods such as the decomposition of phylobetacliversity indices [START_REF] Leprieur | Quantifying phylogenetic beta diversity: distinguishing between 'true' turnover of lineages and phylogenetic diver sity gradients[END_REF] and Bayesian phylogeographic diffu sion models [Eme et al., 2013]. The growing availability of molecular datasets over large spatial scales offers a unique opportunity to understand spatial variation in subterranean species diversity [Mammola et al., 20196]. It also paves the way for exploring factors causing the diversification [START_REF] Ribera | -Ancient origin of a Western Mediterranean radiation of subterra nean beetles[END_REF], as well as the morphological and molecular evolution of subterranean species [START_REF] Saclier | -Bedrock radioactivity influences the rate and spectrum of mutation[END_REF]Fiser et al., 2019]. sequences en parallele. A toutes ces etapes, mais particu lierement au niveau de !'amplification, Jes protocoles ne sont, par definition, pas optimises pour les groupes peu con nus, ce qui peut diminuer leur cl ecouverte et accentuer leur mauvaise caracterisation.

Ill. Au-dela de l'inventaire des especes

Les approches moleculaires permettent egalement de mieux apprehender les causes des variations geographiques du nombre cl 'especes. Le nombre cl'especes clans une region depend en fin de compte de trois processus : la speciation, !'extinction et la dispersion. Les inferences phylogenetiques realisees a partir des sequences d'ADN sur les relations de parente entre Jes especes pennettent de clemeler !'influence relative de ces trois processus sur la crete de forte richesse en especes observee en Europe (voir figure 3) [START_REF] Zagmajster | -Subterranean biodiver sity patterns from global to regional scales[END_REF]. Sur Jes arbres phylogenetiques, la presence a l'interieur de cette crete de nombreuses especes « jeunes » et etroitement apparentees peut refleter une forte speciation in situ. lnversement, la presence d'especes presentant de longues histoires evolutives independantes peut refleter un faible taux d'extinction qui favorise la preservation de la diversite au cours du temps. Enfin, une richesse specifique elevee peut aussi provenir de la dispersion au sein de la crete cl'especes en provenance des regions voisines. le role de'la dispersion peut etre explore a !'aide de methodes telles que la decomposition des indices de phylobetadiversite [leprieur et al., 2012] et des modeles bayesiens de diffusion [Eme et al., 2013]. la disponibilite croissante de donnees moleculaires a de grandes echelles spatiales offre une occasion unique de comprendre Jes causes des variations regionales de la diversite des especes souterraines [Mammola et al., 20196]. Elle ouvre aussi de nombreuses possibilites comme celles de dater plus precisement !es evenements de diversification en milieu sou terrain [START_REF] Ribera | -Ancient origin of a Western Mediterranean radiation of subterra nean beetles[END_REF] et de comprendre Jes mecanismes a l'origine de !'evolution morphologique et moleculaire en milieu sou terrain [Sacher et al., 2020;Fiser et al., 2019].

Remerciements

Nous remercions taus les collecteurs, panni lesquels de nombreux speleologues, qui contribuent a l'inventaire de la biodiversite souterraine. Cette etude a beneficie du soutien financier de l'Agence Nationale de la Recl1erche et de /'EUR . Nous remercions A. Faille et R. Centelles Bascuas pour la revision de eel article.

Figure 2 :KARSTOLOGIAFigure 3 :

 23 Figure 2 : Inventory of Subterranean biodiversity. 1. Sampling of specimens and environmental DNA from water and/ or sediment. 2. Delimitation of species taxa

  moleculaire de la richesse et de l'etendue des aires de distribution des crustaces aquatiques souterrains (Aselloidea, lsopoda et Niphargidae, Amphipoda) en Europe. A. Stations d'echantillonnage (points rouges) des crustaces en Europe. B & C. Richesses specifiques des crustaces par bande latitudinale de 0,9 ° obtenues par delimitations morphologiques (Bl et moleculaire (C) des especes. D & E. Etendues maximales de l'aire de distribution des crustaces obtenues par delimitations morphologique (D) et moleculaire (E) des especes. La delimitation moleculaire est effectuee par la methode du PTP (Poisson tree processes model) sur des sequences du gene COi. Les barres horizontales noires, les points rouges et les boites indiquent respectivement la mediane, la moyenne et l'intervalle interquartile pour les bandes latitudinales. Les lignes rouges continues et en pointille (intervalles de confiance a 95 %) representent l'ajustement d'un modele additif generalise aux moyennes des bandes latitudinales. Modifiee d'apres Erne et al. [2018].

F

  . MALARD, H. VERDIER et al., Subterranean species inventory in the molecular era/ lnventaire des especes souterraines dans l'ere moleculaire 76 KARSTOLOGIA n• 79, 2022 • 73-78 r

FKARSTOLOGIA

  . MALARD, H. VERDIER et al., Subterranean species inventory in the molecular era/ lnventaire des especes souterraines dans l'ere moleculaire 77 2014 -From barcoding single individuals to metabarcoding biological commu nities: towards an integrative approach to the study of global biodiversity. Trends in ecology & evolution, 29, 566-571. https://doi.org/10.1016/j. tree.2014.08.001 DE QUEIROZ K., 1998 -The general lineage concept of species, species criteria, and the process of speciation and terminological recom mendations. pp. 57-75, in D. J. Howard and S. H. Berlocher (Eds.), Endless forms: Species and speciation. Oxford University Press, Oxford, UK. DE QUEIROZ K., 2007 -Species concepts and species delimitation. Systematic Biology, 56, 879-886. https://doi. org/10.1080/10635150701701083 EME D., MALARD F., KONECNY-DUPR� L., LEF�BURE T. & DOUADY C. J., 2013 -Bayesian phylogeographic inferences reveal contrast ing colonization dynamics among European groundwater isopods. Molecular Ecology, 22, 5685-5699. https://doi.org/10.1111/mec.12520 EME D., ZAGMAJSTER M., DELIC T., Fl$ER C., FLOT J.-F., KONECNY-DUPR� L., PALSSON S., STOCH F., ZAK$EK V.

lnventai• re des especes souterraines dans l'ere moleculaire

des especes une longue tradition

  Figure 1 : Occurrence records (red dots, n=5311) of obligate groundwater crustacean species (n=295) in France. Black and purple patterns show quaternary deposits and karsts, respectively. Occurrences (points rouges, n=5311) d'especes de crustaces souterrains (n=295) en France. Les trames noires et violettes representent respectivement des depots quaternaires et des roches karstiques.

	En France, !es premieres decouvertes d'organismes sou
	terrains remontent aux annees 1860 [Ferreira et al., 2003;
	syn these bibliographique clans Ferreira, 2005]. Depuis, des
	efforts soutenus ont ete deployes pour echantillonner et
	identifier !es especes souterraines tant aquatiques I Gibert
	et Culver, 2009; Zagmajster et al., 2014] que terrestres
	Uuberthie et Ginet, 1994 ; Faille, 2006 ; Mammola et al.,
	2019a). Pourtant, l'inventaire des especes est encore lar
	gement incomplet. Panni la faune aquatique, les crustaces
	sou terrains ont beaucoup retenu !'attention des chercheurs
	en raison de leur grande diversite. Le nombre d'especes

. A total of 247 described species and subspecies of groundwater crustacean species are known forFrance and 48await description: species occurrences are more numerous in Southern and Eastern France but this distributional pattern partly reflects a sampling bias linked to the geographical location of research laboratories lnventaire de crustaces souterrains en Europe depasse celui des especes de crustaces vivant en surface

[START_REF] Stoch | Stygobiotlc crustacean species richness: a question of numbers, a matter of scale[END_REF]

. Au total, 247 especes et sous-especes de crusta ces souterrains sont connues en France et au moins 48 attendent d'etre decrites. Les occurrences d'especes sont plus nombreuses cl ans le Sud et l'Est de la France mais cette F. MALA RD, H. VERDIER et a/., Subterranean species inventory in the molecul ar era/ lnventaire des especes souterraines dans l'ere moleculaire

  , DOUADY C. J. & MALARD F., 2018 -Do cryptic species matter in macroecology? Sequencing European ground water crustaceans yields smaller ranges but does not challenge biodiversity determinants. Ecography, 41, 424-436. https://doi.org/10.1111/ ecog.02683 FAILLE A., 2006 -Endemisme et adaptation a la vie cavernicole chez /es Trechinae pyreneens (Co/eoptera : Caraboidea) : approches mo/ecu laire et morphometrique. PhD Thesis, Museum national d'histoire naturelle, Paris.

F. MALARD, H. VERDIER et al., Subterranean species inventory in the molecular era/ lnventaire des especes souterraines dans l'ere moleculaire KARSTOLOGIA n ° 79, 2022 • 73-78

F. MALARD, H. VERDIER et al., Subterranean species inventory in the molecular era / lnventaire des especes souterraines dans l'ere moleculaire KARSTOLOGIA n ° 79, 2022 • 73-78

Acknowledgments

We thank all the people, including many speleologists, who contribute to subterranean biodiversity inventory. Our aclmowledgement goes to the French National Research Agency and EUR for their financial support. We are also grateful to A. Faille and R. Centelles Bascuas for revising an early draft of that paper. densite d'organismes cl ans les habitats souterrains, l'ADNe est susceptible d'etre present en tres faible concentration. Ce qui necessitera la collecte de volumes cl 'echantillons environnementaux tres importants. Deuxiemement, !'ele ment cle cl'une approche d'iclentification moleculaire est l'exhaustivite de la base de donnees de reference (figure 2). lorsque de nombreuses especes, voire meme des groupes en tiers, ne _ sont pas references, !'identification peut echouer. le manque de clonnees moleculaires de reference pour de nombreux groupes souterrains peut egalement biaiser les methocles moleculaires. Une fois collectes, Jes ADNe sont extraits, puis les barcodes sont amplifies avant d'etre

Contents


Anderson S Nunes
email: anderson@esss.com.br

Olivier Chadebec

Gérard Meunier

Patrick Kuo-Peng

Magnetic devices analysis by Face FEM coupled with standard reluctance network

Keywords: Reluctance network, Finite Element Method (FEM), face shape functions, Face FEM

A Finite Element Method (FEM) mesh is converted to a reluctance network through an original magnetostatic formulation based on face shape functions. This meshed reluctance network is coupled with an standard one, characterizing a 0D system. Both approaches are fully-compatible and the hybridized problem can be solve with a single circuit solver. The approach is tested in 2D on a magnetic circuit with an air gap and compared to classical FEM nodal formulation.

I. INTRODUCTION

Optimizing electromagnetic devices can require a large amount of data that might be provided by numerical simulations. Many numerical methods are used to model electromagnetic devices, but the RNM (Reluctance Network Method) and FEM (Finite Element Method) are the most widely used for magnetostatic modeling.

The FEM is well known by its flexibility and generality, once the mathematical formulations are based on a mesh. Furthermore, it is noticeable the knowledge base available for this method, for instance [START_REF] Bastos | Magnetic Materials and 3D Finite Element Modeling[END_REF]. However, it leads to an high number of degrees of freedom and so quite long computation times.

On the other hand, the RNM is one of the most primitive methods for magnetic modeling and its application is based on a reluctance network. This method has remained useful due to its coherent results obtained with low computational effort and low computational simulation time [START_REF] Sykulski | Computational Magnetics, ser. International journal of numerical modelling[END_REF] and has been largely applied to model power transformers [3] [4]. This method is also largely applied to model rotating electrical machines [START_REF] Petrichenko | Contribution à la modélisation et à la conception optimale des turbo-alternateurs de faible puissance[END_REF][6] [START_REF] Tang | Analytical modeling of flux-switching machines using variable global reluctance networks[END_REF] and transmission lines [START_REF] Bormann | Reluctance network treatment of skin and proximity effects in multi-conductor transmission lines[END_REF]. Nevertheless, it is important to notice that these applications are based on a reluctance network defined manually, that might imply an hard, long and non-general task.

In [START_REF] Dular | Dual finite element formulations for lumped reluctances coupling[END_REF] is presented a methodology that couples nodal/edge FEM with external reluctances network. In this paper, we go a step forward by proposing a formulation fully-compatible with both numerical approaches and solved with a single 0D circuit solver.

Finally, the results of modeling a actuator with the classical nodal FEM and with the proposed methodology are compared.

II. MAGNETOSTATIC FACE FEM FORMULATION

The magnetostatic fields might be described by the Ampère's (1) and Gauss (2) laws and the constitutive relation for magnetic materials [START_REF] Lopez-Fernandez | Stray losses control in shell type transformers. part ii: Method of threedimensional network solution[END_REF].

curl H = J (1) div B=0 (2) B = µH (3)
where the magnetic field H, in A/m, is composed by the fields H 0 , due to an imposed current density source J 0 , and H m that can be obtained from the gradient of the magnetic scalar potential V r , as presented in [START_REF] Petrichenko | Contribution à la modélisation et à la conception optimale des turbo-alternateurs de faible puissance[END_REF].

H = H 0 + H m (4) H m = -∇V r (5)
Thus, applying (5) in ( 4) and integrating the resulting equation along a domain Ω, it is possible to obtain [START_REF] Fukuoka | A method for optimizing the design of spm type magnetic gear based on reluctance network analysis[END_REF], for which W is the face interpolation function.

Ω W i • H dΩ + Ω W i • ∇V r dΩ = Ω W i • H 0 dΩ (6)
Since the magnetic induction B, in T , is given by [START_REF] Tang | Analytical modeling of flux-switching machines using variable global reluctance networks[END_REF], the left side of ( 6) is rewritten as [START_REF] Bormann | Reluctance network treatment of skin and proximity effects in multi-conductor transmission lines[END_REF].

B = n f j=1 W i Φ i (7)
Ω W i • HdΩ = n f j=1 Ω υW i • W j dΩ Φ j (8)
where j and i are the faces index, nf is the number of faces, υ is the magnetic reluctivity, and Φ is the magnetic flux.

Applying the divergence theorem in the second term of (6) and evaluating it along two adjacent elements Ω a and Ω b , leads to [START_REF] Dular | Dual finite element formulations for lumped reluctances coupling[END_REF]. Taking into account that the normal component of the function W is constant along the face shared by Ω a and Ω b , its second term becomes null and it can be rewritten as (10).

Ωa+Ω b W i • ∇V r dΩ = Γa-Γ b V r W i • ndΓ- Ωa+Ω b V r ∇ • W i dΩ (9) Ω W i • ∇V r dΩ = - Ωa+Ω b V r ∇ • W i dΩ (10)
This equation can be split in terms of Ω a and Ω b , where the flux direction is defined from Ω a to Ω b . Then, considering that ∇ • w i is the inverse of the element volume, it is possible to obtain (11).

Ω W i • ∇V r dΩ = V r a -V r b (11)
The third term of ( 6) is the magnetic field due to an imposed current density J 0 , which can be obtained using Biot-Savart law, for instance.

Finally, lets rewrite ( 6) as the following matrix system

[ ] [Φ] -[∆V M ean ] = [H 0 ] (12)
where [ ] is a reluctance matrix,[Φ] is the unknown flux matrix, [∆V M ean ] is the magnetic potential jump between two reluctances and [V 0 ] contains the f mm sources. Besides the matrix system is stated, ( 2) is not solved yet. Nevertheless, the B free divergence is constrained solving (12) as a circuit system , where the Kirchhoff's current law is imposed. It suggests that the solution can be obtained by the use of a 0D circuit solver.

III. MODEL AND RESULTS

In order to compare the results obtained with the proposed methodology, a simple magnetic circuit with an air gap is analyzed using the well established nodal first order FEM, based on magnetic vector potential A. The magnetic induction distribution is presented in Fig. 1. Then, this model is reduced to a 0D model, i.e most of the magnetic circuit is modeled by an classical reluctance network and just the air gap region is meshed (i.e. the region where the equivalent reluctance network is not so easy to define), as presented in Fig. 2. Then the problem is solved by a circuit solver and magnetic induction is interpolated using [START_REF] Tang | Analytical modeling of flux-switching machines using variable global reluctance networks[END_REF], resulting in Fig. 3.

The maximum magnetic induction in the air gap obtained with the proposed methodology is 0.625 T and with the classical FEM is 0.606 T, that represents a difference of 3.13%.

IV. CONCLUSION

This paper has presented a methodology capable to convert a FEM mesh into a reluctance network, allowing its easy coupling with a classical network obtained analytically. The original model was reduced to a 0D system keeping coherent results.

Fig. 1 .

 1 Fig. 1. Magnetic induction distribution obtained with classical FEM

Fig. 2 .

 2 Fig. 2. Networks coupling


Douglas Martins Araujo

Diego Pereira Botelho

Yves Marechal

Jean-Louis Coulomb

Brahim Ramdane

Olivier Chadebec

Higher Order Natural Element Method -Boundary Element Method Coupling

Keywords: boundary element method, coupled numerical methods, natural element method

The advantage of a hybrid approach coupling natural element method (NEM) and boundary element method (BEM) is the conjunction of the inherent NEM accuracy and the BEM ability in modeling linear and deformable domains without any mesh. In a previous work, the coupling between standard NEM (1 st order consistency) and a 0 th order BEM was introduced. In the present work the performance of this coupling is tested for both methods in higher order versions. Results are presented in terms of convergence analysis and compared to a classic FEM-BEM coupling in terms of accuracy. The approach developed yields more accurate results with better convergence behavior, proving to be a good alternative in terms of accuracy for the simulation of unbounded problems.

I. INTRODUCTION

ir regions in computational electromagnetics are almost ubiquitous. In magnetostatics, for example, the accuracy of field quantities in air gaps is essential. In this context the Boundary Element Method (BEM) takes advantage over others methods, once it can model linear and deformable domains requiring only the discretization of theirs interfaces [START_REF] Liu | Fast multipole boundary element method: theory and applications in engineering[END_REF].

Others regions, like magnetic materials, must be modeled by other numerical methods. Traditionally Finite Element Method (FEM) is used in these domains. However, in this work, the Natural Element Method (NEM) [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] is employed. It has been previously shown that this method can provide highly accurate and computationally efficient solutions compared to FEM [START_REF] Maréchal | Computational Performances of Natural Element and Finite Element Methods[END_REF].

In a previous work [START_REF] Pereira Botelho | Hybrid Natural Element Method -Boundary Element Method for Unbounded Problems[END_REF] the coupling NEM-BEM was first introduced. In that paper, the scheme was based on a 1 st order consistent NEM and a 0 th order BEM. For the present work the consistency of both methods was improved to the 2 nd order. The performance of this improved scheme is tested and compared to the previous one. The tests presented in this abstract have been performed on a 2D magnetostatic problem.

II. SCHEME FOUNDATIONS

In this section we present briefly the main concepts involved on the implementation of the higher order NEM-BEM coupling.

A. Natural Element Method NEM standard shape functions are evaluated through geometrical relations based on a discretization scheme called Voronoï diagram [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. At the same time that these shape functions provide smooth interpolation they keep important properties, like partition of unity, positiveness, and strict interpolation [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. Among other things, these properties allow easy FEM-like implementation of resolution algorithms. However, as NEM shape functions present originally linear consistency, some additional technique must be applied in order to get higher order approximations. This issue was addressed in [START_REF] González | Higher-order natural element methods: Towards an isogeometric meshless method[END_REF], when a higher order NEM method based on the de Boor algorithm was developed. The main idea of the technique is the combination of different levels of linear interpolation that generates approximations with controllable consistency and precision. The method was first applied to the electromagnetic domain in [START_REF] Pereira Botelho | Higher Order NEM and FEM Accuracy Comparison[END_REF] with good results.

B. Boundary Element Method

If no source field in the free space is taken into account, the problem can be represented by a Laplacian equation in terms of a total scalar magnetic potential. The main idea is to solve Laplacian equation transforming the volume integral equation over the deformable domain into a surface integral over the boundary. The classic boundary integral equation is obtained by using the third Green's identity [START_REF] Liu | Fast multipole boundary element method: theory and applications in engineering[END_REF] and solved using the point collocation method. The scalar potential and the normal component of the magnetic induction are both interpolated using 2 nd order functions.

C. Coupling both methods

The results presented here are obtained through a reduced magnetic scalar potential. The coupling between NEM and BEM systems is taken into account by imposing the conservation of the normal component of the magnetic induction and the unicity of the potential on boundary. It is important to notice that the BEM part of the matrix is fully populated and the NEM part is sparse. By solving the coupled formulation, we get the reduced scalar magnetic potential on all nodes and the normal component of the magnetic induction on the boundary.

III. TEST CASE: U SHAPED ELECTROMAGNET

A. Model Fig. 1 shows a U shaped electromagnet. In order to test the proposed scheme the magnetic induction was computed on the A path 1, in the air. The reduced magnetic field was computed on the second path, inside the electromagnet. A. Interpolating the total magnetic field in the unbounded domain Fig. 2 presents the comparison between the new formulation and the 2 nd order FEM solution (reference) for the evaluation of the magnetic induction vertical. This reference was obtained through a mesh with more than 400,000 degrees of freedom (DoF). Fig. 2. Vertical component of the magnetic induction on path 1: comparison between FEM 2 nd order and the approach presented in this work. The number of DoF for each NEM-BEM simulation case can be consulted in the Fig. 3. Fig. 3 presents the relative error on the vertical component of the magnetic induction for both NEM-BEM A (results from the previous work [START_REF] Pereira Botelho | Hybrid Natural Element Method -Boundary Element Method for Unbounded Problems[END_REF]) and NEM-BEM B (results from this new formulation) approaches.

B. Interpolating the magnetic induction in the magnetic material

Fig. 4 presents the comparison between the tested methods in terms of the horizontal component of the reduced magnetic field. It can be seen that the new scheme is smoother and closer to the reference if it is compared with the previous NEM-BEM and FEM-BEM approaches addressed in [START_REF] Pereira Botelho | Hybrid Natural Element Method -Boundary Element Method for Unbounded Problems[END_REF]. IV. CONCLUSION Since coupled approaches using BEM are still currently being applied [START_REF] Qaseer | Hybrid Finite-Element-Boundary Element Analysis of a Single-Sided Sheet Rotor Linear Induction Motor[END_REF], we propose an improvement of these schemes. In this paper, a hybrid NEM-BEM 2 nd order formulation has been presented. It was able to provide accurate solutions even with very coarse discretizations. Regarding the next steps, we plan to compare the performance of this new formulation with the FEM using Whitney elements of 2 nd degree. We also intend to compare the results presented in this paper, using a reduced magnetic scalar potential, with those from the coupled reduced and total scalar potentials or t 0 -φ formulations [START_REF] Kladas | A new scalar potential formulation for 3-D magnetostatics necessitating no source field calculation[END_REF].

Fig. 1 .

 1 Fig. 1. U shaped electromagnet. Material with linear permeability and a nonmeshed coil. Two different paths are considered.

Fig. 3 .

 3 Fig. 3. Relative error (L2 norm) versus the number of DoF (three mesh densities: coarse, medium and dense).

Fig. 4 .

 4 Fig. 4. Horizontal component of the reduced magnetic field on path 2: comparison between the reference and different hybrid approaches.


Jonathan Siau

Gérard Meunier

Olivier Chadebec

Jean-Michel Guichon

Rémy Perrin-Bit

Volume Integral Formulation Using Facet Elements For Electromagnetic Problem Considering Conductors And Dielectrics

Keywords: Volume Integral Equation, Dielectric materials, Unstructured-PEEC, Magneto-harmonic

A new integral formulation is presented, enabling the computation of resistive, inductive and capacitive effects considering both conductors and dielectrics in the frequency-domain. The considered application here allows us to neglect any propagation effects and magnetic materials. In this paper, we will show how to improve the unstructured-PEEC (U-PEEC) approach to consider dielectric materials, keeping the same benefits. Results obtained with this formulation are compared to experimental data.

I. INTRODUCTION

Electromagnetic problems coupled with external circuit can be efficiently solved using the partial element equivalent circuit (PEEC) [START_REF] Ruehli | Equivalent Circuit Models for Three-Dimensional Multiconductor Systems[END_REF]. It is an integral equation based method, transforming a meshed electromagnetic device into a circuit of lumped elements R-L-C and sources. The classical PEEC is limited, because only structured mesh are supported. This limitation has recently been overcome by using facets element [START_REF] Nguyen | An Integral Formulation for the Computation of 3-D Eddy Current Using Facet Elements[END_REF] for general meshes. This major improvement enables the treatment of more complex geometries.

The PEEC is used for a wide range of frequency, from low-frequency application [START_REF] Gope | Solving Low-Frequency EM-CKT Problems Using the PEEC Method[END_REF] to telecommunication devices as antennas [START_REF] Freschi | Unstructured PEEC formulation by dual discretization[END_REF]. Most of these applications only consider conductive materials, and deals with different type of geometrical regions (volume or surface).

The contribution of this paper is to extends the PEEC method taking into account capacitive effects in presence of dielectric materials. In this paper, the dielectric material is associated to a volume mesh. We also have developed a variant of the formulation where the dielectric mesh can be limited to its boundary. This alternative will be presented in the full paper.

II. FORMULATION

Let us define the following regions: Ω J and Ω D contains respectively the conductors and the dielectrics, and Ω = Ω J ∪ Ω D . Their borders are defined as Γ J = ∂Ω J and Γ D = ∂Ω D .

Let us consider a problem with conductors, dielectrics and without any magnetic material. We have the following behavior laws

     J = σE in Ω J , D = E = 0 E + P in Ω D , B = µ 0 H
and according to Maxwell's equation we get

rotE = - ∂B ∂t = -µ 0 ∂H ∂t (1)
rotH = J + ∂D ∂t = J + ∂P ∂t + 0 ∂E ∂t (2)
and thanks to the Lorentz gauge, we have

A = µ 0 4π Ω 1 r J + ∂P ∂t dΩ (3)
∂V ∂t = 1 4π 0 Ω J + ∂P ∂t ∇ 1 r dΩ (4)
where A and V are the magnetic vector potential and the scalar electric potential respectively. The formulation limited to resistive and inductive conductors [START_REF] Nguyen | An Integral Formulation for the Computation of 3-D Eddy Current Using Facet Elements[END_REF] uses J as unknown. Its solenoidality was imposed with a circuit solver. In this new formulation considering dielectrics, a new unknown J t is choosen:

J t = J + ∂D ∂t = (σ + jω ) E (5)
and thanks to (2), we have divJ t = 0 in the whole problem. Then, we are discretizing J t using the same facet elements interpolation as in [START_REF] Nguyen | An Integral Formulation for the Computation of 3-D Eddy Current Using Facet Elements[END_REF] to overcome the structured mesh restriction.

The integral equation consists in matching electrical behavior law everywhere in the materials. From (3) and ( 4), we can write

E(P) = - ∂A ∂t (P) -∇V (P) (6)
and using the behavior law defined in (5), we have in the frequency domain:

J t (P) σ + jω = -jωA(P) -∇V (P). (7)
Thus, the use of a standard Galerkin projection in the frequency-domain leads to the following matrix system with [R t ] and [L t ] are the resistive sparse-matrix and the inductive dense-matrix respectively. The right-hand-side vector is the voltage of the dual mesh branches (see Fig. 1), where the dual mesh represents the equivalent electrical circuit. We can notice that V i is the averaged voltage on element i(or averaged voltage on border facets of Γ ext ) [START_REF] Nguyen | An Integral Formulation for the Computation of 3-D Eddy Current Using Facet Elements[END_REF]. I j is the current flowing trough the facet j.

([R t ] + jω[L t ]) {I} = {δV } (8) [R t ] i,j = Ω w i w j σ + jω dΩ (9) [L t ] i,j = jω Ω w i Ω σ + jω( -0 ) σ + jω w j r dΩ, (10)
The equation ( 8) corresponds to incomplete circuit equations. Actually, the capacitive effects are yet not taken into account, so we have to add some capacitive branches linking the border facets of Γ ext to a common node ∞, as shown on the Fig. 1. We suppose V ∞ = 0. So,

V Γext i = 1 S i Γext i V i dΓ
Using (4), we can express

V Γext = 1 4π 0 1 jω [P t ]{I} with [P t ] i,j = Γext i 1 S i Ω w j σ + jω( -0 ) σ + jω 0 ∇ 1 r dΩdΓ (11) Finally, the equation to solve is [R t ] + jω[L t ] + 1 jω [P t ] {I} = {δV }. (12)
In practise, we use a circuit solver to ensure the current conservation. This solver does not need any geometrical information, so the conservation is ensured for any geometry and discretization.

For high frequencies, we used an equivalent thin shell model for the conductors [START_REF] Nguyen | 3-D Integral Formulation Using Facet Elements for Thin Conductive Shells Coupled With an External Circuit[END_REF] to facilitate the consideration of the skin effect in the thickness.

III. NUMERICAL EXAMPLE

We have tested our formulation on a real device (see Fig. 2), made of two layers of copper 35µm-thick separated by a layer of dielectric FR4 with a thickness of 1.47mm (Fig 2).

The device is meshed using 1, 143 hexahedrons which leads to 3, 940 degree of freedom. Here, the mesh is light enough to compute the solution using a direct solver. The Lower Upper (LU-factorization) solution takes about 15 seconds on a Dell Precision M4800 with an Intel(R) Core(TM) i7-4800MQ CPU @2.70GHz and 32Go RAM. In addition, this formulation allows us to solve large scale problems using compression techniques, as the H-matrices [START_REF] Hackbusch | A sparse matrix arithmetic based on H-matrices. I. Introduction to H-matrices[END_REF], on the dense integral matrices [L t ] and [P t ].

The impedance measurements without any connection between conductors (opened-circuit) and simulation are compared on the Fig. 3.

We can see that the frequency of almost all peaks match pretty well, until high frequency effects appear. On the other hand, the amplitude of the peaks does not match exactly.

IV. CONCLUSION

We have presented a new formulation to consider conductors and dielectrics with a very few restrictive hypothesis. Thanks to the circuit solver, it can be applied to any geometry and external electrical circuit can easily be added.

Fig. 1 .

 1 Fig. 1. Example of a primal and dual mesh with capacitive external branches.

Fig. 2 .

 2 Fig. 2. Picture and dimensions of the studied device.

Fig. 3 .

 3 Fig. 3. Impedance modulus versus frequency. In blue: The reference (measurements), in red: the simulation results.


Lotfi Hedjazi

Sofiane Belhabib

Angélina D'orlando

Sofiane Guessasma
email: sofiane.guessasma@inrae.fr

Breaking material symmetry to control mechanical performance in 3D printed objects

Keywords: Additive manufacturing, 3D printed ceramic material, tensile performance, microstructure

This study aims at quantifying the effect of material arrangement in 3D printed ceramic on the mechanical performance. Starting from two main printing parameters, namely part orientation and printing angle, twelve different printing configurations are considered. These configurations differ in terms of filament arrangement in the building direction and within the plane of construction. Tensile testing is performed using vertical, lateral and longitudinal orientations for a printing angle of 0°, 15°, 30°, and 45°. Scanning Electron Microscopy is considered to study how the material symmetry affects the fractured patterns. In addition, optical imaging is used to monitor the deformation sequences up to the rupture point. The experimental results show a wide variety of deformation mechanisms that are triggered by the studied printing configurations. The interpretation of the observed trends in terms of mechanical load transfer is related to the lack of material connectivity and the relative orientation of the filaments with respect to the loading directions. This study concludes on the possibility to tune the tensile performance of 3D printed ceramic material by adjusting both the part orientation and the printing angle.

Introduction

Additive manufacturing is one of the breaking through technologies that has attracted many research areas for designing and manufacturing products with a large degree of freedom in design and architectural complexity. This processing route for material processing can be defined as a manufacturing technique layer by layer from a digitalised model. One of the distinctive features of AM is its ability to control the material laying down at a relatively fine scale. Several manufacturing routes fall within the definition of AM. One of the most known one is the fused deposition modelling, which is used to process polymeric structures. In this process, the material is extruded as a filament from a nearly 1.75 mm in diameter wire down to 0.1 -0.4 mm extruded filament. The motion of the printing nozzle that contains the extrusion space and the building platform on which the filament is laid down creates the 3D shape of the material by exploiting the rheological properties of the polymeric filament to quickly solidify after passing its thermal transition. This process generate a particular symmetry in the material due to the crossing of filaments at particular angles. The material is only continuous in one direction, which leaves two directions in which the material discontinuity materialises as a defect in cohesion between adjacent filaments called necking but also through the building direction. The presence of these discontinuities influences the mechanical performance of the parts. For instance, tensile properties measured in the direction parallel to the building direction are several orders of magnitudes lower than the same properties measured in the normal directions. The study of the sensitivity of process parameters on the mechanical performance of printed structures has received a large attention among the community of AM users. One of the main process parameters that contributes to tune the mechanical performance is the part orientation, which is found to significantly determine the amount of loss in properties depending on the loading pattern. In addition, the way the filaments are arranged within the building plane is also an influential factor.

Generally speaking, the filaments are organised in layups with a crossing configuration at -45°/+45°. These layups are complemented by an external frame to improve the structural stability of the printed part.

In this study, the combined influence of the part orientation and the filament crossing angle are studied to reveal the effect of symmetry break in material organisation at the filament scale on the mechanical performance of the printed part under tensile loading.

Experimental layout

The printing concept using FDM technology is illustrated in Figure 1a. Two main parameters determine the 3D filament arrangement within the part (Figure 1b). The first one is the inplane filament arrangement exemplified by the raster. The second one is the out-of-plane filament packing materialised by the building direction. The geometry of the part to be printed is a typical dogbone structure. Depending on the relative orientation of the dogbone structure with respect to the building direction, three main orientations can be selected: vertical (VE), longitudinal (LO) and lateral (LA). The filament arrangement within the plane of construction can be selected with the help of the printing angle parameter. Four different levels of printing angles are selected: 0°, 15°, 30° and 45°. These levels correspond to filament crossings at -45°/+45°, -30°/+60°, -15°, +75°, and 0°/+90°, respectively. By combining the part orientation and the printing angle, 9 different filament arrangements can be achieved, which are illustrated in Figure 2.

The printing process is conducted using a commercially available printer (Raise 3D Pro2 Plus) allowing high printing temperatures as large as 330°C. The feedstock material is a polylacticceramic filament (CER) from Frontierfila company. The recommended printing temperatures for this filament are above 200°C with a mandatory heated bed. The fixed printing conditions are shown in Table 1. As mentioned earlier, the varied printing parameters are the part orientation and the printing angle. Table 1 also summarises the levels considered in this study. The printed dogbone structures are subjected to tensile loading using a universal machine (Zwick Roell Group, Ulm, Germany). This machine is equipped with a load cell of 10kN. The rate of displacement is adjusted to 5 mm/min. The experiment is conducted up to the rupture point of the sample. The number of tested samples is 24 with 2 replicates per printing configuration. Out of the testing, the engineering constants such as Young's modulus, tensile strength, elongation at break are extracted and related to the printing conditions.

All experiments are monitored using optical observation. A high-speed camera (Phantom V7.3, Photonline, Marly Le Roi, 78-France) operating at a full frame of 800 × 600 pixels is used where the frame rate is adjusted between 100 and 5,000 fps (frames per second).

The fractured surfaces of the tested samples are characterised using scanning electron microscopy (SEM) to reveal the relationship between the symmetry of the material and the loading configuration.

Results and discussion

Microstructural features

Figure 2a illustrates the modification of the PLC filament shape from nearly circular to elliptical due to the strategy of laying down. Indeed, the layer height is adjusted to half of the nozzle diameter, which forces the layers to have a larger contact along the building direction.

The consequence is the change in the filament shape factor that affects the packing of the filaments, and in turn the generated porosity network. Within the filament, the ceramic particles are not affected by the increase of the temperature during the printing stage (Figure 2b). These appear as globular particles with a typical dimension of the order of 1 µm and a volume content of 15%. The examination of the interface between PLA and ceramic particles does not reveal particular interfacial damage, which puts more weight on the effect of filament arrangement in triggering the mechanical performance.

Figure 3 shows the out-of-plane layered structure corresponding to the vertical orientation, which is taken at an area far from the rupture region. Continuity in material structuring goes normal to the building direction (Figure 3a). The distance between adjacent filaments is nearly 200 µm, which also refers to the layer height. The central part of the filament appears smooth. However, the junction between adjacent filaments highlights a gap called necking. A close view of the necking (Figure 3b) reveals a textured white large structure of about 100 µm in width. This area emerges as a result of the deformation induced by the packing as the filament is forced to reduce its lateral dimension for 400 µm to 200 µm in the building direction. The cohesion between the filaments in this direction is only ensured by the limited contact area, which makes the 3D printed structure vulnerable to uniaxial loading along the same direction.

The in-plane arrangement of filaments for a sample printed according to the longitudinal orientation is shown in Figure 4. In order to highlight the filament crossing at +45°/-45° within the raster, the observed area is selected near the ruptured region where the stretching of the filaments allows both bottom and top layers to be viewed. Besides the limited connectivity between adjacent filaments within the same layer, filaments are cohesively connected between layers, where the contact area is substantially larger than for printed samples according to the vertical orientation.

In addition to the raster that is defined by the in-plane filament arrangement within the core of the sample, an external frame is added to confer more structural stability. Figure 4b shows the structure of this external frame along the edges in the case of a sample printed according to the vertical orientation. At the edges of the printed structure, the filament dimension is affected by both the compression of the top layers and the deceleration of the printing nozzle that occurs when trajectory changes drastically close to the edge. In addition, after the completion of the frame, the printing nozzle moves to the next position, which leaves a discontinuity as the one shown in Figure 4b.

Effect of part orientation

Figure 5a shows the effect of the vertical orientation on the fracture patterns of PLC printed samples. As the loading is performed in the same direction of material building, the load transfer is directed into a large area where the cohesion between filament is limited. As a consequence, this building configuration results in a limited mechanical strength and the instable cracking along the necking occurs substantially causing the ruin of the material (Figure 5b). The uniaxial loading is also transferred to the sample edge, which leads to the filament decohesion observed in Figure 5c. it can be concluded that for the vertical orientation, the structural role of the external frame is limited because the load bearing does not exploit the longitudinal performance of the filament.

Figure 6a shows the effect of the lateral orientation on the fracture patterns of PLC printed structures. In this case, the loading direction is no more aligned with the building direction. This means that the uniaxial stretching of the filament along its length represents most of the load transfer. Because the volume fraction of PLA matrix is predominant, the amount of stretching induces a severe plasticity behaviour, which is followed by a localisation phenomenon as shown in Figure 6a. The strain localisation induces significant reduction in filament cross-section prior rupture. A closer view at the junction between original and reduced sections of the filament highlights the presence of series of transverse cracks organised in bands (Figure 6b). The filament section reduction allows a large filament volume to be heterogeneously stretched along the length. This section reduction has also another consequence on the filament cohesion in the transverse direction. As shown in Figure 6b, a decohesion front develop in the direction opposite to the loading direction. Between adjacent filaments, the decohesion can be complete or partial, and can possibly leads to the formation and detachment of microfilaments as shown in Figure 6b.

Far from the localisation area, the overstretching of filaments involves mostly the filaments belonging to the external frame (Figure 6c). When the amount of stretching reaches the ultimate elongation at break of the filaments, the rupture pattern captured in Figure 6d eventually occurs. This pattern is typical of a ductile behaviour associated with the PLA matrix behaviour.

Figure 7a shows the last configuration among the possible part orientations, which is the longitudinal orientation. In this configuration, the filaments are stretched in a predominant +45°/-45° arrangement. This means that only a partial load transfer is obtained in the loading direction. This partial load transfer limits the stretching of the filaments compared to the previous filament arrangement (Figure 7a). A zoom-in view on the ruptured filaments shows the change of orientation of the filament observed in Figure 7b, which is a distinctive feature of the longitudinal orientation. In addition, out-of-plane decohesion between successive layers generates an inhomogeneous deformation for which the extent is guided by the width of the contact area. In facts, the filaments are overstretched between the contact points which act as anchoring points. Figure 7c provides a proof of the development of inhomogeneous stretching, for which the overstretched areas can be as large as 2 mm. A closer view of the overstretched filament does not show the presence of significant cracking bands. In addition, the transition between overstretched and understretched areas is abrupt.

Effect of printing angle

Depending on the part orientation, the effect of the printing angle can be amplified or inhibited. For a sample printed in the longitudinal direction, the in-plane arrangement of the filaments significantly influences the load transfer. The printing angle associated to the longitudinal orientation can alter significantly the magnitude of the load transfer as well as the nature of the deformation mechanisms. Figure 8 illustrates such influence for a part printed according to the longitudinal direction with a printing angle of 45°. The filament arrangement corresponding to this printing angle generates layups of 0°/90°. This arrangement allows half of the filaments to be oriented along the loading direction, which contributes significantly to the load bearing. It also generates a particular deformation mechanism that trigger interfilament cracking due to the retraction of the filament in the transverse direction. This retraction induces interfilament decohesion, which travels in the loading direction (Figure 8a).

This situation occurs at multiple spots driven by the lack of connectivity between the filament on top of each other. The longitudinal cracks that are created along the loading direction grow in number and possibly connect to each other as suggested in Figure 8b. Below these cracks, another mechanism of interfilament decohesion occurs in the transverse direction. However, the resulting crack propagates faster because the filament necking is normal to the loading direction. Figure 8c shows a closer view of the residual necking, which hosts a large number of microcracks that are formed by the reduction in filament section. According to this process, the crack bridging in the longitudinal direction adds more complexity to the deformation mechanisms induced by uniaxial loading.

Figure 9a shows the effect of an intermediate printing angle between 0° and 45° for a sample printed according to the longitudinal orientation. This figure exhibits partial detachment of the filament, which is induced by a large contact area between two successive filaments. The amount of damage observed is directly related to the extent of connectivity that is induced by the printing angle. In figure 9b, the combination of lateral orientation with a printing angle of 30° allows more filaments to be more aligned with the loading direction. As a consequence, the localisation phenomenon introduced earlier results in the partial detachment of filaments according to the same mechanism shown in Figure 6b. However, the larger load transfer triggered by a larger printing angle (30°) allows more the microfilaments to act transversely against longitudinal cracking. This is further enhanced when the printing angle is increased to 45° (Figure 9d).

Consequence on mechanical behaviour

Figure 10 illustrates the consequence of printing configurations on the overall tensile behaviour. Figure 10a demonstrates that irrespective of the printing angle, the deformation behaviour of printed structures according to the vertical orientation share the same cracking behaviour and the limited stretching prior rupture. It can be stated that when the building direction is parallel to the loading direction, the printing angle has no significant effect on the tensile response. However, this printing orientation leads in all cases to limited performance.

The opposite situation is shown in Figure 10b for samples printed according to the longitudinal orientation. In this case, the plane of construction contributes significantly to the overall deformation. Indeed, according to this orientation, the plane defining the filament arrangement contains the direction of loading. This means that the overall tensile behaviour depends on the printing angle as shown in Figure 10b. For instance, cracking behaviour is directly related to the printing angle. In addition, shearing develops for particular printing angles such as 0°, 15° and 30° whereas for 45°, uniaxial deformation prevails. Also, the large deformation observed in the case of longitudinal orientation compared to Figure 10a is directly related to the stronger contribution of the filament stretching in the load transfer. The most contrasting case is related to the lateral orientation (Figure 10c). This orientation allows the filaments belonging to the frame to be fully oriented in the loading direction. As a consequence, a large stretching prior ruptured is observed for all printing configurations. In addition, similarly to Figure 10b, there is a correlation between the printing angle and the amount of deformation that the printed structures can ensure prior rupture.

Table 2 summarises the quantified effects of both the printing angle and printing orientation on the mechanical parameters extracted from the mechanical responses. There is no cleat trend between the density of the printed structures and the printing angle. However, in average, it is found that the printed structures according to the lateral orientation exhebits the lowest densities (= 1.06 ± 0.009 g/cm 3 ). The same parts printed according to the longitudinal orientation exhibit the highest densities (= 1.09 ± 0.013 g/cm 3 ). This can be attributed to the amount of support needed to stabilise the structure during the printing process. The removal of this support can alter the density of the material. The intermediate case corresponds to the vertical orientation (= 1.07 ± 0.012 g/cm 3 ). Considerations about large discontinuities in the filament trajectories and the limited length laid down within the plane of construction can explain the lower density associated with this printing orientation.

Regarding the mechanical performance, although it appears that there is no strong trend with respect to the printing angle, the reported results in Table 2 suggest that higher performance is expected from a large printing angle. This is, for instance, true for Young's modulus at = 30° and  = 45°, irrespective of the printing orientation. In the case of the longitudinal orientation, the stiffness performance of the wire is even restored for these angles. In average, samples printed according to the longitudinal orientation are the top ranked ones (E Y = 730 ± 72 MPa). The lowest ranked ones correspond to the vertical orientation (E Y = 474 ± 78 MPa), which corresponds to a 40% of loss in the stiffness with respect to the wire performance. The analysis of the yield performance shows that the same rationale holds for the effect of the printing angle. However, the yield stress is found to be the largest one for the lateral orientation ( Y = 35 ± 2.32 MPa) whereas it decreases to  Y = 6 ± 1.4 MPa for the vertical orientation. For both cases, the yielding performance of the CER wire are not restored. The loss in the yielding performance represents 6% and 82% respectively. The tensile strength of the printed material follows the same tendency, where the gap between the yield stress and the tensile strength is limited for the longitudinal (

) and lateral ( ) orientation. The samples printed according to the vertical orientation do not present a plasticity behaviour. This means that the plasticity behaviour is significantly affected by the printing orientation. The best performing structures are those printed according to the lateral orientation ( Y = 36 ± 2.47 MPa) whereas the less performing ones are those printed according to the vertical orientation ( Y = 6 ± 1.41 MPa). This corresponds to a 28% and 87% of loss in the tensile strength with respect to the wire performance, respectively. The low values of ultimate stress ( R ) compared to the tensile strength observed for the lateral and longitudinal orientation demonstrates the development of localisation behaviour as shown in Figures 6 and7. There is no clear difference between the printing orientations and all configurations rank the same ( R = 7 ± 0.84 MPa). The ultimate stress is significantly reduced by 83% compared to the wire properties. The elongation at break follows the same trend as the tensile strength. It is significantly reduced for the vertical orientation ( R = 1 ± 0.12 %), which represents 98% of reduction with respect to the wire properties. The same elongation at break increases to its maximum level for the lateral orientation ( R = 29 ± 11.71 %), allowing to limit the performance loss to 62%.

Conclusions

This study demonstrates modulation of the mechanical performance of 3D printed ceramic materials can be obtained by combining both the part orientation and printing angle with respect to the loading direction. The SEM micrographs of fractured patterns concludes on varieties of deformation mechanisms that are activated and enhanced by playing on the relative orientation of the filaments with respect to the loading direction. These mechanisms are explained by the lack of material connectivity both in the building direction and within the plane of construction. The worst performing printing configurations are those printed according to the vertical orientation, in which the limited load transfer affects severely the stiffness, strength and ultimate performance. The lateral orientation appears to be the best compromise if the strength and ultimate properties are the main target. The maximum modulation of the tensile performance is found for the longitudinal orientation where it is possible to significantly change the resulting performance by adjusting the printing angle.

This configuration also fully restores the stiffness of the ceramic wire.

Figure captions

List of Tables

Figure 2 .Figure 4 .

 24 Figure 1. (a) Illustration of the fused deposition modelling technique, (b) filament arrangement

  Figure 5 (continued).

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Fractured surfaces of a sample printed according to the lateral orientation: (a) zoom-out at the stretched area, (b) a close view at the localisation area, (c) a global view of the overstretched filaments, (d) a zoom-in view on a ruptured filament.

Figure 9 .

 9 Figure 9. Effect of the printing angle on the deformation mechanisms : (a) longitudinal orientationprinting angle of 15°, (b) lateral orientationprinting angle of 30°, (c) lateral orientationprinting angle of 45°.