Damien Lemaitre
email: damien.lemaitre@g2elab.grenoble-inp.fr

Benoit Sarrazin
email: benoit.sarrazin@g2elab.grenoble-inp.fr

Yohan Wanderoild
email: yohan.wanderoild@edf.fr

Higher Winding Efficiency in the WPT System Thanks to Homogeneous Current Distribution in the Litz Wire

The design of the Litz wire has a strong influence on the performance on wireless power transfer systems. The number and arrangement of the bundles are important factors, whose influence on the distribution has so far been complex to study. This paper presents a new method to easily analyse the current distribution in the Litz wire as a function of these factors. The effect of frequency on the observed phenomena will be explored in simulations and validated experimentally.

INTRODUCTION

The development of electric vehicles (EVs) is leading to a growing need for chargers. EVs are largely charged with plug-in chargers. However, a large amount of research is being done on wireless chargers, which are reported to be over 95% efficient [START_REF] Onar | Modeling, Simulation, and Experimental Verification of a 20-kW Series-Series Wireless Power Transfer System for a Toyota RAV4 Electric Vehicle[END_REF]- [START_REF] Sarrazin | Bidirectional Wireless Power Transfer System with Wireless Control for Electrical Vehicle[END_REF]. One of the interests of inductive charging is the vehicle-to-grid operation (V2G), allowing a bi-directional energy transfer to charge and discharge the EV battery, or to support the grid with energy stored in the battery. The development of intermittent renewable energy sources requires storage resources, which can be provided by EVs equipped with V2G chargers. Maximising the efficiency of these systems becomes more important than that of unidirectional chargers because of the higher amount of energy transferred. A V2G wireless power transfer system (WPTS) for EV charging consists of converters, resonant capacitors and a coupler as shown in Figure 1. The coupler itself is very often composed of windings, magnetic materials and component to ensure the magnetic shielding of the system. Conductors constituting the windings are the most expensive and loss-making of all these components for systems operating in the 20 kHz-200 kHz frequency range [START_REF] Shekhar | Economic Viability Study of an On-Road Wireless Charging System with a Generic Driving Range Estimation Method[END_REF]. These windings are often made of multi-stranded conductors called Litz wire to overcome the skin effect. Figure 2 shows a Litz wire with a 5x5x16 structure, composed of 5 bundles, each composed of 5 sub-bundles, each composed of 16 elementary strands.

The present paper proposes a new method to study the impact of the Litz wire design at the bundle level on the current distribution.

Part 1 establishes the Litz wire model showing the role of the inter-bundle mutual inductances on the impedance of the bundles. Part 2 introduces the proposed method using simple electromagnetics (EM) simulation with finite element method and electrical circuit simulations. The mutual inductances calculated with the EM simulation are used in the simulation of the Litz wire electrical circuit to observe its current distribution. The 2D simulations will allow to understand the effect of twisting on the current distribution without having to perform 3D simulations. Part 3 presents the results obtained on the distribution for different study cases. The influence of the twist, the number of bundles, the presence of bundle in the centre of the Litz wire, and the working frequency on the current distribution is presented. Part 4 presents the experimental results on two 20 kW WPTS using different Litz wires. The influence of the frequency on the phase shift of the current between bundles is measured.

LITZ WIRE MODEL

Litz wire is an assembly of conductors insulated from each other and arranged electrically in parallel. The electrical circuit developed here is on the bundle scale. The approach is possible at the scale of a sub-bundles or strands.

A. Notations

Let us consider the inductance and the resistance of bundle i of an n-bundle Litz wire, noted 𝐿 and 𝑅 respectively, where n is an integer greater than or equal to 2, i is an integer from 1 to n. The impedance 𝑍 of bundle i is therefore described by equation ( 1) if no coupling is considered.

𝑍 _ = 𝑅 + 𝑗. 𝜔. 𝐿 (1)
In the case of electric vehicle chargers, the reactive part is much higher than the resistive part due to the relatively high inductance and frequency. At low frequencies, the impedance is resistive, with a resistance close to the DC resistance. With these considerations and to simplify the analysis, the bundle resistance considered is the DC resistance, noted 𝑅 . The bundles also have the same selfinductance noted 𝐿 due to their identical geometry. Equation (2) below repeats equation ( 1) with these notations.

Z _ = R + j. 𝜔. L (2)
Let us now consider the inter-bundle couplings. The mutual inductance and the coupling between bundles i and j are respectively noted 𝑀 and 𝑘 , i and j being two different integers between 1 and n. The expression for the mutual inductance 𝑀 is given by equation (3). The mutual flux shared by the bundle i and j noted 𝜙 is defined by the equation (4). It is generally agreed that 𝑀 = 𝑀 . In order to simplify the model, a three-bundle Litz wire is considered, i.e. n = 3. The notations T1, T2 and T3 refer to bundles n°1, n°2, and n°3 respectively.

𝑀 = 𝑘 . 𝐿 . 𝐿 = 𝑘 . 𝐿 (3)
𝜙 = 𝑀 . 𝐼 (4)

B. Bundle impedance

First, the total impedance of the T2 bundle is calculated. The impedances induced by bundles T1 and T3 on bundle T2, are respectively noted 𝑍 and 𝑍 . The expression of 𝑍 is given by the Lenz-Faraday law in receiver convention. The mutual flux 𝜙 emitted by the current 𝐼 from bundle T1 and seen by bundle T2 induces a voltage 𝑉 in bundle T2, whose the complex expression is given by the equation [START_REF] Lemaitre | Efficiency-Oriented Design of Litz Wire for Several kW Power Experimented on 20 kW Prototype[END_REF]. The expression for the voltage induced 𝑉 by the current 𝐼 is given by equation [START_REF] Sullivan | Simplified design method for litz wire[END_REF]. The electrical diagram of this three-bundle Litz wire is shown in Figure 3. The influence of the mutual inductances on the bundle impedances and thus on the current distribution in the Litz wire is established. We saw that different coupling coefficients and thus different mutual inductances cause different impedances in the bundles. As the voltage is the same across the bundles, current heterogeneities appear.

The following section presents a method to study the relationship between the current distribution and the number of bundles, their assembly and frequency.

PROPOSED METHOD

We will first look at the steps in the study of the current distribution in the Litz wire. The geometrical considerations will then be presented.

Principle of the method

The EM simulation of the Litz wire allows the calculation of the mutual inductances between the bundles. The number of simulations required is equal to the number of bundles. It is indeed necessary to impose 1 A in each bundle, one after another, in order to be able to deduce from the magnetic fluxes obtained the mutual inductances between the fed bundle and the others. An example of the

Paper n°1346

CIRED 2022 Workshop 3/5

EM and electrical circuit simulations of the Litz wire is given in Figure 4 (a) and (b) in the case of three bundles.

Once the mutual inductances between the bundles have been obtained, a simulation of the Litz wire electrical circuit is required to obtain the current distribution.

Considering a Litz wire with n bundles, the electrical circuit is composed of n coupled inductances, each in series with a resistor. A current source is used to simulate the total current in the Litz wire at the working frequency.

The self-inductances and mutual inductances obtained from the previous simulations define the introduced coupled inductances.

(a) (b) Figure 4 : EM (a) and electrical circuit (b) simulations of a Litz Wire composed of 3 bundles.

Geometric considerations

The resulting current distribution depends on the couplings between the bundles, and therefore their representation in the EM simulations. Three cases can be distinguished and are shown in Figure 5 using the example of a 3-bundle Litz wire. The case Figure 5 (a) represents an ideal Litz, the bundles are equidistant and undeformed. The coupling coefficients are therefore equal in this case: 𝑘 = 𝑘 = 𝑘 . Figure 5 (b) introduces the case of a different distance between the bundles, which results in different coupling coefficients. If the deformation of the bundle cross-section during twisting is considered, case Figure 5 (c) could be one of the possible results obtained after manufacturing. The bundles are however compacted during the twisting phase, so studying the influence of the distance between bundles would not be relevant. The experimental results obtained in part 4 show that the current distribution in the Litz wire respecting the state-of-the-art design recommendations remains rather homogeneous despite the deformation of the bundles. It is also difficult to define which deformations the bundles have undergone. The simulations carried out in part 3 thus follow Figure 5 (a) type constructions.

(𝑎) (𝑏) (𝑐) Figure 5 : Arrangement of the three bundles constituting the Litz wire in the ideal case (a), the case of nonequidistance between the bundles (b) and the case of bundle deformation (c).

Twisting considerations

Only a 3D simulation can represent a twisted Litz wire, which is however time consuming and would increase the complexity of the proposed method. As the objective of twisting is to allow each bundle to be subjected to the same magnetic field, a 2D plane simulation allows the same conditions to be applied. This assumes that the twisted pitch has been carefully chosen. Under this assumption, it is possible to compare a perfectly twisted Litz wire simulated in 2D plane with a non-twisted Litz wire represented in an axisymmetric 2D simulation.

SIMULATION RESULTS

The influence of twist, number of bundles and frequency on the current distribution of the Litz wire is studied in this section. These studies are carried out with the objective of providing a better understanding of which wire to select for a WPTS application for EV charging. Therefore, the working frequency of 85 kHz recommended by the SAE standard is chosen, except for part 3.3 on frequency dependence. Litz wires consisting of 2 to 7 bundles are studied.

