Benoît Piranda

Frédéric Lassabe

Julien Bourgeois

DisCo: A Multiagent 3D Coordinate System for Lattice Based Modular Self-Reconfigurable Robots

Localizing each module in a modular selfreconfigurable robot (MSR) is of paramount importance. In MSR, the communication graph is directly mapped to the real topology which makes the localization problem easy to solve. However, some types of connectors can lose the orientation of the modules, making the problem intractable. In this work, we propose to build a coordinate system for 3D lattice-based modular robots using a multiagent system. We present DisCo algorithm, that uses one agent per module which can only communicate with its connected neighbors and that does not need a central coordination system. We show that the agents can tackle any kinds of 3D lattice and we illustrate it with a Face Centered Cubic lattice (12 neighbors) and a cubic lattice (6 neighbors). Using communications and only four states, DisCo can also deduce the orientation of modules if the connectors do not provide this information.

I. INTRODUCTION

Building a robot from a set of interconnected modules has been an active research field for more than forty years now. Each module forming such a modular robot can also be enhanced by movement capabilities allowing reconfiguration of the modules in space. Such kind of systems have been called metamorphic robotic system and later on modular self-reconfigurable robots (MSR). Lattice-based MSR are one category of MSR in which the modules are aligned on a regular lattice. To ensure scalability, these systems use neighbor-to-neighbor communication such that adding one more module cannot overload the network capability compared to a bus. From a logical point of view, each module can be seen as an agent communicating with its neighbors and forming a multiagent system. An interesting characteristics of MSR is that the network connection which is a logical information, can also be interpreted as a physical information: the location of each neighboring module.

This location of the connected neighbors is a crucial information as it can be used to derive the shape of the ensemble of robots. Building a global coordinate reference system which can be applied to any kind of robots, consists in locating the position of all modules in the lattice and identifying their orientation.

If every module can have the relative location of its neighbors and their orientation, the problem of building the coordinate system is quite straightforward but it can be impossible if the orientation information is missing. Figure 1 shows a set of Blinky Blocks [START_REF] Kirby | Blinky blocks: A physical ensemble programming platform[END_REF], which are robots that can deduce the orientation of robots placed on the same This work was partially supported by the ANR (ANR-21-CE33-0009) and the EIPHI Graduate School (ANR-17-EURE-0002).

University of Franche-Comté, FEMTO-ST institute, CNRS, Montbéliard,

France [name].

[surname]@femto-st.fr plane but cannot get vertical orientation of their top and bottom neighbors. On the left picture, the reference robot is the rightmost Blinky Block and it is able to communicate the relative coordinates and orientations to all green robots. Yellow robots on the second floor and above can be placed but not oriented. The orientation and location of blue Blinky Blocks could not normally be retrieved as they are connected to yellow robots which do not know their orientation. However, there is an algorithmic solution to retrieve their position and orientation. The red Blinky Block can not be placed and oriented at all and it is therefore not possible to differentiate the 4 red robots connected to the top yellow one (cf. Figure 1b). This can cause a collision when moving without a sensor or can cause an incorrect representation of the global shape of the robot.

In this paper, we propose DisCo, a robust multi-agent architecture that builds a 3D coordinate system for latticebased modular robots. This algorithm uses one agent per module which is able to communicate only with its direct connected neighbors with no central coordination system. We demonstrate that the agents can tackle any kinds of 3D lattice even if the relative orientations of neighbors is not strongly defined. We show experiments with a face centered cubic lattice (the 3D Catoms [START_REF] Piranda | Designing a quasi-spherical module for a huge modular robot to create programmable matter[END_REF], 12 neighbors) and a cubic lattice (the Blinky Blocks [START_REF] Kirby | Blinky blocks: A physical ensemble programming platform[END_REF], 6 neighbors). Using communications and only four states, they can also deduce the orientation of modules if the connectors do not give this information.

In the remainder of this paper, we study existing localization algorithms for MSR and swarm. Second, we propose a global generic positioning model which aims at solving the positioning problem for any robot model (2D, 3D, mobile, static). Third, we estimate our model performances by implementing it on specific MSR platforms, through simulations with VisibleSim [START_REF] Thalamy | Behavioral simulations of lattice modular robots with visiblesim[END_REF], and through an implementation on real hardware, the Blinky Blocks.

II. RELATED WORK

The first thing to localize modules into an ensemble is to define a coordinate system which adapts to the different kind of lattice. There are many of them like Bravais lattices [START_REF] Bravais | Mémoire sur les systèmes formés par les points distribués régulièrement sur un plan ou dans l'espace[END_REF] but the most used in robotics is the Cartesian coordinate system augmented with the orientation information. In swarm robotics, modules are not aligned on a lattice and most of the time the space is a 2D space, simplifying the problem compared to a 3D one. Robots use wireless communications (infrared, ZigBee, etc.) which are used to evaluate the distance between them and build an approximated coordinate system. In [START_REF] Mclurkin | Algorithms for distributed sensor networks[END_REF], McLurkin defines the basic algorithms for sensor networks to determine physical positions from network topology: distance estimation, error estimation, neighbor counts and edge detection to cite a few. In the amorphous computing project, first ideas of a 2D coordinate system for sensor networks is found in [START_REF] Harold Abelson | Amorphous computing[END_REF], [START_REF] Coore | Establishing a coordinate system on an amorphous computer[END_REF], [START_REF] Nagpal | Organizing a global coordinate system from local information on an amorphous computer[END_REF] using communication propagated waves and three anchors and it is refined in [START_REF] Nagpal | Organizing a global coordinate system from local information on an ad hoc sensor network[END_REF]. Minimizing communications while having a sufficient accuracy is a complicated trade-off. In [START_REF] Klingner | Fault-tolerant covariance intersection for localizing robot swarms[END_REF], the authors propose to use co-variance intersection [START_REF] Julier | A non-divergent estimation algorithm in the presence of unknown correlations[END_REF] with moving robots.

Although there are similarities between swarm computing localization algorithms and modular robots coordinate systems, the main difference lies in the fact that the communication graph in modular robots is directly mapped to the real topology. Chirikjian proposed in [START_REF] Chirikjian | Evaluating efficiency of self-reconfiguration in a class of modular robots[END_REF] the first coordinate system for 2D hexagonal lattice without orientation. A different and more detailed method is proposed in [START_REF] Reshko | Localization techniques for synthetic reality[END_REF] to localize robots in 2D and 3D for hexagonal lattices but orientation works only for peripheral nodes as it needs movement from the robots to reach a consistent orientation. In [START_REF] Bhalla | A framework for analyzing and creating self-assembling systems[END_REF], the authors use a coordinate system in 3D without orientation to build self-assemble shapes. When the number of modules and the imprecision of lattice increase, a precise localization is difficult to achieve. In [START_REF] Funiak | Distributed localization of modular robot ensembles[END_REF], the authors propose an approximated 3D localization based on probabilistic evaluation of localization using a Maximum Likelihood Estimate (MLE) and a partition of the ensemble between dense areas and sparse ones. The results don't give a precise coordinate system but rather a rough estimate of the location, especially in sparse areas. In [START_REF] Hołobut | Internal localization algorithm based on relative positions for cubic-lattice modular-robotic ensembles[END_REF], the authors propose an algorithm for localizing modules in a 2D grid, without needing the orientation of the modules. This is the closest work to what we propose and the main differences are that they only consider a single topology and moreover only in 2D. Furthermore, the algorithm is synchronous, meaning that it needs a global clock or orchestration, and global halting condition has not been implemented.

III. THE DisCo MULTI-AGENT ARCHITECTURE

A. Theoretical model

In this work, we consider lattice-based MSRs with robots communicating with their connected neighbors using con-nectors disposed on the surface of the robot.

In a first approach, we consider that we can express the position and the orientation of the neighbor connectors relatively to the local coordinate system of the robot. Under this hypothesis, when a robot is connected to a neighbor, it can get the relative orientation of this neighbor. This constraint is very strong but can be solved by the latching system. We will show in Section III-B that it is possible to relax it using the network of robots.

In this article, we consider robots with 6 to 12 connectors in 3D -which is the maximum number that can be reached on spherical robots aligned in a Face Centered Cubic lattice (FCC) lattice. Specialized robots may for instance have 12 connectors for 3D Catoms [START_REF] Piranda | Designing a quasi-spherical module for a huge modular robot to create programmable matter[END_REF] or Datoms [START_REF] Piranda | Datom: A deformable modular robot for building self-reconfigurable programmable matter[END_REF], 6 connectors cubic Blinky Blocks [START_REF] Kirby | Blinky blocks: A physical ensemble programming platform[END_REF], Roombots [START_REF] Spröwitz | Roombots: A hardware perspective on 3d self-reconfiguration and locomotion with a homogeneous modular robot[END_REF], MBlocks [START_REF] Romanishin | 3d m-blocks: Self-reconfiguring robots capable of locomotion via pivoting in three dimensions[END_REF], Pebbles [START_REF] Gilpin | Robot pebbles: One centimeter modules for programmable matter through self-disassembly[END_REF] or Miche [START_REF] Gilpin | Miche: Modular shape formation by self-disassembly[END_REF], 6 connectors for 2D hexagonal robots (like mm catoms [START_REF] Naz | Network characterization of latticebased modular robots with neighbor-to-neighbor[END_REF]), 4 for Crystalline [START_REF] Rus | Crystalline robots: Selfreconfiguration with compressible unit modules[END_REF] or Smart blocks [START_REF] Benoit Piranda | A new concept of planar self-reconfigurable modular robot for conveying microparts[END_REF] in square lattice. Our method is also applicable to heterogeneous robots combining several shapes of robots in a same set. The only constraint is the position and orientation of each connector in the local coordinate system of the robot.