Twisting influence on current distribution

In order to evaluate the influence of the twist on the current distribution in the Litz wire, two simulations are performed. The first simulation in 2D plane of straight bundles is based on the hypothesis of a Litz wire with a sufficiently small twist pitch to neglect the influence of the external magnetic field on the current distribution. This assumption is valid in the windings on condition that the twist pitch is suitable for the application. The second 2D axisymmetric simulation allows to simulate a coil made of a non-twisted Litz wire. Figure 7 (a) and (b) show the current distribution of a fivebundle Litz wire. Its representation in the 2D axisymmetric simulation is shown in Figure 6. The same representation is used for the 2D plane simulation. The values are available in Table 1.

(𝑎) (𝑏) Figure 7 : Current distribution in the Litz wire simulating bundle twist (a) and without twist (b) Similar results are obtained for different numbers of bundles as long as they are arranged at the periphery as shown in Figure 6 and not in the centre of the Litz wire. Increasing the number of bundles often means that the bundles end up in the centre of the Litz wire, which causes amplitude and phase heterogeneities between the bundle currents, as shown in the following section.

High Numbers of Bundle

This study is carried out on two 7-bundle Litz wires, arranged as shown in Figure 8 (a) and (b). The second structure uses a non-conductive material in the centre so that the bundles can be arranged on the outside as proposed by some manufacturers. The simulations are performed in 2D axisymmetric.

(𝑎) (𝑏) Figure 8 : Litz wire consisting of 7 bundles with (a) and without (b) a bundle in the centre of the wire.

The current distribution obtained for the Litz wire illustrated in Figure 8 (a) shows a phase shift between the current in the central bundle and those in the peripheral bundles of 180°. The waveforms are available in Figure 9 (a). This phase shift is not apparent in the absence of a central bundle as shown in the waveforms in Figure 9 (b). The total current in the Litz wire is imposed at 27 A in both cases. However, the absolute sum of the bundle currents is 31.8 A , compared to 27 A in case Figure 9 (b). An increase in the bundle current is indeed necessary to guarantee the total current set at 27 A . Case Figure 9 (a) thus presents 39 % more losses than case Figure 9 (b), considering the difference of the total current only, at equivalent resistance.

(𝑎) (𝑏) Figure 9 : Current distribution of a 7-bundles Litz wire with a bundle in the centre (a) and without bundle in the centre (b).

The effect of these amplitude and phase heterogeneities are related to the different induced impedances in the bundles of the Litz wire as shown in the model presented in part 2. The frequency is therefore a key parameter, and the following section is dedicated to its influence on the current distribution.

Frequency dependence

The 7-bundle Litz wire with a central bundle is used in this study. The current distribution is studied for a working frequency varying from 1 Hz to 1 MHz. The natural and induced reactances of the bundles are proportional to the frequency, which is why the resistive part noted 𝑅 is dominant at low frequency. For this reason, the current distribution is perfectly homogeneous at low frequencies. The phase shift between the current of the central bundle and those at the periphery is represent as a function of frequency in Figure 10. The phase shift quickly becomes very pronounced, reaching 76° at 100 Hz and 164° at 1 kHz. The distribution shown in Figure 9 (a) is therefore very close to the maximum phase disparity, being at a frequency of 85 kHz. In order to compare these theoretical results with reality, the current distribution of two different Litz wires was measured. The results are analysed in the following section.

EXPERIMENTAL VALIDATION

Setup

The tests in this section were carried out on two 20 kW WPTSs, designed for a resonant frequency of 81 kHz. These WPTSs differ only in the Litz wire used to make up their coupler windings. The characteristics of Litz wires are presented in previous work [START_REF] Lemaitre | Efficiency-Oriented Design of Litz Wire for Several kW Power Experimented on 20 kW Prototype[END_REF].

Paper n°1346

CIRED 2022 Workshop 5/5

Results on current distribution

The current distribution was measured in paper [START_REF] Lemaitre | Efficiency-Oriented Design of Litz Wire for Several kW Power Experimented on 20 kW Prototype[END_REF] for Litz wires n°1 and n°2 of 16 bundles and 5 bundles respectively. The current distributions obtained in measurements and in simulation for the Litz wire n°2 are presented in Table 2. The measurements and simulations of the currents show a homogeneity of phase for the Litz wire n°2. Even if the heterogeneity of the amplitudes is similar between the measurements and the simulation with an average relative error of 4.1 %, the causes are different. These causes are the non-consideration of twisting in simulation, while the mechanical deformation of the bundles seems to be the reason for this heterogeneity of the measured current amplitudes. [START_REF] Lemaitre | Efficiency-Oriented Design of Litz Wire for Several kW Power Experimented on 20 kW Prototype[END_REF]. These phase shifts reach near phase opposition with values of 162°, 178° and 165°. The bundles 9, 10 and 14 seem to be effectively located in the centre of the Litz wire. This confirms the results obtained in simulation on the phase shift obtained between the current of the central and peripheral bundles.

Results on the frequency influence

Further measurements were made to observe the frequency dependence of the phase shift of the bundle currents. It should be noted that the power supply we used did not allow for more than 2,600 Hz. The measurements of the phase shift between the bundle current n°14 and the Litz wire n°1 current are available in Table 3. An evolution similar to the results in Figure 10 can be observed.

LITZ WIRE DESIGN

The recommended structures [START_REF] Sullivan | Simplified design method for litz wire[END_REF] have at the first stage of twisting a number of elementary strands fixed by the frequency and the diameter of a strand. The number of bundles at the higher stages of twisting should be less than or equal to 5, in particular to avoid the phenomena seen previously. For instance, at the 85 kHz frequency recommended by the SAE standard, wires with a 3x4x20, 5x5x20, or 4x5x5x20 design would be suitable, depending on the copper cross-section of interest. As each stage of twisting involves an additional cost, a structure with a nonconductive material in the centre allows the number of bundles to be increased at the last stage of twisting without the appearance of a phase shift between the currents. This design can be an intermediate solution between a design with n stages of twisting lacking copper cross-section and a more expensive wire with n+1 stages of twisting. The possibility and economic interests of making such a structure remains to be studied with the manufacturer as it is not very common until now. The additional crosssection and volume induced by the non-conductive material is also to be considered.

CONCLUSION

A new method that is easy to implement dedicated to the analysis of the current distribution in the Litz wire was presented. The amplitude and phase heterogeneities of the currents in the bundles were simulated and experimentally validated for several designs. Based on these results, an additional design to the existing recommendations is proposed.

Figure 1 :

 1 Figure 1: Electrical schematic of a WPTS between the grid and the vehicle battery with a Series-Series resonant topology.

Figure 2 :

 2 Figure 2 : Example of a LW design that consists of 5 bundles, each consisting of 5 sub-bundles, each consisting of 16 strands: 5x5x16.

  𝑉 = 𝑗. 𝜔. 𝜙 = 𝑗. 𝜔. 𝑀 . 𝐼 (5) 𝑉 = 𝑗. 𝜔. 𝜙 = 𝑗. 𝜔. 𝑀 . 𝐼(6)By dividing the induced voltages 𝑉 et 𝑉 by the current 𝐼 , the expressions for the induced impedances 𝑍 and 𝑍 appear: of the T2 bundle given by equation (9) is then obtained by summing 𝑍 and 𝑍 to the expression (2). Impedances 𝑍 and 𝑍 of bundles T1 and T3 are obtained in similar ways. 𝑍 = 𝑅 + 𝑗. 𝜔. 𝐿 + 𝑗. 𝜔. 𝑀 .

  Figure 3 (a) represents the electrical circuit of the Litz wire at the bundle scale, considering their couplings. The latter are translated using the induced inductances in Figure 3 (b).

Figure 3 :

 3 Figure 3 : Electrical circuit of a three-bundle Litz wire at the scale of the bundle showing couplings (a) and induced impedances (b).

Figure 6 :

 6 Figure 6 : 2D simulation of a 5-bundle Litz wireIn order to be able to compare with the experimental results, the total current in the Litz wire imposed in the simulation of the electric circuit is identical, i.e. equal to 27 A root mean square (rms), i.e. 38.2 A in amplitude. A perfect current distribution results in a current of 7.6 A in each bundle. This is shown in the waveforms in Figure7(a) for a properly twisted Litz wire. A non-twisted Litz wire constituting a coil shows instead a heterogeneity of current in amplitude as shown in Figure7 (b). The magnetic field generated by the coil differs slightly depending on the bundle position.

Figure 10 :

 10 Figure 10 : Phase shift between central and peripheral bundle currents as a function of frequency.In order to compare these theoretical results with reality, the current distribution of two different Litz wires was measured. The results are analysed in the following section.

Table 1 :

 1 Simulated current distribution in a 5-bundles Litz wire. The differences with the total current divided by the number of bundles is given.

		T1	T2	T3	T4	T5
	Current	7.6 A	6.2 A	6.7 A	8.5 A	9.1 A
	Difference 0.7 % 18.6 % 12.1 % 11.7 % 19.7 %

Table 2 :

 2 Comparison of the measured and simulated current distribution of Litz wire n°2 (5x5x51) at 82 kHz.

	Bundle	Measures	Simulation	Error
	T1	4.77 𝐴	4.43 𝐴	7.2 %
	T2	6.19 𝐴	6.40 𝐴	1.5 %
	T3	4.77 𝐴	5.21 𝐴	9.2 %
	T4	5.48 𝐴	5.40 𝐴	1.5 %
	T5	5.80 𝐴	5.86 𝐴	1.1 %
	Significant phase shifts were measured at 82 kHz in 3 of
	the 16 bundles of Litz wire n°1

Table 3 :

 3 Phase shift measured as a function of frequency of bundle n°14 of a 16-bundles Litz wire.

	Frequency	0.1 kHz	1 kHz	2.6 kHz 82 kHz
	Phase shift	3.1°	19.5°	47.9°	165°


Kimon Keramidas

Silvana Mima

Adrien Bidaud

Estimation approches for materials

come

> RECYCLE:

-with recycling strategies, steel scrap price ,  scrap competitiveness with climate policies: electricity price ,  electric arc competitiveness -rec+ scenario with increased steel scrap recovery factor (75%100% of recoverable; 25% considered unrecoverable)  scrap available rises to 60-80% of annual demand  EAF covers 50-70% of steel production

PROJECTIONS of materials demand

Estimation approaches for materials demand, recycling and substitution using the POLES model Kimon KERAMIDAS, Silvana MIMA, Adrien BIDAUD (Université Grenoble-Alpes)

kimon.keramidas@univ-grenoble-alpes.fr -Building codes: cement use can be cut by as much as half without losses in structural integrity -Material substitution: timber products can partially substitute concrete (as much as 80% for multi-storey buildings) and even steel; glass can cover a larger share of outer shell > REUSE:

-Design with end-of-life in mind: wall panels can be reused instead of sent to demolition landfill -Circularity beyond CO 2 considerations: concrete can be ground to substitute gravel for new concrete  Large potential for material efficiency still untapped

  Scenarios modelled with POLES, derived from JRC GECO 2021 https://ec.europa.eu/jrc/geco | All results from 1.5°C scenario: world carbon price logistic over 2021-2100, at 1200 USD/tCO 2 in 2050 MAIN DRIVERS of materials demand FURTHER RESEARCH -Quantify material consumption related to urban planning design, land use patterns, buildings types -Further justify country differences on specific material consumption levels -Estimate material substitution trade-offs WORK IN PROGRESS findings not final + TRADE-OFFS on land-use -PV needs for hydrogen electrolysis for steel direct reduction: x14 land needs compared to a current steel plant -Timber instead of cement in buildings construction (lead to 1/3 less cement demand): plantation surfaces for timber up to x4 vs present In addition, not modelled above: > REDUCE:



  All drivers push towards more materials demand

	METHODS for materials demand
	Method	Pros	Cons	Method	Explicitly modelled	Implicit/not modelled
	Econo-	Based on readily	Resulting stock	Buildings	m 2 /cap, cap/dwelling, shift from	construction codes, building
	metric,	available data on annual	differences across		informal to modern materials,	typology, cities density, real estate
	Intensity	demand	countries unexplained;		renewal & renovation rates	cyclical activity, secondary houses,
	of Use		no end-use detail				domestic migration
	Intensity of Stock	Annual demand coherent with provision	Stock data recreated from annual demand	Road transport	car ownership ratio, paved network size	cities density, population distribution
		of services	data; no end-use detail
	Bottom-	Annual demand	Resulting stock	Power gen	technologies mix	sources distance to grid network
	up	coherent with provision	differences across
		of services	countries unexplained		Buildings	Road vehicles	Road network	Power gen.
			This study
	World steel demand	World cement demand
		Resulting material stock

+ main production-side mitigation options: electrification, CCS, biomass, hydrogen


J.-M Gaudry

J.-M Veillé

Férec

M Riou

G Ausias

Y Grohens

T Gaudry

J Férec

Thermoplastic foaming with thermo-expandable microcapsules: Mathematical modeling and numerical simulation for extrusion process

Keywords: Thermo-expandable microcapsules, Growth model, Thermoplastic, Foaming, Extrusion, Numerical simulation, Die swell

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, foam materials are competitive materials in industry market and may fulfill some product requirements. Although introduced by Dow Chemical ® in the early 1970s [START_REF] Morehouse | Expansible thermoplastic polymer particles containing volatile fluid foaming agent and method of foaming the same[END_REF], Thermo-Expandable Microcapsules (TEM) have recently gained in interest due to their performance and affordability. A TEM is made of a polymeric shell which encapsulates saturated hydrocarbons. Upon heating, the shell softens while internal pressure increases, leading to blow the microcapsule. Compared to other foaming technologies, microstructures created by TEM do not involve cell nucleation and gas diffusion phenomena. It results in a much more controllable process and product quality. In fact, the cell density is directly managed by the dozing system and the cell size is tailored by TEM composition and process parameters. Furthermore, TEM foaming technology can be performed on every standard extrusion system. However, a more extensive understanding of the foaming behavior during extrusion is necessary to optimize and prevent possible issues during industrial production.

Bubble growth and cavitation problems have been studied for one century. Rayleigh [START_REF] Rayleigh | on the pressure developed in a liquid during the collapse of a spherical cavity[END_REF] carried out the pioneer work by proposing a simplified bubble radius evolution into an incompressible inviscid liquid. Later, Plesset [START_REF] Plesset | The dynamics of cavitation bubbles[END_REF] generalized this relationship by taking into account the surface tension, the liquid viscosity and a non-constant internal versus external pressure difference. Then, Barlow et al. [START_REF] Barlow | Diffusion of gas from a liquid into an expanding bubble[END_REF] combined mass and momentum transfer equations to simulate diffusion phenomena during bubble growth. Street et al. [START_REF] Street | Dynamics of phase growth in viscous, non-newtonian liquids. initial stages of growth[END_REF] modeled a single bubble growth in a non-Newtonian fluid, whereas Patel [START_REF] Patel | Bubble growth in a viscous newtonian liquid[END_REF] simplified the single bubble growth in a Newtonian fluid by taking into account Rosner and Epstein assumptions on concentration profile [START_REF] Rosner | Effects of interface kinetics, capillarity and solute diffusion on bubble growth rates in highly supersaturated liquids[END_REF]. Amon and Denson [START_REF] Amon | A study of the dynamics of foam growth: Analysis of the growth of closely spaced spherical bubbles[END_REF] introduced a polymer shell in their model to consider the limited amount of available gas when multiple bubbles are growing simultaneously. Arefmanesh and Advani [START_REF] Arefmanesh | Diffusion-induced growth of a gas bubble in a viscoelastic fluid[END_REF][START_REF] Arefmanesh | Nonisothermal bubble growth in polymeric foams[END_REF] proposed a complete model based on mass, momentum and energy transfer equations. Furthermore, an upper convected Maxwell (UCM) model was used to describe the bubble growth in a viscoelastic fluid. Recently, Kawaguchi et al. [START_REF] Kawaguchi | Development of thermally expandable microcapsule and their mathematical models for polymer foaming[END_REF] and then Fujino et al. [START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF] developed models for TEM growth surrounded by air and by a polymer. Up to now, the later work seems to be the most advanced model. They considered mass and transfer equations including diffusion, and used the UCM constitutive equation to represent shell and fluid viscoelasticity.

Only few papers deal with polymer foaming phenomenon during continuous process. The complexity of blowing mechanisms restricts the number of contributions. Ranganathan and Advani [START_REF] Ranganathan | A non-isothermal process model for consolidation and void reduction during in-situ tow placement of thermoplastic composites[END_REF] modeled a composite, which contains voids in compression molding. Combining Patel growth model [START_REF] Patel | Bubble growth in a viscous newtonian liquid[END_REF] and Baldwin fluid model [START_REF] Baldwin | Microcellular polymer processing and the design of a continuous sheet processing system[END_REF], Tsujimura et al. [START_REF] Tsujimura | A study of bubble nucleation in foam extrusion die[END_REF] proposed a way to estimate the onset position of foam formation in the extrusion die for a PP-isobutane system. Shafi et al. [START_REF] Shafi | Bubble size distributions in freely expanded polymer foams[END_REF] and Joshi et al. [START_REF] Joshi | Prediction of cellular structure in free expansion of viscoelastic media[END_REF] described the nucleation and growth dynamics for batch foaming systems. Shimoda [START_REF] Shimoda | Polymeric foaming simulation for extrusion processes[END_REF] adapted Shafi et al. model to the extrusion process. Recently, Uchio et al. [START_REF] Uchio | Effect of viscoelastic properties of the base polymer on extrusion foaming with thermally expandable microcapsules[END_REF] investigates the TEM foaming phenomenon by extrusion for different PP matrices. For low shear rates and highly viscous polymers they found that expansion can be restricted.