1) The robot model: Our goal is to deduce the (x i , y i , z i ) coordinates for every robot of a connected ensemble (i ∈ [0..n]). Although most of the robots are placed in a 3D environment, our model can process 2D environments by simply assuming the z coordinate to be 0. Only knowing its 3D coordinate is not sufficient enough for a robot to be part of a global coordinate system. It must also know its orientation defined by its angles related to the ( X, Y , Z) vectors of the global coordinate system. Then, each position and orientation of a robot k may be described by an homogeneous transformation matrix M k . For each robot, we consider a list of connectors C i (i ∈ [0..m]) which are placed and oriented in the local coordinate system of the robot applying a matrix W i .

2) Computing the position of a neighbor: If we consider a robot A, already placed in a referential, and connected to another robot B using the connectors C i from A and C j from B, then, the position and orientation of this second robot is given by the following relation:

M A • W i • S = M B • W j ( 1
)
Where S is the symmetry operator which must be applied to C i to get the orientation of C j . We can express M B from M A in the general case by:

M B = M A • W i • S • (W j ) -1 (2)
Equation 2 shows that robot B can deduce its position and orientation using the matrix M A • W i received from A. We can then deduce the following algorithm:

• First, elect the Coordinate System Origin (CSO), which set its local matrix M CSO to the identity and flood (M CSO • W i ) to all its direct neighbors. • Then, when a module B receives the matrix for the first time, it computes its position and orientation M B applying Equation 2, and sends the (M B • W i ) to all its neighbors but the sender.

B. Real case, an algorithm to relax constraints

When a robot is connected to a neighbor, the neighbor fill one of the cell of the associated lattice, the position of the neighbor can be deduced with confidence but its orientation around the connector is not always known.

In the following, we present DisCo, a full distributed algorithm to deduce the position and orientation of robots even if connection information are partial.

We first use the Blinky Block robot shape, Blinky Blocks are not fully symmetrical robots (see Figure 2 and1), these Lego like robots are placed vertically on the floor or over a neighbor, but rotations around the vertical axis cannot be detected by the hardware. As a consequence, we must consider two kinds of connectors: lateral connectors {East, N orth, W est, South} sharing the horizontal plane and the {top, bottom} connectors. For the lateral connectors, only one orientation of connection is possible (the two neighbors must be vertically aligned), then all robots placed in the same horizontal plane of an already placed and oriented robot are placed and oriented also as presented in section III-A .

But for neighbors placed on the top face (or the bottom face), there are 4 possible rotations around the vertical axis. Then a Matrices message to a neighbor placed on one of its {T op, Bottom} faces will embed 4 matrices.

On Blinky Block robots, the {East, N orth, W est, South} connectors C i (with i ∈ [0..3]) are placed at the center of the lateral faces. We express W i by the following homogeneous transformation matrices:

W i =         cos iπ 2 -sin iπ 2 0 1 2 cos iπ 2 sin iπ 2 cos iπ 2 0 1 2 sin iπ 2 0 0 1 0 0 0 0 1         (3)
And for the upper and lower faces:

W U pper =     1 0 0 0 0 1 0 0 0 0 1 1 2 0 0 0 1     (4)
and

W Lower =     1 0 0 0 0 1 0 0 0 0 -1 -1 2 0 0 0 1     (5)
For lateral connectors the symmetry matrix is unique and equal to:

S lateral =     -1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 1     (6)
But for top and bottom faces, the symmetry matrix may be one of the 4 following matrices S i with i ∈ [0..3]:

S i =     cos i π 2 -sin i π 2 0 0 sin i π 2 cos i π 2 0 0 0 0 -1 0 0 0 0 1     (7)
We can notice that in the case of Blinky Blocks, when several modules are placed one on top of the other, the combinations of matrices can be simplified to obtain only 4 different matrices. But of course, the number of stored matrices can become larger when modules are connected on the side faces of these Blinky Blocks.

1) Orientation Resolution algorithm: The model starts by choosing the CSO, which can be arbitrarily set (based on an identifier value for instance), or elected with some constraints in terms of connectivity. This CSO's orientation is set to be the global orientation.

Each agent placed on each robot stores a list of candidate matrices encoding potentials positions and orientations. When a robot floods geometrical data, it sends the list of all possible matrices to a connected neighbor. When the neighbor receives this list, it merges it with its local current list in order to add new possible matrices or to reduce the list by confirming possibilities and/or by deleting others. Each agent can take four different states, associated with the size of the list of matrices.

• UNKNOWN is the initial state of the agents except the CSO which is in PLACED AND ORIENTED state, it means the agent has never received any matrix. This state corresponds to an empty list of matrices. • UNTRUSTED: More than one possible position have been received. This situation appears when the list contains more than one matrix and the matrices do not correspond to the same position.

• PLACED ONLY: The list of matrices contains many matrices but only the orientation differs.

• PLACED AND ORIENTED: The list of matrices is reduced to only one matrix which gives a sure position and orientation of the robot. As soon as an agent receives a Matrices Message it goes from UNKNOWN to one of the three other states depending on the size of the sent list of matrices.

In order to change their state, agents must receive a message that embeds a list of matrices. Merging these matrices with their local matrix reduces the number of possible matrices. The merge function compares the current list of matrices and the one received in the message in order to retain only compliant ones. If the receiving robot has an empty list, it copies all matrices from the source. Else, it retains only the intersection of both matrices lists.

2) Messages and agent attributes: Each agent stores the list of matrices in myMatrices. In our implementation, to reduce the memory size of embedded variables, we compress each matrix in 7 bytes (6 bytes for the coordinates and one for the orientation, coded by the number of the connector aligned with -→ x ).

3) Breadth first search algorithm: The main problem of the algorithm is its termination, indeed an agent cannot know when the process of matrices exchange is finished. Only agents in PLACED AND ORIENTED state have finished the process but all the others may wait for a new matrix forever.

To tackle this problem, we propose to use a Breadth First Search (BFS) that activates step by step the diffusion of the matrices to all children with the same distance to the root module. When the farthest modules have finished their diffusion and return an acknowledgment to the CSO, the positioning process is finished. The DisCo algorithm requires three main steps and 4 agent states as shown in the right side of Figure 3. Initially, all agents are in the "Waiting" state:

1) The first treatment is the election of the CSO, the elected module changes to "CSO" state and the others go to the "Gradient" one.

2) The second step consists in creating a spanning tree from the CSO, storing at each node the list of children (myChildren attribute) and the distance to the CSO (myDistance attribute), which corresponds to the depth of the node in the tree. At the end of the process, modules send an acknowledgment to the CSO, and change to the "Solving" state. Steps 1 and 2 may be merged into one step when using an election algorithm such as ABC-Center [START_REF] Naz | Abc-center: Approximate-center election in modular robots[END_REF]. 3) During the third step, the CSO agent sends its matrices to all its children and activates all the nodes at the same distance from itself and wait for an answer. It activates all modules at distance 1 first, then modules at distance 2, and so on. When activated, an agent sends its matrices to all its neighbors. These matrices are merged with the others in order to reduce the list of possible positions. If the number of matrices is reduced, it re-sends the new list to its neighbors. 4) Finally, the CSO broadcasts a message such that each module will have its final state: "PLACED AND ORI-ENTED ", "PLACED ONLY " or "UNTRUSTED ". Each module is then ready to begin the execution of another application which would need the coordinates and orientations of the modules.

For more details, this algorithm uses 6 different messages presented in the right part of the Figure 3. Depth and AckDepth are used for building the spanning tree and inform the CSO at the end of this step. Solve and AckSolve are needed for activating nodes and for informing the CSO at the end of the activation. The message Matrices is sent by activated nodes to exchange matrices between neighbor modules. Finally, the Finish message is in charge of the termination of the algorithm.

4) Generalization to 3D Catoms geometry and connections: The 3D Catom is a quasi-spherical module placed in a FCC lattice and connected to up to 12 neighbors. Considering two 3D Catoms, latched using two connectors, there are two possible orientations called up and down as shown in Figure 5. The IDs of a connector are written in black for the up orientation and in red for the down orientation, the axis represent the local coordinate frames. Figure 5b shows the two possible orientations of a latched 3D Catom for the same couple of connectors. The gray 3D Catom (left) is latched using its connector #1, and the right 3D Catom uses the same connector labeled #6 on the top (yellow) and #18 on the bottom (blue) depending on its orientation, up (top, yellow) and down (bottom, blue). As a consequence, each 3D Catom placed in the FCC lattice can be oriented in 24 different ways. To implement the computation of neighbor position on 3D Catoms, we define the twelve W i matrices (i ∈ [0..11]) deduced from the coordinate systems of the connectors C i . And we consider two symmetry matrices S up and S down which corresponds respectively to a rotation of 180°around the -→ z and the -→ y .

We implement DisCo in VisibleSim simulator to show the different steps of the DisCo Figure 4a presents our model made of 45,189 3D Catoms randomly oriented, and the position of the elected CSO. In Figure 4b, the colors represent the distance to the CSO after the gradient step. With this configuration all modules are well placed and oriented by the DisCo algorithm as it can be seen from the green color of the Figure 4c. 5) Time complexity of the algorithm: As the computation time is negligible in regards of the communication time, the complexity of the "Gradient" and the "Solving" steps is expressed in function of the time of sending and receiving a message between two neighbors, named t c .

The "Gradient" step consists in constructing a spanning tree from the CSO module whose complexity is O(d × t c ) [START_REF] Lynch | Distributed Algorithms[END_REF]. The "Solving" step consists in activating sequentially all the modules at a distance r = k, with k ∈ [1.. ]. The complexity of this step is therefore

O (2 × r=1 r) = O( 2 ).
Whatever the position of the CSO in the configuration, we have ≤ d where d is the ensemble diameter. Therefore the time complexity of the whole algorithm is

O((d + 2 )t c ) ≈ O(d 2 × t c ).
This result can be checked in Figure 7, drawing the simulation time in number of messages produced during the simulation of the DisCo algorithm depending on the eccentricity. It shows two graphs for two different models made of a wide number of modules (up to 45,000) and an eccentricity from 10 to 111. The two curves are aligned, which shows that the time of the simulation does not depend on the shape (nor the number of modules) of the configurations but on the eccentricity of the CSO only.

IV. EXPERIMENTS AND RESULTS

The algorithm has been implemented both in real robots Blinky Blocks and in the VisibleSim simulator. This last tool allows to scale up the number of modules. In Figures 7 and8, we use VisibleSim to compare the duration times for two different shapes: a cube and the LegoGuy presented Figure 4. These two objects being vectorial (CSG object as presented in [START_REF] Thadeu | Efficient scene encoding for programmable matter self-reconfiguration algorithms[END_REF]), it is easy to increase the number of modules in resizing the objects. We simulate 25 configurations made of up to 46,656 3D Catoms.

We observe in Figure 8 that the shape of the configuration has an effect on the simulation time for the same number of modules. The cube being a dense object, other shapes will have a larger diameter for the same number of modules, which will produce a longer time of computation.

About memory used in modules during the algorithm, it needs a fixed number of local variables (myDistance, myParent, myChildren, etc.), plus the list of matrices. The matrices are encoded in modules using 4 integer (3 for coordinates and one for the orientation). The number of matrices depends on the size of "one module large" loops. When the ensemble is dense, for example in the shape of Figure 4, the maximum number of matrices stored in a agent is 4. But in the disadvantageous case presented Figure 6 The video 1 shows the algorithm in action for the two implementations. First in VisibleSim, we can see the propagation of several messages in large sets of modules (up to 45,000). Then, we present many different executions the algorithm implemented on 17 to 450 Blinky Blocks. We also show that the algorithm reacts to dynamical updates of the configuration by adding and removing modules to a set of modules that already have their coordinates.

V. CONCLUSION

When a latching connector in a lattice based MSR does not give the orientation of the connected module, setting up a coordinate system is not trivial. In this paper, we propose a multiagent system which uses communication to solve this loss of orientation in the general case of 3D lattice based MSR, no matter the kind of lattice being used. The complexity of the algorithm is in O(d 2 ), d being the diameter of 1 Youtube video link: https://youtu.be/WH lP3kgpB8 the system. We test this algorithm in two different contexts, a simulator, to show the scalability and a real modular robot, to prove our algorithm works in real conditions. We show our method scales up to 45,000 robots in simulation and using 450 real Blinky Blocks, we validated the dynamical use of DisCo. Indeed, Blinky Blocks are connected/disconnected during the execution of the algorithm without causing any issue.

We envision several future works to enhance DisCo and different applications that would benefit from a coordinate system on MSR. First, we observed that it is possible to enhance the algorithm by adding a validation process that uses a larger neighborhood. It would reduce the number of matrices by checking if a possible orientation of a robot induces collisions or connections that does not exist. Then, the choice of the CSO could be enhanced by choosing a central node, however, there is a trade-off to be studied in the benefit of running a center election algorithm, which will also use resources. Many uses stem from having a consistent coordinate system over a MSR, like for example, Fault detection and ID assignment. Adjacent robots with faulty network connections can detect the faults by combining information from the coordinate system and from local communications and a consistent coordinate system can concatenate x,y,z coordinates into a unique software identifier.

Fig. 1 .

 1 Fig. 1. Illustration of the issue of the coordinates definition on Blinky Blocks robots.(Better printed in color)

Fig. 2 .

 2 Fig. 2. Two pictures of experiments with Blinky Blocks extracted from the video.(Better printed in color)

Fig. 3 .

 3 Fig.3. On the left: the global graph of agents executed in the robots during the algorithm with transition events. On the right, the DisCo algorithm detailed by agent.

Fig. 4 .

 4 Fig. 4. 3 steps of the algorithm illustrated on the LegoGuy model (45,189 modules) using VisibleSim [26]. a) Election of the CSO module. b) Gradient from the CSO, the colors represent the distance to the CSO. c) Merging of matrices. (Better printed in color)

Fig. 5 .

 5 Fig. 5. a) Zoom on a 3D Catom connector with its IDs: black for the 'up' orientation and red for the 'down' orientation. b) Two 3D Catoms linked by the same connectors #1 and #6. On the top, right 3D Catoms (in yellow) is oriented 'up', on the bottom, right 3D Catoms (in blue) is oriented 'down'.

Fig. 6 .

 6 Fig. 6. Worst case configuration of 3D Catoms. Green modules are well placed and well oriented. Orange one are well placed but many orientations are possible. The position and orientation of the red module cannot be deduced by the algorithm.(Better printed in color)

Fig. 7 .

 7 Fig. 7. Duration of the simulation (expressed in function of tc) depending on the eccentricity of the CSO for two shapes: a cube and the Lego Guy model.(Better printed in color)

Fig. 8 .

 8 Fig. 8. Duration of the simulation (expressed in tc) in function of the number of modules for two shapes: a cube and the Lego Guy model.


Wei Guo

Yu-Ke Ma

Yan-Feng Wang

Vincent Laude

Yue-Sheng Wang

Dual-tunable phononic waveguides for manipulation of guided Lamb waves

In this paper, we design and fabricate dual tunable waveguides in a twodimensional periodic plate with threaded holes. Dual tunability is realized by using rods held with nuts as well as assembly prestress of the nuts. A straight waveguide, a bent waveguide and a wave splitter are designed by changing the distribution of rods and nuts in different circuits. The experimental and numerical results show that the frequencies of guided waves can be tuned by the assembly prestress. By increasing the amount of prestress, the frequency range of the passing band can be shifted upward. Confinements, guiding and splitting of Lamb waves are clearly observed in both experimental measurements and numerical simulations. This work is essential for the practical design of reconfigurable phononic devices.

I. INTRODUCTION

 . Bragg bandgaps are mainly dependent on the structural periodicity. The corresponding wavelength is in the same order with the lattice constant. While in contrast resonant bandgaps are mainly dependent on the resonance of the scatterers. The wavelength inside the bandgap is much larger than the lattice constant. The nature of bandgaps can be used for vibration reduction 5-7 , and sound isolation 8-10 .