Otherwise, extrudate swell problem for Newtonian fluid is well-known [START_REF] Park | Birefringence studies in die flows of an HDPE melt[END_REF][START_REF] Georgiou | The compressible newtonian extrudate swell problem[END_REF][START_REF] Ganvir | Prediction of extrudate swell in polymer melt extrusion using an arbitrary lagrangian eulerian (ALE) based finite element method[END_REF][START_REF] Choi | Simulation of extrudate swell using an extended finite element method[END_REF]. Die swell is a physical phenomenon due to the velocity profile rearrangement of a polymer flow subjected to shear stress. Investigations of non-isothermal [START_REF] Konaganti | Nonisothermal extrudate swell[END_REF] and compressible flows [START_REF] Taliadorou | Numerical simulation of the extrusion of strongly compressible newtonian liquids[END_REF][START_REF] Mitsoulis | A study of various factors affecting newtonian extrudate swell[END_REF] have been treated. For low Reynolds numbers, Talidorou et al. [START_REF] Taliadorou | Numerical simulation of the extrusion of strongly compressible newtonian liquids[END_REF] and then Mitsoulis et al. [START_REF] Mitsoulis | A study of various factors affecting newtonian extrudate swell[END_REF] showed that strand swelling is not significantly affected by a highly compressible flow, in the sense of pressure dependent. Recently, 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] studied and compared 2D and 3D swelling resolution with a PTT model for different width/height die ratio. For high width/height die ratio, they found less deviations between 2D and 3D die swell simulation due to a larger axial velocity than for a low width/height ratio.

In this paper, a microscopic Newtonian thermo-expandable microcapsule growth model is combined with a macroscopic Newtonian flow field simulation. Non-isothermal flow and die swell are considered to estimate the foaming influence on process and material properties. The flow is assumed compressible in the sense of foam dependent. The model is detailed in section 2. The numerical implementation is explained in section 3. Section 4 presents methods and material parameters used in this work. Microcapsule shell viscosity tuning is described in section 4. Die swell, temperature and matrix viscosity influences on foaming flow characteristics are exposed and analyzed in section 5.

TEM thermoplastic extrusion foaming model

The aim of this model is to represent the foaming process of TEM embedded inside a thermoplastic matrix during an extrusion process. The global model can be divided into two separate models.

Microscopic model for a TEM expansion

This model describes the TEM growth at the microscopic scale in a surrounding infinite medium.

Model statements

Fig. 1a shows the top view of a TEM obtained by optical microscopy. A cross section view obtained by scanning electron microscopy is illustrated in Fig. 1b. An outer diameter of 13.4 µm and a shell thickness of 3.0 µm can be measured. The TEM growth is found to contain three phases: an outer polymeric shell, an inner liquid phase and an inner gaseous phase. A schematic representation for an embedded TEM is shown in Fig. 2, where the infinite surrounding fluid is either air or a polymer. For the later case, it is worth mentioning noting that the surrounding medium and the TEM shell are made of two different polymers. The microcapsule growth is driven by the internal pressure, which is related to the nature of the blowing agent and the microcapsule dimensions. Subjected to elevated temperatures, the internal pressure increases whereas the shell softens. Once the internal pressure is high enough to counteract the external pressure and the shell resistance, the expansion phenomenon occurs. Several assumptions are required to describe the TEM expansion. Most of them are usual assumptions proposed by [START_REF] Kawaguchi | Development of thermally expandable microcapsule and their mathematical models for polymer foaming[END_REF][START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF]. The model considers one microcapsule with an average behavior, which represents the real expansion for a batch of microparticles inside an infinite fluid. h1: A TEM is assumed to be a concentric sphere of external radius R ext and internal radius R int , and remains concentric during the expansion. Therefore, a microcapsule grows uniformly along the radial direction r and the velocity vector expressed in a spherical coordinates (r, θ , ϕ) simplifies to V = V r (r,t) e r , where t denotes the time. The interaction between neighboring TEM is neglected.

h2: Gravity and inertial effects are neglected since highly viscous polymers and small particles imply low Reynolds numbers.

h3: Both the shell material (index sh) and the surrounding polymer matrix (index pm) are supposed to be incompressible Newtonian fluids.

h4: ρ sh represents the shell mass density at room temperature, which is assumed to be temperature independent. The shell material is not damaged during the expansion. h5: The surface tension γ is neglected. Indeed, shell interfacial tension orderof-magnitude is around γ ∼ 10 -2 N/m [START_REF] Owens | Estimation of the surface free energy of polymers[END_REF][START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF] and therefore 2γ/R int and 2γ/R ext terms are around 10 3 Pa, whereas pressure difference is around 10 5 to 10 6 Pa.

h6: The gas diffusion through the shell is not considered.

h7: The blowing agents are miscible in the liquid phase and behave as an ideal gas in the gaseous phase.

h8: The heat flow within the TEM is fast enough to consider the same temperature as the surrounding fluid, and therefore, expansion is assumed isothermal.

Blowing dynamics is ruled by the continuity and the momentum balance equations. Following hypothesis h2 and h3, these governing equations reduced to

∇ • V = 0, (1)
-∇P + ∇ • σ = 0, ( 2
)
where P is the pressure, σ = 2η ε is the deviatoric stress tensor, η the Newtonian viscosity of a material and ε is the strain rate tensor. In term of spherical coordinates (r, θ , ϕ), the continuity equation, Eq. (1), yields after simplification

∂ ∂ r r 2 V r (r,t) = 0, (3)
which can be rearranged to give after integration

V r (r,t) = F(t) r 2 , (4)
where F(t) is a r-independent function. The inner velocity V r (R int ,t) ≡ Ṙint and the outer velocity at the shell interface V r (R ext ,t) ≡ Ṙext can be expressed in function of F(t) such as

F(t) = Ṙint R 2 int = Ṙext R 2 ext . (5)
From assumption h1, the deviatoric stress tensor becomes a diagonal tensor with components σ rr and σ θ θ = σ ϕϕ (for symmetry condition). Therefore, Eq. ( 2) can be written as

- ∂ P ∂ r + ∂ σ rr ∂ r + 2 (σ rr -σ θ θ ) r = 0 . (6)

TEM expansion in polymer

As suggested by Kawaguchi et al. [START_REF] Kawaguchi | Development of thermally expandable microcapsule and their mathematical models for polymer foaming[END_REF] and Fujino et al. [START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF] to tackle the problem of the microcapsule expansion in a polymer matrix, Eq. ( 6) is integrated between r = R int and r = R ext for the shell domain (index sh) and between r = R ext and r → ∞ for the surrounding medium (index pm). Thanks to the Chasles theorem, one obtains -P(R ext ) + P(R int ) + σ rr,sh (R ext )σ rr,sh (R int ) + ∆σ sh -P(R ∞ ) + P(R ext ) + σ rr, pm (R ∞ )σ rr,pm (R ext ) + ∆σ pm = 0, [START_REF] Rosner | Effects of interface kinetics, capillarity and solute diffusion on bubble growth rates in highly supersaturated liquids[END_REF] where

∆σ sh = 2 R ext R int σ rr,sh -σ θ θ ,sh r dr, (8)
and

∆σ pm = 2 ∞ R ext σ rr,pm -σ θ θ ,pm r dr. (9)
The stress continuity condition at the inner microcapsule surface requires that P(R int )σ rr,sh (R int ) = P int as the surface tension is neglected (h5). At the interface between the outer microcapsule surface and the polymer matrix, a force balance leads to -P(R ext )+σ rr,sh (R ext ) = +P(R ext )-σ rr,pm (R ext ). The condition of stress continuity at the far-field boundary is P(R ∞ )σ rr,pm (R ∞ ) = P ∞ , where P ∞ is the bulk pressure. In the case of extrusion, P ∞ represents the pressure inside the canal. Also, at an infinite distance from the TEM, the stress σ rr, pm (R ∞ ) can be neglected since the polymer deformation is close to zero. Thus, Eq. ( 7) simplifies to

P int -P ∞ + ∆σ sh + ∆σ pm = 0. ( 10
)
∆σ sh and ∆σ pm are obtained from their rheological models. Since a Newtonian behavior is considered for both the shell and the thermoplastic polymer surrounding the TEM (assumption h3), their expressions for the normal stress components are σ rr = -2σ θ θ = -4ηF(t)/r 3 . After substituting the relationship for F(t), the integration of Eq. ( 8) is performed to yield

∆σ sh = -4η sh Ṙint R int R 3 ext -R 3 int R 3 ext . ( 11
)
On the same way, one obtains the polymer matrix contribution from Eq. ( 9)

∆σ pm = -4η pm Ṙext R ext = -4η pm Ṙint R 2 int R 3 ext . (12)
When Eqs. [START_REF] Kawaguchi | Development of thermally expandable microcapsule and their mathematical models for polymer foaming[END_REF] and [START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF] are inserted in Eq. ( 10), a governing equation for the growth dynamics of a TEM made of a Newtonian shell material inside an infinite Newtonian thermoplastic matrix is

Ṙint = R int (P int -P ∞ ) 4 η sh R 3 ext -R 3 int R 3 ext + η pm R 3 int R 3 ext . ( 13
)
It is worth mentioning that the upper dot symbol means the total time derivative. When air is considered as the suspending fluid, η pm = 0 and P ∞ = P atm with P atm the atmospheric pressure, one recovers almost the dynamics growth equation for a bubble surrounded by a shell in air Ṙint = R int (P int -P atm )