When one or some scatterers are removed or replaced with different shapes or materials, the periodicity of PCs is broken. Defected states are then induced inside the bandgap [START_REF] Sigalas | Elastic wave band gaps and defect states in two-dimensional composites[END_REF][START_REF] Sigalas | Defect states of acoustic waves in a two-dimensional lattice of solid cylinders[END_REF][START_REF] Li | Acoustic confinement and waveguiding in two-dimensional phononic crystals with material defect states[END_REF] and can be used to design novel wave devices. Point, linear and planar defects can be obtained by modifying one, one line or one layer of scatterers. Wave or energy can be trapped around the point defect, so point defects might be applied for the design of novel cavity [START_REF] Torres | Sonic band gaps in finite elastic media: Surface states and localization phenomenain linear and point defects[END_REF][START_REF] Khelif | Trapping and guiding of acoustic waves by defect modes in a full-band-gap ultrasonic crystal[END_REF] or energy harvester [START_REF] Jo | Designing a phononic crystal with a defect for energy localization and harvesting: Supercell size and defect location[END_REF][START_REF] Wu | Acoustic energy harvesting using resonant cavity of a sonic crystal[END_REF] . A linear defect can be used to confine the energy around the defect and guide the propagation of energy. So they can be used to design waveguides [START_REF] Khelif | Guiding and bending of acoustic wavesin highly confined phononic crystal waveguides[END_REF][START_REF] Li | Bending and branching of acoustic waves in two-dimensional phononic crystals with linear defects[END_REF] . This can also be realized by using a linear chain of point defects [START_REF] Wang | Channeled spectrum in the transmission of phononic crystal waveguides[END_REF] . Planar defect introduces modes with energy localized in the plane [START_REF] I E Psarobas | Phononic crystals with planar defects[END_REF] . Combination of different defects can be used to design various wave devices, such as filter, coupler, splitter and frequency demultiplexer [START_REF] Pennec | Tunable filtering and demultiplexing in phononic crystals with hollow cylinders[END_REF][START_REF] Pennec | Acoustic channel drop tunneling in a phononic crystal[END_REF][START_REF] Wang | Guid-ing and splitting lamb waves in coupled-resonator elastic waveguides[END_REF] .

Despite the existence of the above wave behaviors, practical applications of PCs are rarely demonstrated because the operating frequency ranges of the PC-based devices are generally fixed. The concept of tunable PCs has thus been developed as a separate research field [START_REF] Wang | Tunable and active phononic crystals and metamaterials[END_REF] . Generally, there are two main regulation mechanisms. On the one side, tunable PCs are realized by physically adjusting the material parameters. Such systems generally consist of multiphysics coupled media, which are sensitive to the electric field [START_REF] Li | A self-adaptivemetamaterial beam with digitally controlled resonators for subwavelengthbroadband flexural wave attenuation[END_REF][START_REF] Ning | Broadband square cloak in elastic wave metamaterial plate with active control[END_REF][START_REF] Zhang | Tunability of band gaps in two-dimensional phononic crystals with magnetorheological and electrorheological composites[END_REF] , magnetic field [START_REF] Yang | Significant tuning of band structures of magneto-mechanical phononic crystals using extraordinarily small magnetic fields[END_REF][START_REF] Gao | Investigation on tunable low-frequency property of magnetic field induced phononic crystal with archimedean spiralbeams[END_REF][START_REF] Gu | Research on the tunability of point defect modes in a two-dimensional magneto-elastic phononic crystal[END_REF] , optical field [START_REF] Baumgartl | Tailoring of phononic band structures in colloidal crystals[END_REF] and thermal field [START_REF] K C Chuang | A tunable elastic metamaterialbeam with flat-curved shape memory alloy resonators[END_REF][START_REF] Yang | Tunable band gap and wave guiding in periodic grid structures with thermal sensitive materials[END_REF][START_REF] Liu | Thermally tunable topological edge states for in-plane bulk waves in solid phononic crystals[END_REF] . On the other hand, tunable PCs are realized by mechanically changing the geometrical topology. In addition to filling fluid in PCs with blind holes [START_REF] Wang | Tunable fluid-filled phononic metastrip[END_REF][START_REF] Wang | Wave propagation in coupledresonator acoustoelastic waveguides[END_REF] , prestress is commonly used to tune wave propagation. Apart from granular PCs [START_REF] Kim | Direct measurement of superdiffusive energy transport in disordered granular chains[END_REF][START_REF] Allein | Tunable magneto-granular phononic crystals[END_REF] , the related investigations on the continuous prestressed PCs generally include two parts.

First, some investigations focus on PCs with traditional elastic materials. Gei et al investigated flexural wave propagation in the prestressed periodic beam on an elastic foundation [START_REF] Gei | Band-gap shift and defect-induced annihilation in prestressed elastic structures[END_REF] .

They showed that bandgaps can be shifted upward (downward) for tensile (compressive) prestress. Wang et al. calculated the effects of prestress on the bandgaps of three-dimensional pizeoelectric phononic crystals [START_REF] Wang | Band gaps of elastic waves in three-dimensional piezoelectric phononic crystals with initial stress[END_REF] . The location of bandgaps was found to increase with increasing prestress, while the width of bandgap remains almost unchanged. Zhu et al studied the bandgap adjustability under prestretch strain [START_REF] Zhu | Tunable characteristics of low-frequency bandgaps in two-dimensional multivibrator phononic crystal plates under prestrain[END_REF] . Applying prestretch strain to the matrix can realize active realtime control of low-frequency bandgap under slight deformation and broaden the low-frequency bandgap. Li et al showed that the frequency ranges of nonreciprocal transmission can be shifted down for elastic waves propagating in the combination of nonlinear material and elastic metamaterials with increasing prestress [START_REF] Li | Tunable nonreciprocal transmissionin nonlinear elastic wave metamaterial by initial stresses[END_REF] . Large stresses, however, are required for significant tunability of the elastic materials [START_REF] Feng | Tuning the band-gap of phononic crystals with an initial stress[END_REF] .

Second, some investigations considered PCs with hyperelastic materials. Compared to their elastic counterparts, hyperelastic materials generate large deformation even under a small stress. Bertoldi's group has investigated wave propagation in several soft porous PCs under prestress [START_REF] Bertoldi | Mechanically triggered transformations of phononic band gaps in periodic elastomeric structures[END_REF][START_REF] Bertoldi | Wave propagation and instabilities in monolithic and periodically structured elastomericmaterials undergoing large deformations[END_REF] . Novel configuration patterns based on the instability and bulking provides rich room to tune the wave behaviors [START_REF] Shan | Harnessing multiple folding mechanisms in soft periodic structuresfor tunable control of elastic waves[END_REF] . Sharma et al. introduced the optimum topologies of the soft compressible unit cells for varying levels of the applied prestretch as well as the longitudinal prestress [START_REF] Kumar Sharma | Topology optimization of soft compressible phononic laminates for widening the mechanically tunable band gaps[END_REF] . They observed that increasing the lateral prestretch together with applying a compressional longitudinal prestress have favorable impact on widening the bandgaps. Bayat and Gordaninejad formed wrinkles in 1D prestressed PC slab composed of a thin soft film attached to a thick substrate [START_REF] Bayat | Switching band-gaps of aphononic crystal slab by surface instability[END_REF] found a remarkable shift of dispersion curves for inertial amplified PCs compared to Bragg or resonant ones [START_REF] Miniaci | Effect of prestress on phononic band gaps induced by inertial amplification[END_REF] . More related work may be found in the recent review [START_REF] Bertoldi | Flexible mechanical metamaterials[END_REF] .

Although there is a lot of literature on prestressed PCs, investigations on the defected PCs are rarely reported [START_REF] Wang | Tunable guided waves in a soft phononiccrystal with a line defect[END_REF] , especially for Lamb waves. Manipulation of Lamb waves by PCs combining prestress and defects remains a difficult task. In this paper, we design and fabricate dual tunable waveguides in a two-dimensional periodic plate with threaded holes. Dual tunability is realized by using rods held with nuts as well as assembly prestress of the nuts.

Guided waves can be tuned by adusting the distribution of rods and nuts in different circuits and controling the amount of assembly prestress. Straight waveguide, bent waveguide and wave splitter are designed by changing the distribution of rods and nuts in different circuits. The frequencies of guided waves are tuned by the assembly prestress. Confinement, guiding and splitting of Lamb waves are clearly observed in both experimental measurements and numerical simulations. This work is important for the design of reconfigurable or programmable elastic wave devices. Numerical simulations are conducted by using the finite element software COMSOL. The assembly prestress is taken into account through prestrain in the Solid Mechanics Module.

Before solving for the band structures or transmission properties of the structures, a sta-tionary study is performed to calculate the static deformation of the structure under the prestress. The geometrical and material parameters of the plate, rods and nuts are exactly the same with those in Ref. [START_REF] Wang | Guided lamb waves in reconfigurable phononic crystal waveguides[END_REF] . Fig. 1(b) illustrates the supercell of the waveguide used for the calculation of band structures. The plate is divided into 4 parts in the thickness direction. Quadratic serendipity element is applied to approximate the displacement field.

Periodic Bloch bondary conditions are applied along the x and y directions, and the other surfaces are set as traction-free. The band structure is obtained by sweeping the wave vector along the irreducible Brillouin zone of the supercell. Eigenmodes at the marked points are also obtained by choosing the suitable wavenumber and frequency. In addition, the polarization amounts for the out-of-plane modes are calculated based on the eigenmodes. They represent the contribution of out-of-plane displacements to the total polarization of elastic waves.

Transmission properties are evaluated by considering a 3D plate model. An asymmetric wave source is applied to the top of the plate. The response is collected in the region S 1

(with the area A 1 ) on the end of the waveguides. We define the numerical transmission as

F = 20log 10 ( S 1 |U 1 |/A 1 dS |U 0 | ), (1)
where U 0 = 1 is the amplitude of the z-polarized wave source and U 1 is the total displacement received on the end of the waveguides (S 1 ).

III. RESULTS AND DISCUSSIONS

In this section, we discuss the band structures and transmission properties of dual tunable waveguides with different prestress.

A. Defect modes of prestressed waveguides

First we investigate the band structures and defect modes of prestressed waveguides. The supercell is shown in Fig. 1(a). For comparison, we also consider the phononic waveguides assembling by rods and nuts with perfect bonds. Band structure of the waveguide with perfect bond is presented in Fig. 2(a), together with the vibration modes of the defected bands in the marked points. It is observed that the band structure is different from that in [START_REF] Wang | Guided lamb waves in reconfigurable phononic crystal waveguides[END_REF] because two nuts are added on both sides of the plate. This also indicates the reconfigurability of the proposed waveguide. The defect mode at point P1 is flexural-dominated but asymmetric with respect the wave propagation direction. So it is a deaf mode and can not be excited by a symmetric plane wave. This also holds for the defect mode at point P3.

So only the defect mode at point P2 contributes to the transmission. ilar to the case with perfect bond. The defect band around 65 kHz is almost independent of the prestrain because the threaded rod feels almost no vibrations, see the corresponding vibration modes at points P4 and P7. The vibration is asymmetric with respect to the wave propagation direction, so these two bands are deaf. This also holds for points P6 and P9.

Only the guided modes P5 and P8 contribute to the transmission.

It is noted that the adjustment of prestress has a range. The stress of the rods cannot exceed the yield limit. and 6.3Ö10 8 N/m 2 , respectively. They are all within the yield limit of the material, which is 6.9Ö10 8 N/m 257-59 .

The above analysis shows that the defect bands appearing in the band gaps are sensitive to the prestrains. Waveguiding at different frequencies thus can be tuned by changing either the prestrains applied on the nuts or the geometric pattern of the defects.

B. Straight waveguide

In this subsection, we focus on wave propagation in straight waveguides with a length to the coupling of the modes, the causes of frequency shift are complicated. In addition, there are other differences between numerical results and experimental measurements. When screwing the nuts in the experiments, it cannot be guaranteed that the force exerted by the nut on the rod is uniformly acting on the force area of the rod, which is not accounted for in the numerical simulation. Besides, the thread structures are not considered in the simulation. The observed differences may be attributed to all these contributions.

C. Bent waveguide

In addition to the straight waveguide, tunability is also feasible for the design of other waveguides, such as waveguides with 90

• bend. The bent waveguide is composed of 7 defected units. Transmission of the bent waveguide with perfect bond is given in Fig. 6a. The results with different prestrains are given in Figs. 6b and6c. It is found that there is a passing band above 65 kHz for the waveguide with perfect bond, the frequency range of which is identical to that of the straight waveguide in Fig. 4a. With the increase of prestain, passing band in the transmission spectrum shift upwards. phenomenon is also observed around 55 kHz, identical to the experimental results in Fig. 5 for the straight waveguide.

D. Wave splitter

In addition, we also fabricated a wave splitter in the plate with the aid of reconfigurability. Figs. 8 and9 give the simulated and experimental transmissions and displacement This paper is relevant to the practical application of elastic acoustic devices. With the aid of control elements, the geometric pattern and applied prestress between nuts and plates can be tuned in a programmable way. Active and even smart control of elastic waves are expected.

  . The bandgaps were adjusted by surface instability. Chen et al. realized a tunable soft acoustic diode in the bandgap by changing the applied force 50 . De Pascalis et al. employed genetic algorithm to optimize prestressed phononic media for the low frequency bandgaps of antiplane elastic waves 51 . Miniaci et al.

Figure 1 :

 1 Figure 1: Dual tunable phononic splitter. (a) Photograph of the wave splitter and the Digital Torque Wrench used to evaluate the assembly prestress. (b) Supercell of waveguide used to calculate the band structure

Figure 2 :

 2 Figure 2: Band structures and vibration modes of the phononic waveguides with perfect bonds (a) or different prestrains: (b) 0.0083 and (c) 0.0167. The color scale represents the polarized amount from 0 (blue) to 1 (red). The light-gray areas indicate the passing band for the out-of-plane polarized waves in the perfect PC slab. The dark-gray areas indicate the considered frequency range of guided bands. Vibration modes at marked points are shown on the right. The color scale represents the normalized amplitude of out-of-plane displacements from -1 (blue) to 1 (red).

Figs. 2

 2 Figs.2(b) and (c) show the band structures and vibration modes of the phononic waveguides with different prestrains. The two coupled defect bands shift upward with increasing prestrain. It seems, however, that the band structures and the bandgaps do not change significantly when prestrain is applied. This can be quantitatively identified by observing the eigenfrequency of the marked points. The vibration modes of the defect bands are sim-

Figure 3 :

 3 Figure 3: Von Mises stress distribution of the supercells with different prestrains: (a) 0.0083 and (b) 0.0167. The structures around the rods are hidden. The color scale represents the stress amount from 0 (green) to 6.3Ö10 8 (red).

Figs. 3

 3 (a) and Figs. 3(b) show the von Mises stress distribution of the structures with different prestrains. It is found that the stress is mainly concentrated on the rods. Under the two prestrains, the maximum von Mises stress in the rods are 3.4Ö10 8

of 7 Figure 4 :

 74 Figure 4: Simulated transmissions of the straight waveguides with perfect bonds (a) or different prestrains: (b) 0.0083 and (c) 0.0167. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to 1 (red).

Figure 5 :

 5 Figure 5: Measured transmissions of the straight waveguides with different torques: (a) 2.3 N • m, (b) 2.7 N • m, and (c) 3.0 N • m. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to maximum (red).

Figure 6 :

 6 Figure 6: Simulated transmission properties of the bent waveguides (7 units) with perfect bonds (a) or different prestrains: (b) 0.0083 and (c) 0.0167. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to 1 (red).

Figure 7 :

 7 Figure 7: Measured transmissions of the bent waveguides (7 units) with different torques: (a) 2.3 N • m, (b) 2.7 N • m, and (c) 3.0 N • m. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to maximum (red).

Figure 8 :

 8 Figure 8: Simulated transmission properties of the wave splitter with perfect bonds (a) or different prestrains: (b) 0.0083 and (c) 0.0167. Transmissions collected on the right and left outputs are marked by blue and red, respectively. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to 1 (red).

Figure 9 :

 9 Figure 9: Measured transmissions of the wave splitter with different torques: (a) 2.3 N • m, (b) 2.7 N • m, and (c) 3.0 N • m. Transmissions measured on the right and left outputs are marked by blue and red, respectively. Vibration modes at the corresponding labeled points are given in panel (d). The color scale represents the normalized amplitude of out-of-plane displacements from 0 (blue) to maximum (red).

Acknowledgments

The authors wish to thank Dr. Shuai Yang for his help in the numerical simulation. Financial supports by National Key R&D Program of China (2022YFB3806101), National Natural Science Foundation of China (12122207, 12021002, and 11991032) and the EIPHI phononic crystal circuits formed by coupled acoustoelastic resonators. Physical Review Applied, 8:014006, 2017.

Data Availability Statement

Data in this article is available upon request.

Conflicts of Interest

The authors declare that they have no known competing financial interests.


H Wang

F Rabhi

A Cherouat

A Gilbin

T Barriere

Experimental and numerical investigation of extrudate swell of polylactic acid via Extrusion-based Additive Manufacturing process

Keywords: EAM, Extrudate swell, Polylactic Acid, Level Set method, Experiment

Extrusion swell is a rheological phenomenon of polymers that occurs after the extrusion die exit due to relaxation of the residual molecular stress. This phenomenon is essential to ensure the accuracy and stability of components manufactured by Extrusion-Based Additive Manufacturing (EAM) process. Extrudate swell can be affected by multiple factors including material properties and processing parameters that can be coupled and it is difficult to fully understand their effects, especially in absence of accurate online measurement devices. In this study, we investigated experimentally the extrudate swell as a function of the extrusion rate, melt temperature and nozzle diameter for polylactic acid (PLA) material. The Computational Fluid Dynamics (CFD) and Level Set (LS) method in COMSOL Multiphysics were used to simulate the polymer flow out of the extrusion nozzle. The simulation results matched well with experimental results and showed that the swell rate decreases with increasing temperature and nozzle diameter and decreasing extrusion rate. These results can be used to optimize the parameters of EAM processes.

Notations a

Transition factor cp Heat capacity J K -1 kg -

Introduction