4η sh R 3 ext -R 3 int R 3 ext . ( 14
)
This equation is equivalent to Eq. (57a) presented by Elshereff et al. [START_REF] Elshereef | Comparison and analysis of bubble growth and foam formation models[END_REF], and similar to Eq. ( 12) in terms of the volume of the shell of Amon et al. work [START_REF] Amon | A study of the dynamics of foam growth: Analysis of the growth of closely spaced spherical bubbles[END_REF], when the interfacial tension is neglected.

TEM Thermodynamics

According to assumption h7, the blowing agent mixture follows an ideal behavior, therefore each component participates to the global internal pressure. The internal pressure P int is given by the Raoult law such as P int = ∑ P k , where P k is the partial pressure of the k blowing agents. Since the hydrocarbons are mostly in the liquid state at the initial state, the pressure inside a close cavity is determined by the saturated vapor pressure equation. For an ideal mixture, each component contributes to the internal pressure up to their relative partial pressure, which is given for the liquid phase by

P k = P sat,k X M,k . (15)
In the above equation, X M,k is the molar fraction of for the k th component and P sat,k is its vapor pressure. The saturated vapor pressures for the different components are obtained by the Antoine equation [START_REF] Antoine | Tensions des vapeurs : nouvelle relation entre les tensions et les températures[END_REF]

P sat,k = exp A k - B k T +C k , (16)
where A k , B k and C k are Antoine equation constants. The partial pressure can be also expressed in function of the vapor phase by the ideal gas law

P k = m vap,k R T M k V vap , (17)
where R = 8.314462 J/(mol • K) [31] is the ideal gas law constant, M k represents the molar mass of the k th component, V vap and m vap,k denotes the whole volume and the mass of the k th component for the vapor phase, respectively. When there is presence of liquid phase, a liquid/vapor equilibrium exists, thus Eqs. ( 15) and ( 17) are equal. Therefore, m vap,k can be estimated from this equivalence.

m vap,k = P sat,k X M,k M k V vap R T . ( 18
)
As the diffusion phenomenon is neglected (hypothesis h6), the balance for a blowing agent material is given by

m k = m liq,k + m vap,k , (19)
where m k and m liq,k are the total mass and liquid mass of the k th component, respectively. The density of the hydrocarbon mixture is defined by their individual density weighted by their own volume fraction such as

ρ hc = k ∑ i=1 ρ liq,k X M,k . (20)
In the above equation, ρ liq,k is the density of the k th component. Initially, most of the mixture is in the liquid state, the following approximation can be made

m k | t=0 ≈ k ∑ i=1 m liq,k | t=0 , (21)
where m k | t=0 and m liq,k | t=0 are the initial total mass and the initial liquid mass of the k th component at the initial state, respectively. Eqs. ( 20) and ( 21) are used to estimate the liquid volume V liq and therefore to approximate the gaseous phase volume V vap with the following conservation relationship V int = V liq +V vap , where V int the internal volume of a TEM. By solving iteratively Eq. ( 19) using Eqs. ( 20) and ( 21), m liq,k and m vap,k can be determined at each time.

Blowing dynamics resolution

The total volume of a TEM can be defined as

V T EM = V int +V sh , (22)
where

V T EM = 4πR 3 ext /3, V int = 4πR 3 int /3 and V sh = 4π(R 3 ext -R 3
int )/3 are the external, internal and shell volume of the TEM, respectively. If R ext | t=0 and R int | t=0 are the initial external and internal radii of a TEM, R ext is (by assuming that the volume of the shell V sh remains constant)

R ext = 3 R 3 ext | t=0 + R 3 int -R 3 int | t=0 . (23)
The initial contents of hydrocarbons can be evaluated by solving Eq. ( 19) with the help of Eqs. ( 18) and [START_REF] Park | Birefringence studies in die flows of an HDPE melt[END_REF][START_REF] Georgiou | The compressible newtonian extrudate swell problem[END_REF][START_REF] Ganvir | Prediction of extrudate swell in polymer melt extrusion using an arbitrary lagrangian eulerian (ALE) based finite element method[END_REF]. Once the internal pressure is known, the growth dynamics, Eq. ( 13), is used to determine the expansion velocity. The internal radius R int is evaluated iteratively to predict the expansion under specific system conditions defined by P ext , T , η sh and η pm . When the liquid is fully evaporated, the internal pressure is only given by the ideal gas law Eq. ( 17). As the microcapsule continues to grow, the internal pressure begins to decrease until an equilibrium is reached depending on system conditions.

Macroscopic flow model for a syntactic foam

This section aims to simulate the TEM blowing phenomenon within a thermoplastic for a steady-state flow. A set of hypotheses is addressed to simulate the flow for this syntactic foam (specified by subscript "s f ") H1: The syntactic foam material is considered as an isotropic material. H2: Gravity and inertial effects are neglected. Indeed, dealing with highly viscous polymer and small particles implies a low Reynolds number.

H3: The density of the suspending matrices is presumed to be temperature independent.

H4: The polymer matrice is assumed to be a slightly compressible Newtonian fluid. The density of the syntactic foam (ρ s f ) is dependent on the TEM expansion (R ext ). Details are presented in Appendix B.

H5: The viscosity of the fluid is assumed to be not affected by the evolution of ρ s f . H6: As the Mach number is low for the flow, the work done by pressure is neglected in the energy equation.

H7: Heat capacity (C p ) and thermal conductivity (λ T ) are assumed constant. Neither temperature nor foaming affect both parameters. Regarding the porosity dependency, the extrusion speed is relatively high enough for the outer boundary length and for the low thermal conductivity of polymer melt, to limit the effect of the foaming phenomena on the temperature profile induce by conduction in y-direction.

H8: This set of equations gives the quasi-stationary solution.

Following these hypothesis, the governing equations for a non-isothermal flow for a syntactic foam is given by

∇ • ρ s f u = 0, (24)
-∇P + ∇ • σ s f = 0, (25)
σ s f = 2η app ε - 2 3 η app (∇ • u) I, (26)
ρ s f C p ( u • ∇T ) -λ T ∆T = σ s f : ∇ u, (27)
ρ s f = X v,T EM | t=0 • R ext R ext | t=0 3 • (ρ T EM -ρ pm ) + ρ pm , ( 28
)
where u is the fluid velocity vector in a 2D planar geometry with u and v the velocity components in the xand y-directions, respectively. I is the identity tensor and η app is the apparent shear viscosity determined experimentally through a rheometric slit die. The last term in Eq. ( 27) describes the heat flux arising from the viscous dissipation. Eq. ( 28) clearly shows the dependency of syntactic foam density ρ s f on the TEM size evolution R ext .

Numerical implementation

Geometry and boundary conditions

A 2D geometry with a Cartesian coordinate system (x, y) is used to examine the melt flow in a slit die. The channel geometry, shown in Fig. 3, is split in two domains. The left part represents the die, whereas the right one acts for the extrudate at the die exit. Non-isothermal conditions are applied to simulate auto-heating inside the die and cooling at the die exit. The length, the thickness and the width of the slit die are L die = 60 mm, t die = 2 mm and W die = 40 mm, respectively. In the real process, the extrudate enters inside a cold mechanical shaper in a cooled water bath after the die exit. Therefore, the foaming phenomenon is stopped at L ext = 50 mm.

Q m BC in BC out h(x,t) h die BC fs BC die BC sym L die L ext x y

Flow boundary conditions

A symmetric boundary condition BC sym is applied along the line at y = 0. As the result, the geometry displays only half thickness (h die = t die /2). The inlet boundary condition at BC in considers a fully developed flow carrying a fixed mass flow rate Q m =< u > W die h die ρ s f , where < u > is the average velocity of the melt flow. A no-slip condition, defined by u = 0, is applied on the die wall boundary (BC die ). A zero-pressure condition is imposed at the outlet boundary (BC out ), which provides for a zero fluid stress along this surface and the backflow is also suppressed. The BC f s condition consists in a free surface that must satisfy three conditions simultaneously. The conditions of zero shear and normal stresses are prescribed and the zero normal velocity condition is imposed by the following kinetic equation,

∂ h / ∂t + u ∂ h /∂ x = v [32]
, where h(x) the height of swelling along the x-direction. This kinematic condition defines the free surface as a streamline where each element fluid particle remains part of the free surface, and therefore, the particle velocity should always remains tangent to the free surface. One way to handle this issue is to use the arbitrary Lagrangian-Eulerian (ALE) formulation [START_REF] Hughes | Lagrangian-eulerian finite element formulation for incompressible viscous flows[END_REF][START_REF] Ganvir | Prediction of extrudate swell in polymer melt extrusion using an arbitrary lagrangian eulerian (ALE) based finite element method[END_REF][START_REF] Choi | Simulation of extrudate swell using an extended finite element method[END_REF]. This method consists in two steps in the sense that the Lagrangian mode allows the mesh to move with the material and then the mesh can freely move as required by the flow domain. The ALE formulation reduces the mesh distortion while tracking material flow. The kinetic equation is implemented in the COMSOL Multiphysics ® software using a coefficient form boundary PDE. An artificial diffusion coefficient of c num = 10 -5 m/s 2 is used to stabilize the computation.