Extrusion-based Additive Manufacturing (EAM) [START_REF] Rane | A comprehensive review of extrusion-based additive manufacturing processes for rapid production of metallic and ceramic parts[END_REF], one of the most widely used Additive Manufacturing (AM) processes, involves extruding materials in their molten state through a die/nozzle and depositing them layer by layer to build complex shapes. The materials can be thermoplastics polymers [START_REF] Van Den Eynde | 3D Printing of Poly(lactic acid)[END_REF] [START_REF] Wickramasinghe | FDM-Based 3D printing of polymer and associated composite: A review on mechanical properties, defects and treatments[END_REF] or feedstocks composed of a binder and metal [START_REF] Charpentier | PIM-like EAM of steel-tool alloy via bio-based polymer[END_REF] [START_REF] Hassan | The effect of printing parameters on sintered properties of extrusion-based additively manufactured stainless steel 316L parts[END_REF] or ceramic powders [START_REF] Gonzalez-Gutierrez | Additive manufacturing of metallic and ceramic components by the material extrusion of highly-filled polymers: A review and future perspectives[END_REF], in filament or pellet shape. Compared to other AM processes, EAM processes use low-cost equipment and allows a wide range of materials to be selected.

Extrudate swell is considered as a flow instability phenomenon describing the rise of extrudate diameter compared to the die diameter and it can be also called die swell or Barus effect. When the polymer melts flow through a contraction regime, they undergo shear and extensional forces and the chains are more stretched and later gain recoiling immediately after existing the flow channel, which results in a decrease of flow velocity along the flow direction and this velocity reconstruction leads to the increase of the extrudate diameter. Thus, the extrusion swell relates to the viscoelasticity of the material. The degree of the extrusion swell is usually quantified by the ratio of the extrudate diameter versus die diameter. Extrudate swell can be affected by multiple factors, such as temperature, die geometry, extrusion stress and molecular weight of the polymer [START_REF] Sinthavathavorn | Melt rheology and die swell of PA6/LDPE blends by using lithium ionomer as a compatibilizer[END_REF]. For example, Liang et al. [START_REF] Liang | Die-swell behavior of PP/Al(OH)3/Mg(OH)2 flame retardant composite melts[END_REF] used a capillary rheometer to measure the die swell ratio of polypropylene (PP). Upon increasing the shear rate, the swell ratio increased sharply at low shear rate and then gradually reached the maximal value at the highest shear rate. When the shear rate was fixed, the swell ratio decreased with a rise of temperature. Anand et al. [START_REF] Anand | Die swell behaviour of polypropylene -An experimental investigation[END_REF] reported the extrudate swell of PP melts with capillary dies of various length to diameter L/D ratios. They found that the swell ratio decreased with increasing the L/D ratio when the shear rate is fixed. Behzadfar et al. [START_REF] Behzadfar | Extrudate swell of HDPE melts: I. Experimental[END_REF] used a capillary die with different L/D ratio and entrance angle to assess the extrusion swell ratio. They found that the swell ratio decreased with a rise of L/D ration and decreased with increasing the entrance angle until to 180°. Thus, it's essential to control the die swell especially in applications that require high precision and dimensional stability.

In order to understand the relation of viscoelasticity, processing conditions and swelling, the experimental data of die swell is essential. Early research attempt to measure the extrudate swell with pinch-off method [11] [12] [13] that the extrudate is frozen for measurement. However, this method is inaccurate as it involves the error resulted from cooling and solidifying. More accurate methods include the use of intime equipment like infrared transmission and optical cameras [START_REF] Béreaux | Optical measurement and modelling of parison sag and swell in blow moulding[END_REF] [START_REF] Tang | Isothermal flow of neat polypropylene through a slit die and its die swell: Bridging experiments and 3D numerical simulations[END_REF].

Numerical modeling is recommended to investigate the extrudate swell as it gives insights into the flow behavior inside and outside the flow channel. In the 1970s, Tanner [START_REF] Tanner | A Theory of Die-Swell[END_REF] developed a Tanner analytical method that allows the estimation of the extrudate swell ratio assuming that the die has infinite length. Tang et al. [START_REF] Tang | Three-dimensional flow simulations for polymer extrudate swell out of slit dies from low to high aspect ratios[END_REF] used finite element method including a differential multimode Phan-Thien-Tanner (PTT) constitutive model to simulate extrudate swell of polypropylene through slit dies. They found that increasing the width/height ratio from 1 to 20 contributed to 10-30% extrudate deformation in different directions and delayed the equilibrium position along the flow direction. Konaganti et al. [START_REF] Konaganti | Extrudate swell of a high-density polyethylene melt: II. Modeling using integral and differential constitutive equations[END_REF] used viscoelastic constitutive models, a 6-mode integral Kaye-Bernstein-Kearsley-Zapas (K-BKZ) model and a differential PTT model, to investigate the extrudate swell of a molten high-density polyethylene (HDPE) through capillary dies. They found that extrudate swell decreases with increasing temperature, apparent shear rate and the L/D ratio. However, the K-BKZ model overestimated and the PTT model underpredicted the experimental results, leading to an unexplained problem. Most research focused on the capillary dies with long L/D ratio compared to nozzles used in EAM processes.

In this study, an online method using a high-definition camera and a thermal camera to capture the extrusion swell and temperature distribution downstream of the nozzle exit was adopted. Processing parameters such as temperature, extrusion rate and nozzle diameter of an EAM printer were studied. The experimental results were fed into the Two-Phase Flow module of the COMSOL Multiphysics software with the Level-Set (LS) method to model the extrusion swell phenomenon. These results and analyses were used to predict and optimize the extrudate swell in EAM process.

Experimental materials and methods

Material used

Polylactic acid (PLA) is the most widely used polymer in the EAM applications [START_REF] Tümer | Extrusion-Based 3D Printing Applications of PLA Composites: A Review[END_REF] due to its biodegradability and environmental properties. The PLA (PLI005, Natureplast, France) was used in this study to characterize extrudate swell during EAM process. Differential Scanning Calorimetry (DSC) was used to analyze the thermal properties of PLA. Their physical and thermal properties are shown in Table 1.

Measurements of rheological behaviors and modeling

The shear viscosity of PLA in the molten state was characterized to determine the rheological behaviour using a plane-cone rotational rheometer (HAAKE MARS 60 Rheometer) for low shear rates and a capillary rheometer (Rosand RH2000) for high shear rates. The experimentally measured viscosity as a function of shear rate for temperature range of 180-220℃ is shown in Fig. 1. The shear viscosity decreases with increasing the temperature and shear rate, illustrating that the PLA has a thermally dependent pseudoplastic type of rheological behaviour in our case.

( )(1 ( ) ) n a a      -  = + - + (2) 0 11 exp ( ) a T E a R T T  =-  (1) 1 0
where 0

 is the viscosity for zero shear rate,   the viscosity at infinite shear rate, λ the characteristic time of the material, a the parameter adjusting the transition regime between the Newtonian plateau and the shear-thinning regime, n the pseudo-plasticity index,  the shear rate, T the temperature, 0 T the reference temperature (here was 190 C  ), a E the activation energy and R the gas constant.

The values of the identified parameters were obtained by the least square method using the experimental data and are summarized in Table 2. The value of the pseudoplastic index n (<1) indicates that the PLA is pseudoplastic.

Experimental characterization of die swell and temperature distribution

A dynamic camera (Basler Ace acA3088) with a macro lens was used to measure the extrudate swell and an infrared camera (Optris PI 450i) was selected to quantify the temperature of the extrudate. Both cameras were installed perpendicular to the extrudate flow direction as shown in Fig. 2.

The EAM (printing) equipment is a screw-based 3D printer in which the material is fed and extruded by the rotation of a single screw from thermoplastic polymer pellets. For each photograph, the outline of the extrudate was detected using ImageJ after calibration by an object of known size placed in the same plane where the extrudate flowed perpendicular to the camera.

Fig. 2 The setup of the dynamic and the infrared cameras.

Before the measurements of the temperature, the emissivity of PLA was determined by adjusting the emissivity of the thermal camera to achieve the same temperature captured by the thermal couple which was placed on the surface of the melted materials. A small control region of 3×3 pixels was selected near the die exit to measure the extrudate temperature and three other areas were selected to evaluate the evolution of temperature along with the extrudate, as shown in Fig. 3(a). When PLA flows out of the nozzle exit, it undergoes a creep phenomenon that it deforms (expansion) to relax the build-up residual stresses generated by the shear and compression forces to the polymer chain network during extrusion through the extruder. With contact to the nozzle surface, the expanded material forms a bulb at the onset of the extrusion. As the extrudate flows, the bulb leaves the nozzle surface and the hanging mass of the extrudate results in an apparent necking. Both the bulb and necking part are not suitable for charactering the extrudate swell. The extrudate swell was measured at the top of the extrudate where the extrudate swelled to the maximal diameter.

As the extrudate flowed downward past the observation points (P1-P4), the temperature at each point as a function of time was recorded in real-time, as shown in Fig. 3(b). The curves can be divided into three domains. In the extrudate arrival domain, the temperature remains the ambient temperature until the arrival of the extrudate, except for P1, which starts at a higher temperature because it is right under the nozzle and has already some residual PLA. The unsteady state represents the variation of temperature resulting from the extrusion process that the newly fed cooler material exchanges the heat or mix with the material that in front of it. As the material of higher temperature heated in the static state is extruded out, the temperature evolution becomes stable and then the measurement of extrudate swell is started in this steady state. The measurement of extrudate swell started after the stabilization of the temperature profile. In order to maintain the consistency of the die swell, before each measurement when the extrudate reached its thermal steady state, the extrudate under the nozzle was clipped without stopping the extrusion process. The maximal diameter of the extrudate within 2mm below the nozzle exit when the extrudate reached a total length of 6mm was measured for determining the die swell. 3 Finite element simulation