Thermal and thermodynamics boundary conditions

In the extrusion process, a temperature T die is imposed at the BC die surface and along the symmetric axis of the die (BC sym ), the surface is insulated from the system. To establish the condition at the flow inlet (BC in ), the velocity component along the x-direction (i.e., u) for a Newtonian fluid in a steady state regime inside a planar slit die is defined such as

u = 3 2 < u > 1 - y h die 2 , (29)
The dissipated heat is given by

σ s f : ∇ u = 2η app ∂ u ∂ y 2 = 18 η app < u > 2 y 2 (h die ) 4 (30)
In case of polymer processing, the heat conduction in the flowing direction (x) is negligible compared to (ρ s f C p u ∂ T /∂ x) which corresponds to the convection term [START_REF] Agassant | Polymer Processing 2E: Principles and Modeling[END_REF]. Then, for an equilibrium regime the convection term becomes equal to zero. Thus, after integrating Eq. ( 27) the fully developed temperature profile at the flow inlet (BC in ) is prescribed [START_REF] Agassant | Polymer Processing 2E: Principles and Modeling[END_REF] T

in = T die + 3η app < u > 2 2λ T 1 - y h die 4 . (31)
A free convective heat transfer is assumed to the BC f s such as q = h cv • (T amb -T die ), where T amb the ambient temperature and h cv the convective coefficient for air convection.

Finally, a global constraint condition is applied on the internal pressure inside a TEM to respect the fact that the liquid phase herein is incompressible

P int = P ∞ for P ∞ > k ∑ i=1 P k , P int = k ∑ i=1 Pk for P ∞ < k ∑ i=1 P k . ( 32
)
As a result, microcapsules are assumed to not collapse under high external pressure.

Solution strategy

The extrusion foaming simulation is illustrated by the flow chart in Fig. 4. A stationary flow solution is first used to determine the velocity and the pressure fields. Then, the macroscopic and microscopic models are simultaneously solved. On the one hand, the ALE formulation computes the free surface whereas the non-isothermal problem Eq. ( 27) gives the temperature field distribution. On the other hand, the microscopic blowing model predicts the TEM foaming, which modifies the density and therefore the velocity field. It consists in solving Eq. ( 13) which is implemented in the software through a convection-diffusion PDE type equation and gives the radius R int and the internal volume V int at each position. Eq. ( 19) is solved using an algebraic differential equation to converge toward a blowing agent state composition and P int is finally determined. An additional step can be appended to stabilize the computation consisting in solving a time-dependent multiphysics model (free surface and non-isothermal flow), before solving simultaneously the macroscopic and microscopic models. The time-dependent solver uses the Backward Differentiation Formula (BDF) algorithm to control the time step and employs a segregated solution approach to solve the system of equations. A non-linear Backward Euler method is applied to approximate the time derivative and each step represents one physics. The first two steps deal with the flow and thermal solutions and then the third step solves the blowing equations Eqs. ( 13) and [START_REF] Uchio | Effect of viscoelastic properties of the base polymer on extrusion foaming with thermally expandable microcapsules[END_REF].

All the computations are carried out on a Dell Precision Tower 7810 workstation, equipped with an Intel ® Core ™ Intel ™ Xeon ™ CPU E5-2620 v4 at 2.1 GHz and 128 Go of RAM. The fluid discretization employed for velocity and pressure field is P1 + P1. The default mesh consists of coarse size elements and 13

Start

Flow input parameters

 app Q m  pm W t Ldie Lext Pext u v
Blowing input parameters R ext (0) t sh (0)  TEM  pm TEM (0)sh Xm,hc Xm,TEM

Stationary flow resolution

Thermal input parameters Step 2

T ext h cv A k B k C k
Step 3

Step 1

Figure 4: Simulation flow chart is generated and adapted by the software following the physics employed in the simulation (i.e., the mesh close to the no-slip boundary layer, the edges and free surface is made of quadrangle elements, in order to keep an acceptable precision facing strong gradients). In total, the geometry surface is meshed with 10167 elements.

Materials and methods

Thermoplastic matrices

A vulcanized thermoplastic elastomer (TPV) and a polypropylene (PP) matrices have been considered as suspending fluids in this study. Matrices have a different rheological behavior. The thermal conductivity is measured by the hot wire method following EN ISO 22007-1:2017 [START_REF]AFNOR EN ISO 22007-1, Plastics -determination of thermal conductivity and thermal diffusivity -part 1 : general principles[END_REF] using Rheograph 20 supplied by Göttfert. The density at room temperature is determined by an AlfaMirage MDS-300 electronic densimeter following ISO standard 1183-1 [START_REF] Afnor Nf | Plastics -methods for determining the density of non-cellular plastics -part 1 : immersion method, liquid pycnometer method and titration method[END_REF]. PVT measurements are done to determine the melt density (ρ T ) of both thermoplastics in the process temperature range (180-220°C). A capillary rheometer Rheograph 20 is also used to perform rheological, thermal conductivity (λ T ) and PVT measurements. Firstly, a Bagley correction are used to obtain the apparent viscosity η app , then a Rabinowitsch correction considering a power-law shear thinning behavior [START_REF] Macosko | Rheology: Principles, Measurements, and Applications[END_REF] is applied to determine the true viscosity η pm . In this simulation, because the Newtonian model is considered, η app is applied in the flow computations to obtain the correct pressure profile and η pm represents the bulk viscosity seen by TEM during expansion. Thermal conductivity (λ T ), viscosities (η app and η pm ) and melt density (ρ T ) are slightly affected by temperature variations in the range of temperature extrusion (180-220°C). Each characteristic of both polymer matrices is approximated by one average value and Table 1 summarizes the data parameters used for the simulation.

Thermo-expandable microcapsules (TEM)

A TEM grade has been selected to match temperatures during an extrusion. Expansion behavior has been measured using a hot plate and an optical microscope coupled with a high speed camera. A temperature heating slope of 10 °C/min has been set, ensuring a long time period at elevated temperature, which is roughly equivalent to the residence time of TEM in an extrusion unit. The acquired images are binarized and then analyzed to determine the average TEM behavior in terms of temperature. Fig. 5 depicts the average expansion ratio λ ext = R ext /R ext | t=0 . The expansion curve consists in three differents step. An abrupt expansion until the liquid phase is evaporated (see Appendix A), then a gently growth due to the pressure increasing along the temperature and a shrinkage due to gas leaking through the TEM shell. To simulate TEM expansion in air, Eq. ( 13) is used to approximate the bubble growth behavior. Since the viscosity measurements η sh for the pure shell material of the TEM is tricky to perform, η sh is obtained by the inverse method with experimental data. The expression for η sh is reported in Table 2. The form is similar to the Kawaguchi's et al. model [START_REF] Kawaguchi | Development of thermally expandable microcapsule and their mathematical models for polymer foaming[END_REF], which considers an Arrhenius dependance of the viscosity with the temperature. The starting temperature and the starting expansion slope are fitted as close as possible to the experimental data and the result (called "Original fit") is shown in Fig. 5. As expected, the model prediction exhibits a higher expansion level attributed to the absence of diffusion (see Fujino et al. [START_REF] Fujino | Mathematical models and numerical simulations of a thermally expandable microballoon for plastic foaming[END_REF]) and to the choice of a Newtonian behavior for the TEM shell. More accuracy is obtained by adjusting the blowing agent amount, which reduces the liquid reservoir and affects P int when only vapor phase is remaining, as represented by the curve called "Final fit" is shown in Fig. 5.

Q m ρ T C p λ T [Pa • s] [Pa • s] [kg/h] [g/cm 3 ] [J/(kg • K)] [W/(m 2 • K)] TPV 30
The effect of viscosity variations on the TEM growth kinetics is shown by Fig. 6. These curves have been computed by considering T = 205 °C as an initial condition. It represents the extreme case where initially P ∞ > P int , and then P ∞ O r i g i n a l f i t F i n a l f i t H e a t i n g r a t e : 1 0 °C / m i n decreases infinitively fast until P atm . Fig. 6 clearly indicates the major role of viscosity in expansion kinetics as shown by Kawaguchi et al. [START_REF] Kawaguchi | Thermally expandable microcapsules for polymer foaming-relationship between expandability and viscoelasticity[END_REF]. When the viscosity is higher the growth is slower and inversely. This shift in time scale can be converted into a shift in temperature behavior. In fact, higher viscosities will apply more resistance to counteract the internal pressure and lead to expansion at higher temperatures. The initial microcapsule dimensions are measured by laser diffraction and microscopy techniques. R ext | t=0 corresponds to the initial average external radius of TEM and t sh | t=0 is the initial average thickness of TEM. The density ρ sh is determined after isolating the polymer shell material. It is done by using N,N-Dimethylformamide solvent to remove blowing agents and then, by evaporating it. Since ρ sh and ρ liq are known, the initial density of a TEM computed from volume considerations is estimated at ρ T EM | t=0 = 0.82 g/cm 3 (see Appendix B).

Gas contents are measured by gas chromatography coupled with a mass spectroscopy. Pentane, isopentane and isooctane are the three gases detected and their properties are reported in Table 3.

For instance the initial state at T = 478,15 K, computed for the blue curve in Fig. 6 is presented in Table 4. m liq and m vap stand for the whole amount in the a Experimental data b Antoine parameters are in SI system units (Pa and K) [31] liquid and vapor phase, respectively. As the blowing agent exists in the liquid phase at the initial state, the internal pressure (P int ) can be estimated by using Eqs. ( 15) and ( 16) with the parameters given in Table 3.

Syntactic foam extrusion

The extrusion system is an extruder equipped by a barrier screw (D=60 mm, full length L/D=26) and a dosing system. A slit die with the channel dimensions of L die = 60 mm, W die = 40 mm and t die = 2 mm is used. The cooling bath is located at L ext = 50 mm from the die exit and after this distance the syntactic foam microstructure is considered as fixed. The TEM mass fraction is imposed to X m,T EM = 0.65 %. Temperature profile is adjusted to avoid premature expansion inside the extrusion unit. In practice, a temperature at ± 2 °C is measured manually over 30 s with an infrared pyrometer, at the center of the upper extrudate surface as close as possible to the die exit. Polymer foam melt characteristics are determined experimentally during extrusion process from one pressure transducer located at 10 mm after the slit die entrance. The pressure transducer uncertain- ty is ± 0.5 bar. Density measurements and optical microscopy assessment are done after sample cooling. An optical microscopy image from a cross-section part of an extruded sample for a TPV syntactic foam is illustrated in Fig. 7. A quite good dispersion is observed for the TEM size distribution. The light reflection dots at the center of TEM reveal the presence of remaining shells. The external radius determined by image processing is equal to R ext = 20.0 ± 0.4 µm.

Die swell validation for a Newtonian fluid

By computing only step 1 and step 2 (see Fig. 4) or applying X m,T EM = 0 %, the extrudate swelling for a Newtonian fluid in both isothermal and nonisothermal conditions by considering the TEM are investigated. A part of the mesh and the final shape of the free surface are presented in Fig. 8. In both cases, the die swell ratio is found to equal 1.186 (see also Fig 16), which is in excellent agreement with the literature value for low Reynolds numbers [START_REF] Park | Birefringence studies in die flows of an HDPE melt[END_REF][START_REF] Ganvir | Prediction of extrudate swell in polymer melt extrusion using an arbitrary lagrangian eulerian (ALE) based finite element method[END_REF][START_REF] Mitsoulis | A study of various factors affecting newtonian extrudate swell[END_REF]. The fact that temperature gradient does not affect the die swell is due to the independency of suspending matrix density (ρ T ). Therefore, it implies the same velocity profile in both cases. Note that the swelling occurs in the first 2 mm after the die exit.

Thermoplastic foaming

Three probes PB in , PB end and PB out are created to evaluate some property values at the flow inlet, the die exit and the outlet problem, respectively. Furthermore, the start points of three streamlines (SL1, SL2 and SL3) are uniformly vertically distributed and are located at y/h die = 0.25, 0.50 and 0.75, as depicted in Fig. 8. For all the tested cases, the heat transfer coefficient h cv is set to 150 W/(m 2 • K) which corresponds to a forced-air convection, the ambient temperature is set to T amb = 25 °C and T die = 190 °C.

TPV foaming

The numerical results for the velocity and relative pressure fields and the temperature field are shown in Fig. 9 and Fig. 10, respectively. A linear pressure drop is observed inside the die, which is characteristic for a laminar flow in a slit die. White arrows in Fig. 9 represent the velocity field and are proportional to the velocity magnitude.

The temperature field in Fig. 10 exhibits the viscous effects which increase by 10 to 15 °C the temperature of the thermoplastic melt along the height. The average temperature is found to equal T = 205 °C at the die exit, and the forcedair convection cools the polymer melt to T = 202 °C at x = L ext .

Around 150 iterations are needed to reach convergence as shown in Fig. 11, where R ext is plotted versus iteration number. Results clearly indicate three different foaming levels for the three probe points. At the die entrance, R ext = 11.7 µm, which means the expansion has already started. It could be explained by the fact that pressure at the die entrance is lower than the vapor pressure inside the TEM (P int ≈ 32 bar). At the die exit, R ext = 13.8 µm whereas at x = L ext , R ext = 19.7 µm. Since the initial TEM radius is 9.4 µm, it means that 45% of the radius growing happens inside the extrusion line and 55% outside.

To explore the TEM growth, the evolution of the external TEM radius along three streamlines (i.e., along the x-direction) is plotted in Fig. 12. It is found that the external radius continues to increase if L ext > 50 mm. In fact at x = L ext , P int is still higher than P atm and therefore the microcapsules continue to grow. Nothing counteracts this expansion because the model does not take into account blowing agent diffusion through the shell, which reduces P int . In practice, the extrudate after x = L ext goes into a cooled mechanical shaper in a water cold bath. Consequently, the temperature decreases sharply and the expansion is stopped. That is why, R ext computed at x = L ext is assumed to be the final TEM radius.

Furthermore, almost no expansion variations are obtained between the three streamlines as shown in Fig. 12. In spite of the temperature gradient along the die height, created by the viscous dissipation, the TEM expansion is not tremendously affected. This is consistent with the low sensibility of the input TEM growth model at elevated temperature (> 190 °C). Additionally, the Newtonian behavior for the shell material leads to an uniform matrix viscosity in the y-direction. Following the simulated results, the foamed material would be ho- mogeneous. Optical microscopy observations on cross-section of extruded samples show no particle size gradient along the y-direction, which is in accordance with the simulation analysis illustrated in Fig. 12. As the growth variation in the y-direction is almost insignificant, only results along the streamline SL2 are used in the following discussions. Table 5 summarizes and compares the main flow characteristics with experimental process data. Overall, the predicted flow characteristics are found to be in reasonable agreement with the available experimental data as shown in Table 5. The relative pressure P ∞ and the temperature at the die exit T are well described and the porosity is slightly under-estimated. This can be explained experimentally by uncertainties in dosing the amount of TEM and in local porosity measurement. However, these results are acceptable. It is worth mentioning that at x = L ext , the average velocity is under-estimated by the simulation (< u >= 8.0 m/min). In practice, the extrusion speed is adjusted to counteract the die extrudate swell phenomenon in order to achieve the required part dimensions. The expansion ratio (λ ext ) and the porosity (φ ) are depicted in Fig. 13. λ ext at x = L ext is clearly limited to 2.1 when compared to 3.5 for air medium surrounding the TEM. In fact, for < u >= 8 m/min, it takes 0.375 for the flow to travel the 50 mm out of the die. At the end of the die R ext = 13.8 µm. Looking to the growth kinetics (see blue curve in Fig. 6), the time to reach R ext = 19.7 µm from 13.8 µm is equal to 0.442 s, which is in accordance with 0.375 s. It means the expansion is limited by the short time spent at high temperature and atmospheric pressure. Then, the porosity evolution is related to the macroscopic volume change due to the TEM expansion. It is found that an increase of 40% of R ext represents only 25% of the final porosity. Under these conditions, it could be concluded that the major part of the expansion occurs after the die exit.

Die swell investigation

The effect of the TEM expansion on the die swell is explored with the numerical model. Fig. 14 depicts the porosity rate in terms of the mass fraction of TEM with values ranging from 0 to 2%. Fig. 15 shows an increase of the average Figure 13: Expansion along the extrusion the SL2 streamline velocity in function of TEM foaming. This statement could be extended to the flow rate. The die swell increases with increasing the flow velocity as presented in Fig. 16, in which a higher die swell magnitude is observed with a larger TEM content. However, it has to be noticed that the relative increase of the die swell magnitude is lower than the total porosity. In fact, for X m,T EM = 2.00 %, φ = 18.3 % when die swell increasing is 8.5 %.

Temperature influence on thermoplastic foaming

The temperature is found to be a key parameter as it drives the TEM growth. Hence, two others temperatures (T die = 170 °C and T die = 180 °C) are tested to assess the predictability of the numerical model. The final average TEM radius along the height (R ext ) at x = L ext is compared with the mean experimental value measured by image processing in Fig. 17 and also an approximation for the TEM growth dynamics (see Fig. 5). Even if the model is not completely accurate quantitatively, qualitative analyses can be conducted.

PP foaming

For PP foaming simulation, the temperature profile is adapted to fit the same average temperature as for TPV foaming (i.e., 202°C) at x = L ext and therefore, T die = 174 °C. Although, this set temperature is lower than the one for the simulation with TPV, the average temperature remains almost the same at the die exit. In fact, since the PP is more viscous, it induces a higher viscous dissipation. This is in accordance with the extrusion experiments and from a practical point of view, a foaming reduction is expected due to this higher viscosity. When compared with TPV results, a higher pressure field is observed for the PP matrix.

Fig. 18 present the temperature along the SL2 streamline for both TPV and PP materials. The temperature is found higher for foaming the PP. This is due to the higher temperature gradient induces by the lower set temperature and the higher auto-heating.