Numerical models

The level-set method LS is a numerical technique for tracking interfaces and shapes using a fixed mesh. The level-set method LS was developed in 1979 by Alain Dervieux [START_REF] Dervieux | A finite element method for the simulation of a Rayleigh-Taylor instability[END_REF], and subsequently popularized by Stanley Osher and James Sethian [21] [22]. It has become popular in many disciplines, such as image processing, optimization and computational fluid dynamics. The interface is represented by an isocontour of a globally defined function, the level set function  . Let a domain

( ) ( ) ( ) ,0 ,0 ,0 t x x t x x t x x    + -      =         (4) If the curve ( , ) tx  moves in the normal direction () nx    =
 with a speed u, then the level-set function satisfies the level-set equation to solve advection equation to propagate 0  = as:

( , ) ( , ) 0 tx tx t   +  =  u ( 5
)
Given at time t = 0, 0 (0, ) x  = for initial domain (0)  is given, and its motion over time, it is possible to know ( , ) tx  at any time using a second-order in time Runge Kutta scheme as:

( ) t t t t t t t given     + = -       u ( 6
)
The numerical simulations were carried out with the Two Phase Flow (TPF) solver of Level Set (LS) method in COMSOL Multiphysics based on finite element method. D. Kazmer et al. [START_REF] Kazmer | Concurrent characterization of compressibility and viscosity in extrusion-based additive manufacturing of acrylonitrile butadiene styrene with fault diagnoses[END_REF] characterized the compressibility of ABS in the material extrusion process with an instrumented hot end design and stated that the compressibility effects can be significant on the transient variances in material extrusion. In this study, the compressibility was not considered and we assumed that the polymer flow is incompressible and the external forces are negligible compared to the internal forces, the Navier-Stocker equations ( 7) and ( 8) governing the viscous flow and the LS transport equation [START_REF] Anand | Die swell behaviour of polypropylene -An experimental investigation[END_REF] tracking the interface between two immiscible fluids are expressed as follows:

( ( )

T p + + t      +      +    u u u) = -I g u u (7) 0  = u (8) ( ) ( )
1 ls t           +   =    - -     u ( 9
)
where  is the fluid density, u the velocity vector, p the pressure,  the dynamic viscosity, g the acceleration of gravity, I the identity vector,  the volume fraction of the fluids varying between 0 and 1 across the free surface and is constant at 0 or 1 in the bulk of the two fluids (the value of 0.5 is assigned to the interface between the polymer and the air), t the time,  the reset parameter and ls  the interface thickness control. In order to model the heat transfer in the flow simulation, the advection-diffusion equation ( 10) is added by using a partial differential equation (PDE) in COMSOL:

( ) 2 p Tk TT t ρc   +  =      u ( 10
)
where T is the temperature, k the thermal conductivity, ( ) ( )

1 air polymer      = + - (11)
where air  and polymer  represent the diffusivity of air and polymer respectively.

The simulation for predicting the extrudate swell at the nozzle exit is realized by coupling different modules in COMSOL Multiphysics software, as illustrated in Fig. 4.

Geometrical modeling, boundary conditions and computational grid

The extrusion simulation of PLA was modeled as 2D axisymmetric and the geometry is shown in Fig. 5. The yellow part is the nozzle channel filled initially by melted polymer and the green part is the air domain where enters then the polymer. The nozzle radius is R=0.38 mm or 0.78 mm and the nozzle length is L=0.6 mm. The boundary conditions are shown in Fig. 5, including no slip condition and constant extrusion temperature m T on inner wall of the nozzle, ambient pressure on the outlet and environment temperature env T (25 C  ) in the air cavity. In fact, the volumetric flow rate is not a constant value for a given rotation speed of the screw, but also denpends on the temperature and nozzle geometry. Besides, the realistic extrusion temperature varies with the rotation speed. Thus, the input boundary conditions of inlet velocity V and extrusion temperature m T were assigned corresponding to each experimental conditions. The 2D computational grid is physics-controlled meshed (67241 elements) with free triangles and refined in the nozzle contraction region, as shown in Fig. 5. Fine mesh size is chosen at the interface polymer-air (55461 elements in the area 2 and 10138 elements in the area 3) and at a very refined mesh (1642 elements) is performed at the exit of the nozzle extruder (area 1).

Results

Effects of extrusion rate, nozzle diameter and extrusion temperature on the exit temperature and volume flow rate

The printing process used (extrusion temperature, nozzle diameter and extrusion rate) to study the die swelling phenomenon are defined in Table 3. Fig. 6(a) shows the evolution of volume flow rate as a function of the extrusion rate. It also showed the effect of nozzle diameter and extrusion temperature. The slope can be considered as a factor indicating the feeding facility. The feeding process can be considered as a confrontation process of the shear force driven by screw against the force due to the pressure drop inside the nozzle. Nozzle with larger diameter at higher extrusion temperature possesses a high slope. Flowing through a smaller nozzle diameter and under lower temperature, the melt undergoes higher shear force, which hinders the feeding of the material. For instance, the slope of PLA with nozzle diameter =1.56mm at 210 ℃ is 0.85, higher than that (0.28) of nozzle diameter =0.76mm at 190 ℃. Hyvärinen et al. [START_REF] Hyvärinen | The modelling of extrusion processes for polymers-A review[END_REF] reported that the extrudate flow rate will tend to decrease with screw extrusion rate due to the increased slip at the screw, which was not detected in our study. Fig. 6(b) shows the evolution of exit temperature as a function of the extrusion rate for different nozzles and extrusion temperatures. The exit temperature tends to decrease as the extrusion rate increases, which can be explained by the reduction of the residence time of the material to heat inside the extrude. A larger nozzle diameter facilitates the material flow, which also reduces the residence time inside the extruder. Reasonably, a higher extrusion temperature results in a higher exit temperature. In contrast to the assumption of the filament-used EAM processes that the volume flow rate of the melt is proportional to the feeding rate in conditions that the friction force of the extruder gear wheels do not exceed the strength of the filament material and the melt is incompressible, the feeding force in pellet-used EAM is driven by the friction and compression of the material, which results the feeding rate of the pellets a nonlinear parameter in relation to the processing parameters. Thus, to reflect the actual feed rate, the extrusion rate is represented by the corresponding volume flow rate. The evolution of exit temperature in Fig. 6(b) is replotted as a function of volume flow rate in Fig. 7.The measured nozzle exit temperatures for PLA were lower than the assigned temperature, which is consistent with what has been reported in [25][26]. This is due to the fact that the heat transfer between the liquefier walls and the material is insufficient for the material to reach its assigned temperature. At the same time, as the flow rate increases, the difference between the extrusion temperature and the outlet temperature increases. This is because a higher extrusion rate results in a shorter residence time for the material to pass through the liquefier, making the temperature difference more pronounced. Different nozzles have different volume flow rate under the same processing parameters, but it is interesting to note that their exit temperatures are close in the volume flow rate overlap region, indicating that the reduced exit temperature is primarily due to the volume flow.

Fig. 7 The effect of flow rate and extrusion temperature on the exit temperature According to the above-mentioned analysis, it turns out to be necessary to account for the real exit temperature and volume flow rate at their stable state for the input of finite element simulation, though the presentation of the results utilize still the assigned values.

Comparison of the experimental and numerical swell

In our work, the experimental results obtained from the die swelling with two nozzle geometries of different diameter were supported with Tanner's model. He modeled the die swelling (B) as a function of the shear stress 𝜏 𝑝 at the wall as described by equation [START_REF] Dutta | A study of parison development in extrusion blow molding[END_REF], which can also be expressed as a function of the shear rate 𝛾 ̇as equation (13):