Fig. 19 compare the TEM radius along the SL2 streamline for both TPV and PP materials. Curiously, the external radius of TEM in PP matrix (R ext = 19.9 µm) is close to the TPV foaming results (R ext = 19.7 µm). This simulation It can be concluded that the viscosity does not affect the expansion level. In practice, for PP foaming around 205 °C, the average external radius of TEM is equivalent to the external radius measured for TPV foaming. This confirms the model predictions, and therefore, the material viscosity appears to be insignificant regarding the final foaming state. These results are summarized in Table 6. Contrary to TPV foaming, the average velocity at the die exit (i.e., x = 0 mm) obtained by computations does not fit well with the experimental data. In fact, a strong die swell was experimentally observed for the PP foaming extrusion. The average velocity at x = L ext (i.e., < u >= 5.3 m/min) match better with the process data. It means that the Newtonian fluid die swell simulation is more consistent for PP foaming case than for TPV foaming. Furthermore, Fig. 19 shows a shift of the TEM expansion. In fact, because of the higher pressure due to the higher PP viscosity, the expansion of TEM only starts at the end of the die. It has to be noticed that the expansion for PP matrix occurs mainly after the die exit, consequently the expansion speed is higher. In fact, for < u >= 5.3 m/min the time to travel the 50 mm after the die exit is equal to 0,57 s whereas the time to grow from 10 µm to 19.2 µm is equal to 0,80 s following blue curve in Fig. 6. This is due to the higher pressure difference P int -P ∞ induced by the sudden P ∞ decreasing. Generally speaking, the expansion is faster for a more viscous polymer but catch up the same final level if enough time at high temperature and atmospheric pressure is let to the system.

Conclusion

In this work, thermoplastic foaming by TEM is investigated by developing a thermo-mechanical blowing model considering Newtonian behaviors. After a viscosity calibration from experimental data of TEM expansion in air, this microscopic growth model is coupled with a non-isothermal macroscopic flow simulation. Representing the real process conditions lead to consistent results compared to experimental measurements. Model predictions indicate that foaming occurs gently inside the die and then it continues after the die for low viscosity thermoplastic whereas for large viscosities the foaming phenomenon is shifted toward the die outside. The die swell phenomenon is promoted by the TEM expansion outside of the die. Overall, the foaming increases the extrusion flow rate when compared with the pure matrix extrusion.

Furthermore, results show that matrix viscosity has no significant influence on foaming for an equivalent temperature profile for the tested conditions. A TPV matrix with a relatively low viscosity produce a syntactic foam with the same porosity than a syntactic based PP matrix with a relatively high viscosity. The explanation comes from TEM shell nature, which is tremendously more viscous than thermoplastic melt and therefore, leads to similar expansion level. From an industrial point of view, material and tooling changes are not key parameters to pay attention in order to ensure a correct foaming. However, TEM are really dependent on temperature and the process temperature should be controlled carefully. It has to be noticed that for higher viscosity of the suspending fluid (i.e., with a more viscous polymer or lower shear rates) and/or lower shell viscosity, following the model the final foaming state can be affected, as experimentally shown by Uchio et al. [START_REF] Uchio | Effect of viscoelastic properties of the base polymer on extrusion foaming with thermally expandable microcapsules[END_REF]. Finally, the model suggests that the time at high temperature and atmospheric pressure after the die seems to be an important parameter to enable the TEM to reach their potential higher expansion level.

Despite the strong assumptions made for this model, predictions are in accordance with experimental results. This model can get more accurate by taking into account more realistic material behaviors such as a UCM constitutive equation for modeling the microcapsule growth, a non-Newtonian model to represent the shear-thinning behavior of flowing polymer including viscosity rise outside the die and the gas diffusion through the shell.

Acknowledgments

We warmly acknowledge Sekisui Chemical CO., LTD. for their support and their cooperation. mass density. From Eq. (B.5), the evolution of the syntactic foam density is finally given by

ρ s f = X v,T EM | t=0 • R ext R ext | t=0 3
• (ρ T EMρ pm ) + ρ pm , (B.7)

Figure 1 :

 1 Figure 1: Microscopic views of a TEM

Figure 3 :

 3 Figure 3: Die swell and foam extrusion model geometry

  p e r a t u r e [ °C ] E x p e r i m e n t a l E x p e r i m e n t a l ّ σ e x p

Figure 5 :

 5 Figure 5: Expansion ratio in term of temperature

Figure 6 :

 6 Figure 6: Viscosity influence on TEM growth kinetics.

Figure 7 :

 7 Figure 7: Cross-section part of an extruded sample

Figure 8 :Figure 9 :

 89 Figure 8: Deformed finite element mesh for the 2D slit die extrusion model

Figure 10 :

 10 Figure 10: Temperature field for the TPV foaming

Figure 11 :

 11 Figure 11: TPV foaming computational resolution

Figure 12 :

 12 Figure 12: Expansion ratio and porosity along the three streamlines

Figure 14 :

 14 Figure 14: Predicted porosity change along the SL2 streamline

Figure 15 :Figure 16 :

 1516 Figure 15: Velocity along streamline for different amount of TEM (data taken along the symmetry axis)

  Ex p e r i m e n t a l T P V S i m u l a t i o nA v e r a g e t e m p e r a t u r e a t t h e d i e e x i t ( T ) [ °C ] E x t e r n a l r a d i u s ( R e x t ) [ µ m ] I n i t i a l d i a m e t e r

Figure 17 :

 17 Figure 17: TEM radius in term of average temperature at the die exit

Figure 18 :Figure 19 :

 1819 Figure 18: Predicted temperature along the SL2 streamline for the two thermoplastics

Table 1 :

 1 Properties of the suspending thermoplastic matrices

	Material	η pm	η app

Table 2 :

 2 TEM shell properties

		ρ sh	a	ρ liq	R ext | t=0	a t sh | t=0	a	η sh
		[g/cm 3 ] [g/cm 3 ]		[µm]	[µm]	[Pa • s]
	TEM	1.20	0.64		9.4	2.4	exp (28000/T -44.3)
					a Experimental data
				Table 3: TEM blowing agent properties
		Substances k	X M,k [%]	a	M k [g/mol] A k [-] Antoine parameters b B k [-] C k [-]
		Pentane	1 51.5	72.15 20.73 2475.06 -39.94
		Isopentane 2 29.7	72.15 20.53 2348.66 -40.05
		Isooctane 3 18.8	114.23 20.58 2894.27 -52.41

Table 4 :

 4 Initial state at 205 °C Location x = 31.4 1.71 • 10 -13 0.74 • 10 -13 0.27 • 10 -15 1.17 • 10 -15

	Parameters	P int [bar]	m liq [kg]	m vap [kg]	V liq [m 3 ]	V vap [m 3 ]

Table 5 :

 5 Comparison between simulation and experimental process data for a TPV foaming

	Parameters	P ∞ [bar]	T [°C] [m/min] < u > Porosity [%]
	Location x = -50 mm 0 mm	0	+50 mm
	Experiment	22.8	205	9.4	10 ± 0.2
	Simulation	21.2	205	9.4	6.1

Table 6 :

 6 Experimental and predicted process data for a PP foaming

	Parameters	P ∞ [bar]	T [°C] [m/min] < u >	Porosity [%]
	Location x = -50 mm 0 mm 0 mm + 50 mm
	Experiment	69.6	207	5	10 ± 0.2
	Prediction	71.1	204	6.2	6.0

Appendix A. Evolution of internal pressure, liquid and vapor amount At the initial state at T = 145 °C the blowing agent is mainly in liquid state. When the liquid phase exists, the internal pressure continues to increase in function of the temperature. Following Fig. 5, the TEM starts to grow at 172 °C, then the liquid is fastly evaporated. As shown in Fig. A.1 the liquid at 175 °C is fully evaporated, then the internal pressure decreases according the growth of the TEM until reaching an equilibrium. This value tends to the external pressure P ∞ which is equal to the atmospheric pressure, as the shell is too thin and too elastic to counteract the expansion.

Appendix B. Calculations of TEM and syntactic foam density

This appendix gives the equations which yield to an expression of the syntactic foam specific mass ρ s f dependent to the TEM expansion and entrance parameters. To begin with, the specific mass ρ s f can be expressed as

where ρ pm and ρ T EM are the specific mass of the thermoplastic and the TEM, respectively. X v,T PV and X v,T EM represent the volume fraction of the thermoplastic and the TEM, respectively. They are related to each others by the following conservation equation X v,pm + X v,T EM = 1. Therefore, ρ T EM can be expressed by

with ρ sh , ρ hc are the mass density of the microcapsule shell and the hydrocarbons, respectively. ρ hc can be expressed as follows

where X m,hc and X v,hc are the mass and the volume fraction of the hydrocarbons, respectively. X m,hc and X v,hc are the mass and the volume fraction of hydrocarbons, respectively. Measurements gave X m,hc = 14,1 % and X v,hc = V int /V ext . Taking initial values of V int and V ext , gives X v,hc at t = 0, denoted by X v,hc | t=0 , which leads to ρ hc | t=0 . Combining Eqs. (B.2) and (B.3) allows us to estimate the initial mass density of TEM

Hence, the initial volume fraction of TEM X v,T EM | t=0 is equal to