2 (1 ) m p B A C = + + (12) 2 (1 ) m B A C = + + (13)
The parameters A, C and m are identified by the least squares method and shown in Table 4. The experimental and Tanner model results of die swell ratio versus shear rate for both nozzles are shown in Fig. 8. The evolution of swelling as a function of shear rates obtained from the experimental point of view correlates well with the results obtained by the Tanner model. For the largest die diameter, the die swell ratios are close to the Newtonian limit of the die swell rate of 1.1, indicating that under the processing conditions, PLA exhibits more non-Newtonian viscous behavior than elastic behavior. For a nozzle diameter of 0.76mm, the die swell ratio is well above 1.1 indicating viscoelastic behavior. A comparison of experimental and numerical results for PLA at extrusion temperature T=200 ℃,extrusion rate of V=20 rpm for nozzle diameter of =0.76 mm is shown in Fig. 9. Fig. 9(a) shows the temperature distribution captured by the thermal camera as well as the temperature field of simulation (in red dotted box), both revealed the cooling behavior of the extrudate along the flow direction while the simulation results gave an insight to the cooling behavior from the center to the edge. However, the temperature on the surface of the extrudate in simulation seems to reduce faster than the experimental results. The temperature contour in both cases cannot reflect the true contour of the extrudate swell. Fig. 9(b) shows the extrudate shape captured by the dynamic camera as well as the field of volume fraction  of PLA (  =1, red) in simulation (in red dotted box). After a sharp swelling, the diameter of the extrudate increases gradually to a stable regime where indicates the arrow. The summary of the results of experiment and simulation is listed in Table 5. The extrudate swell is usually characterized by swell ratio B: Fig. 10 presents the transient behavior of the extrudate at 200 ℃ with nozzle diameter D=0.76mm. The color contour represents here the volume fraction of level set equation, the value of 0.5 is regarded as the interface of PLA and air. The extrudate creeps a bit during the initial flow but didn't form a bulb as obvious as that described in Fig. 3(a) in lack of a precise initialization. Logically, higher volume flow rate leads to longer lengh and more significant die swell at the same time due to the higher shear rate experienced in the nozzle. As the extrudate continues flowing, it undergoes a fast and significant swelling right out of the nozzle exit and tends to equilibrate between the creep behavior and the tension forces from the hanging mass, after that, the diameter of the extrudate decreases in the flow direction due to the gravity effect. Fig. 12 shows that the swell ratio decreases with increasing the temperature. This can be explained by the effect of temperature on the polymer chain mobility that the higher temperature promotes the movement of polymer chains and contributes to their reorientation and release of deformation resulted from the shear and compression forces in extrusion. The deviation between the simulation results and experimental results is within 4%, which proves that the simulation correlates well with the experimental results and can be used to investigate and optimize the processing parameters. According to the above results, it can be concluded that reducing the extrusion swell can be achieved by decreasing the extrusion rate, increasing the extrusion temperature and nozzle diameter, but the printing speed and dimensional resolution should be considered.

Fig. 12 Effect of temperature on the swell ratio

Conclusion

The extrudate swell ratio of PLA in EAM process flowing through a nozzle was experimentally characterized. The processing parameters of temperature and extrusion rate, and the nozzle diameter were studied. The relations among the screw rotation speed, extrusion temperature and volume flow rate were determined, and the real values were taken as the input of numerical simulation. The numerical simulation of extrudate swell of PLA was realized using the TPF module and LS method in COMSOL Multiphysics. The results of experiment and simulation agreed well showing a maximal error of 4%. Both results showed that the extrusion swell decreases with the increase of the extrusion temperature and the nozzle dimeter, and with the decrease of extrusion rate.

Fig. 1

 1 Fig.1 Shear viscosity vs. shear rate curve of PLA. The Carreau-Yasuda model which accounts for both the Newtonian regime (constant viscosity at low shear rate) and shear-thinning regime (viscosity decreases with increasing the shear rate) was used as the rheological model and Arrhenius model was used to account for temperature. The shear viscosity ,) T  ( and the two models are expressed in equations: 0 , ) ( , )

Fig. 3

 3 Fig.3 Temperature measurement (a)Position of the measure points along the extrudate (b) Evolution of temperature for PLA at 200 ℃ and V=20 rpm (zoom sur la figure (a) et enlever la partie noir).(figure b P1 P2) Temperatures of 190℃, 200℃ and 210℃ were tested, and each temperature was related to four extrusion rates (Rotation speed of the screw at 10, 20, 30 and 40 rpm). Two nozzles with different diameters (D=0.76mm and D=1.56mm) were studied for comparison. Unlike the filament-used EAM process where the extrusion rate is explicit, the extrusion rate of the pellet-used EAM process using a single screw to extruder the pellets represents only the rotation speed of the screw. So, during each measurement, the mass of the extrudate was recorded for converting to material's volume flow rate.

.

  This function is assumed to take positive values inside the region delimited by the curve  and negative values outside. In order to investigate the flow in the domain  , let the level set function divide the domain in two parts +- =   , where:

  the thermal diffusivity  within the polymer/air interface and is expressed (same principle for density and viscosity) in (11):

Fig. 4

 4 Fig.4 Schema of simulation implementation for extrudate swell.

Fig. 5

 5 Fig.5 Geometry of simulation model, boundary conditions and computational grid.

Fig. 6

 6 Fig.6(b) Effect of nozzle diameter, extrusion rate and extrusion temperature on the volume flow rate on the exit temperature

Fig. 8

 8 Fig.8 Experimental results and Tanner's model of die swell ration as a function of shear rate for nozzle diameter of (a) D=0.76mm and (b) D=1.56mm.

  extrudate diameter and D  the nozzle diameter.

Fig. 9

 9 Fig.9 Experimental and numerical results for PLA at 200 ℃, 20 rpm (a) Temperature field captured by thermal camera (simulation result in dotted box); (b) Extrudate contour captured by dynamic camera (simulation result in dotted box).

Fig. 10

 10 Fig.10 Transient behavior of the extrudate at 200℃ for nozzle diameter of 0.76mm at volume flow rate of (a) 4.07 mm 3 /s; (b) 8.15 mm 3 /s;(c) 12.22 mm 3 /s;(d) 14.74 mm 3 /s.

Fig. 11

 11 Fig.11compares the swell ratio between experimental and numerical results as a function of the flow rate at 200℃, which shows that the swell ratio increases with increasing the flow rate for both nozzles. The nozzle of larger diameter has lower swell ratio, which is in accordance with literature[START_REF] Serdeczny | Experimental and analytical study of the polymer melt flow through the hotend in material extrusion additive manufacturing[END_REF].

Fig. 11

 11 Fig.11 Comparison of experimental and simulation results of die swell for two nozzles at different flow rate and at extrusion temperature of 200℃

Table 1

 1 Physical and thermal properties of PLA.

	Polymer	PLA
	Density 𝝆 (g.cm -3 )	1.25
	Heat capacity 𝒄 𝒑 (J.kg -1 .K -1 )	2100
	Thermal conductivity 𝒌 (W.m -1 .K -1 )	0.205
	Glass transition temperature 𝑻 𝒈 (℃)	52 -65
	Melting temperature 𝑻 𝒇 (℃)	169 -179

Table 2

 2 Summary of the identified parameters of Carreau-Yasuda law and Arrhenius law.

	𝑻 (℃)	𝜼 𝟎 (𝐏𝐚 𝐬) 𝜼 ∞ (𝐏𝐚 𝐬)	𝝀 (𝐬)	𝒏	𝒂	𝑬 𝒂 (kJ mol -1 )
	180 190	1187 917	0 0	0.04 0.10	0.28 0.34	0.96 1.35	105.83

Table 3

 3 Effects of processing parameters on the extrudate swell phenomenon.

Table 4

 4 Parameters of die swell model of Tanner as a function of shear rate.

	Nozzle	A	B	m	R
	diameter
	D1=0.76mm	0.1000	0.0059	0.0205	0.9879
	D2=1.56mm	-0.0141	0.0687	0.0201	0.9990

Pas coherent avec figures 11 et 12 poru swell ratio B oiu swell ratio

Table 5

 5 Comparison of extrudate swell ratio between experimental and numerical results for PLA at 200℃.

		Volume flow rate	4.07	8.15	12.22	14.74
		(mm 3 /s)
	D=0.76mm	𝐵 of experiment	1.18	1.22	1.23	1.24
		𝐵 of simulation	1.19	1.24	1.27	1.27
		Error (%)	0.85	1.64	3.25	2.42
		Volume flow rate	8.34	16.29	24.44	32.58
		(mm 3 /s)
	D=1.56mm	𝐵 of experiment	1.05	1.07	1.08	1.09
		𝐵 of simulation	1.06	1.09	1.11	1.12
		Error (%)	0.95	1.81	2.64	2.68

Statements and Declarations

Funding and acknowledgements

The authors gratefully acknowledged the financial support received from the China Scholarship Council (CSC). Experimental tests were carried out thanks to SupMicroTech-ENSMM workshop and local plateforms (AMETISTE for metrology and tribology, MIFHYSTO for cutting forces). This work was carried out within the Manufacturing 21 network, which gathers about 20 French research laboratories. The covered topics are the study and modelling of the manufacturing processes, especially the numerical simulation and investigation of extrudate swell of polylactic acid (PLA) via extrusion-based Additive Manufacturing process.

Competing Interest

The authors have no competing interests to declare that are relevant to the content of this article. This work was realized on a 3D printing equipment at the department of applied mechanics of Besançon, the presented results are free of commercial rights.

Author contributions

H. Wang: Design; Methodology; Data analysis; Visualization; Writing-Original draft.

F. Rabhi: Data analysis; Writing-Editing A. Cherouat: Supervision; Review -Editing A. Gilbin: Methodology T. Barriere: Supervision; Review-Editing