S Amir Bahrani

Souria Hamidouche

Masoud Moazzen

Khady Seck

Caroline Duc

Metin Muradoglu

James B Grotberg

Francesco Romanò
email: francesco.romano@ensam.eu

S Amir Bahrani
email: amir.bahrani@imt-nord-europe.fr

Propagation and rupture of elastoviscoplastic liquid plugs in airway reopening model

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

Liquid plugs in conduits and pipes are typical flow regimes encountered in two-phase microfluidics [START_REF] Tsaoulidis | Flow patterns and pressure drop of ionic liquid-water two-phase flows in microchannels[END_REF]. At small scales inertial and gravitational effects play a minor role in the dynamics of the flow, while surface tension and Marangoni stresses gain an increasing importance from millimetric down to nanometric characteristic lengths. Understanding, predicting and controlling liquid plugs in confined flows is of crucial importance for designing micron-size and labs-on-a-chip devices. Moreover, the presence of liquid plugs is commonly reported in pulmonary flow due to disease or closure of distal airways and the propagation of liquid plugs plays a major role in medical treatments based on direct instillation of liquid solutions [START_REF] Kazemi | Surfactant delivery in rat lungs: Comparing 3D geometrical simulation model with experimental instillation[END_REF]. The fundamental understanding of liquid plug dynamics in pulmonary flows represents the main motivation of our study, with a special focus on the effect of non-Newtonian rheological properties of the liquid phase.

The thin liquid film lining the inner airway walls is formed by the mucus layer and the serous layer, with the latter in contact with the airway walls [START_REF] Widdicombe | Regulation of depth and composition of airway surface liquid[END_REF]. The dynamics of such a thin film is of major importance for the formation of liquid plugs in lungs. When the airway themselves, lung injury and disease. The protective impact of surfactant, which reduces the mean surface tension, has been experimentally and numerically demonstrated by Cassidy [START_REF] Cassidy | Liquid Film Dynamics in the Pulmonary Airways[END_REF] and Halpern et al. [START_REF] Halpern | Liquid and surfactant delivery into pulmonary[END_REF], respectively. All these investigations focused on the capillary instability of the liquid film coating the airways. However, simultaneously to the Plateau-Rayleigh instability, also the elastic instability of the deformable airway wall plays an important role in the airway closure (see e.g. Halpern and Grotberg [START_REF] Halpern | Fluid-elastic instabilities of liquid-lined flexible tubes[END_REF]). A detailed study of the interplay between interfacial and elastic instabilities is reported by Heil et al. [START_REF] Heil | The mechanics of airway closure[END_REF], who showed that both such instabilities play in favor of the occlusion of the respiratory airways. Another element of complexity of the liquid plug formation in human lungs is due to the non-Newtonian behavior of mucus, which includes viscoelasticity [START_REF] Yeates | The Lung: Scientific Fundations -Mucus Rheology[END_REF][START_REF] Romanò | The effect of viscoelasticity in an airway closure model[END_REF], yield stress and sheardependent viscosity (shear-thinning [START_REF] Quraishi | The rheology of nasal mucus: A review[END_REF]). The non-Newtonian effects have a significant impact on airway closure as shown by Halpern et al. [START_REF] Halpern | The effect of viscoelasticity on the stability of a pulmonary airway liquid layer[END_REF] and their characterization represents an active field of research in pulmonary fluid mechanics.

Once a liquid plug occludes the airway, the air-blown plug propagates inside the lungs, either rupturing or persisting and splitting when it encounters an airway bifurcation. The former scenario occurs when the precursor liquid film is thinner than the trailing film [START_REF] Halpern | Nonlinear saturation of the Rayleigh instabilitydue to oscillatory flow in a liquid-lined tube[END_REF][START_REF] Fujioka | Steady propagation of a liquid plug in a two-dimensional channel[END_REF][START_REF] Grotberg | Respiratory fluid mechanics[END_REF] . This phenomenon is called airway reopening, and it causes a pressure wave detectable by a stethoscope [START_REF] Piirila | Crackles: Recording, analysis and clinical significance[END_REF]. The reopening of respiratory airways has been extensively investigated because of its damaging effects due to very high stress levels at the airway walls. The capillary wave characteristic of the retracting air finger (see Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF], Taylor [START_REF] Taylor | Deposition of a viscous fluid on the wall of a tube[END_REF]) and the dynamics immediately after rupture [START_REF] Huh | Acoustically detectable cellular-level lung injury induced by fluid mechanical stresses in microfluidic airway systems[END_REF] are responsible of such a high level of wall stresses and of the consequent lethal or sub-lethal response of epithelial cells [START_REF] Huh | Acoustically detectable cellular-level lung injury induced by fluid mechanical stresses in microfluidic airway systems[END_REF][START_REF] Bilek | Mechanisms of surface-tension-induced epithelial cell damage in a model of pulmonary airway reopening[END_REF][START_REF] Kay | Pressure gradient, not exposure duration, determines the extent of epithelial cell damage in a model of pulmonary airway reopening[END_REF]. The strong correlation between plug rupture and high wall stress levels has been confirmed by the simulations of Fujioka and Grotberg [START_REF] Fujioka | Steady propagation of a liquid plug in a two-dimensional channel[END_REF][START_REF] Fujioka | The steady propagation of a surfactant-laden liquid plug in a two-dimensional channel[END_REF], Fujioka et al. [START_REF] Fujioka | Unsteady propagation of a liquid plug in a liquid-lined straight tube[END_REF], Hassan et al. [START_REF] Hassan | Adaptive Lagrangian-Eulerian computation of propagation and rupture of a liquid plug in a tube[END_REF], and Muradoglu et al. [START_REF] Muradoglu | Effects of surfactant on propagation and rupture of a liquid plug in a tube[END_REF], who considered an airway modeled either as a rigid pipe or a rigid channel. Further confirmations are reported by the in-vivo experiments by Muscedere et al. [START_REF] Muscedere | Tidal ventilation at low airway pressures can augment lung injury[END_REF] and Taskar et al. [START_REF] Taskar | Surfactant dysfunction makes lungs vulnerable to repetitive collapse and reexpansion[END_REF]. Other investigations focused on the effect of surfactant that plays a crucial role in reducing the surface tension at the liquid-air interface [START_REF] Halpern | Surfactant effects on fluid-elastic instabilities of liquidlined flexible tubes: A model of airway closure[END_REF]. Fujioka and Grotberg [START_REF] Fujioka | The steady propagation of a surfactant-laden liquid plug in a two-dimensional channel[END_REF] numerically investigated the effect of surfactant on the steady propagation of a liquid plug within a two-dimensional channel. The results showed that the mechanical stresses induced by the plug propagation is reduced even if the preexisting surfactant concentration is small. The same observation has been also made by Tavana et al. [START_REF] Tavana | Epithelium damage and protection during reopening of occluded airways in a physiologic microfluidic pulmonary airway model[END_REF], who performed an experimental study on the effects of surfactant on cell injury. Adding a physiologic amount of Survanta, they showed that plug propagation significantly reduces the lethal response of epithelial cells. This has been confirmed by the extensive simulations of Muradoglu et al. [START_REF] Muradoglu | Effects of surfactant on propagation and rupture of a liquid plug in a tube[END_REF], capable to reproduce the whole rupture process, including the post-rupture phase. The effect of compliance of airway walls was considered by Zheng et al. [START_REF] Zheng | Liquid plug propagation in flexible microchannels: A small airway model[END_REF], who showed an enhancement of the stress level on the airway walls due to wall deformability. The effect of yield stress on plug rupture has been experimentally and numerically investigated by Hu et al. [START_REF] Hu | A microfluidic model to study fluid dynamics of mucus plug rupture in small lung airways[END_REF][START_REF] Hu | Effects of surface tension and yield stress on mucus plug rupture: A numerical study[END_REF], Low et al. [START_REF] Low | Pulmonary airway reopening: Effects of non-Newtonian fluid viscosity[END_REF] and Zamankhan et al. [START_REF] Zamankhan | Steady motion of bingham liquid plugs in two-dimensional channels[END_REF]. They either dealt with liquid plugs on a dry or a pre-wetted channel, focusing on the dynamics of the non-Newtonian liquid plug up to the rupture event. Other numerical and experimental studies focused on a train of plugs and discussed catastrophic events in the airway network [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF] or surfactant-replacement therapies where the liquid plugs are formed due to the liquid instilled into airways for medical treatment [START_REF] Kazemi | Surfactant delivery in rat lungs: Comparing 3D geometrical simulation model with experimental instillation[END_REF][START_REF] Filoche | Three-dimensional model of surfactantreplacement therapy[END_REF]. For literature reviews of respiratory airway closure, liquid plug propagation, and rupture, we refer to Grotberg [START_REF] Grotberg | Respiratory fluid mechanics[END_REF][START_REF] Grotberg | Pulmonary flow and transport phenomena[END_REF][START_REF] Grotberg | Respiratory fluid mechanics and transport processes[END_REF].

In this paper we focus on the impact of the non-Newtonian properties of mucus on the liquid plug rupture. We perform a detailed rheological analysis of various samples of synthetic mucus and employ such samples to study airway reopening. Finally, we derive a reduced-order model to shed some light on the standalone effect of viscoplasticity and viscoelasticity. The paper is structured as follows:

In Section 2 the experimental set-up and the corresponding results are described; In Section 3 we define the framework of our reducedorder model, derive such a model based on up-to-date correlations reported in literature, and presents the corresponding results; finally, the conclusions of our study are drawn in Section 4.

Experiments

Synthetic mucus samples

The synthetic mucus used in this study is a compound of water and two polymers. The Scleroglucan Polysaccharide (Actigum ® ) is employed as a chemical agent to increase the liquid viscosity and Locust Bean Gum (Vidogum ® L-175) is used as a thickening agent. We dissolve 0.25 g of Vidogum in 50 ml of distilled water for 2 h and then different Actigum concentrations ( cact = 0.25, 0.5, 1 and 1.5 g/ml.) are added to such aqueous solutions. Copper sulfate CuSO 4 (0.05 g) is also added as a stabilizer before keeping the mixtures under agitation conditions for 24 h at room temperature (around 21 • C). After that, the aqueous solutions are left at rest for 24 h in order to recover the gel structure before performing any other manipulation. The solutions are stocked in close containers that prevent significant evaporation as they are stored at 4 • C. In the followings, we name them after their characteristic Actigum concentrations, i.e., MUC-0.25%, MUC-0.5%, MUC-1% and MUC-1.5% correspond to cact = 0.25, 0.5, 1 and 1.5 g/ml, respectively.

As aforementioned, we use vidogum and CuSO 4 for preparing our synthetic mucus. This is different from the experimental protocol proposed by Lafforgue et al. [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF], who rather employ viscogum and NaCl. We stress that our choice of CuSO 4 instead of NaCl is intended to avoid the spontaneous development of microbiological culture through the synthetic mucus. As a result, the lifetime of the final aqueous solutions is increased.

Rheological measurements

The rheological behavior of the synthetic mucus is characterized using a cone/plate geometry (50 mm/1 • ) with truncation gap of 𝜃 = 0.104 mm. All measurements are carried out in an Anton-Paar MCR 302 rheometer at a constant temperature of 21 • C. Two flow sweep protocols are applied to measure the rheological characteristics of the synthetic mucus and the relaxation time. At first, in order to find viscosity and fitting Herschel-Bulkley model the flow curve measurement was done. This was happened by following a standard protocol of flow sweep [START_REF] Ahmadikhamsi | Impact of surfactant addition on non-Newtonian fluid behavior during viscous fingering in Hele-Shaw cell[END_REF][START_REF] Bahrani | Intermittency in the transition to turbulence for a shear-thinning fluid in hagen-poiseuille flow[END_REF] within the shear-rate interval of 0.01/s < ̇𝛾 < 1000/s. Then to get the dynamic, shear rate dependent, relaxation time 𝜆, we followed step by step the protocol detailed by Casanellas et al. [START_REF] Casanellas | The stabilizing effect of shear thinning on the onset of purely elastic instabilities in serpentine microflows[END_REF]. This experimental protocol is consisted of several incremental steps. In each step a constant shear rate is imposed with multiple measuring points. Right-after a flow sweep test is conducted where the solution was sheared at a small shear rate of reference ( ̇𝛾 = 1𝑠 -1 ). On each step, subtracting the means of two values of the first normal stress difference 𝑁 1 gives a way to resolve the instrumental drift of the normal force and find real value of 𝑁 1 ( ̇𝛾). This latter will be completed by removing the contribution of fluid inertia to the normal force. Moreover, it can clarify the non-linear rheological behavior of the different samples in which the real rheological response of the samples is shear-rate dependent.

A shear-thinning behavior is observed for all the aqueous solutions. The resulting rheological curves are characterized by fitting the Herschel-Bulkley model, 𝜏 = 𝜏 0 + 𝐾 ̇𝛾𝑛 , where 𝜏, 𝜏 0 , 𝐾, and 𝑛 are the shear stress, yield stress, consistency index, and power-law index, respectively (see Fig. 1). The rheological parameters of the synthetic mucus solutions are summarized in Table 1. The results show that the yield stress for MUC-0.25% is 0.59 Pa and it increases upon an increase of the Actigum concentration ( c𝑎𝑐𝑡 ). A similar trend is also observed for the consistency index, the power-law index, 𝑛, slightly decreases as cact 1.

Table 1

Rheological parameters of the synthetic mucus suspensions derived by fitting the Herschel-Bulkley model 𝜏 = 𝜏 0 + 𝐾 ̇𝛾𝑛 as well as the White-Metzner model [START_REF] White | Development of constitutive equations for polymeric melts and solutions[END_REF] in order to find the relaxation time, 𝜆. All solutions follow the relations of 𝑁 1 ∝ ̇𝛾𝛼 and 𝜆 ∝ ̇𝛾𝛽 , where the exponents 𝛼 and 𝛽 govern the viscoelastic behavior of the solutions.

Sample

Herschel increases. This shows that the solutions become more shear-thinning as cact increases.

The viscoelastic behavior of the synthetic mucus can be described by employing various rheological models. One way is to assume that the relaxation time is independent of the shear rate. This approach is typically employed in the Oldroyd-B, the FENE-P and FENE-CR models, applied generally to solutions with low concentration of polymers. However, in the case of non-linear rheological behavior, likely for higher polymer concentrations, we can use the White-Metzner constitutive equation [START_REF] White | Development of constitutive equations for polymeric melts and solutions[END_REF]. This model (see Eq. ( 1)) assumes that the relaxation time, shear viscosity and first normal-stress difference 𝑁 1 are shear-dependent

𝑁 1 = 2 [ 𝜂( ̇𝛾) -𝜂 s ] 𝜆( ̇𝛾) ̇𝛾2 (1)
where are obtained experimentally using the rheometer. However, as already mentioned, the raw data of 𝑁 1 was corrected by subtracting the inertial contribution. This was done by subtracting the negative inertial effect which tends to pull down the cone due to the presence of liquid between the cone-plate rheometer instrument [START_REF] Macosko | Rheology: Principles, Measurements, and Applications[END_REF]. The corresponding correction is given by: 𝑁 1 = 𝑁 1,𝑡𝑜𝑡 + 0.15𝜌𝛺 2 𝑅 2 , where 𝑁 1,𝑡𝑜𝑡 is the total/row normal-stress difference directly given by the rheometer and 𝛺 = tan (𝜃) ̇𝛾 is the angular velocity. It is worth to mention that these corrections rely on the accurate measurements of 𝑁 1 . This latter depends on resolution of the normal force sensor of the rheometer, especially at the low polymer concentrations and shear rates. As can be seen in Fig. 2(a) the evolution of 𝑁 1 as a function of shear rate shows a power law like 𝑁 1 ∝ ̇𝛾𝛼 . The results show that the slope of 𝛼 decrease by increasing of shear-thinning degree depending on the concentration of Actigum cact (see Table 1). Indeed the value of 𝛼 determines the complexity of non-linear rheological response of the fluids, for 𝛼 = 2, the behavior follows the Oldroyd-B model. The relaxation time as a function of shear rate (1)b is also represented by a power law function as 𝜆 ∝ ̇𝛾𝛽 . The values of 𝛽, summarized in Table 1, indicate larger relaxation times as the concentration of Actigum increases. Fig. 2 also demonstrates that as the shear-thinning degree increases the relaxation time becomes more sensible to shear rate.

Experimental apparatus

The experiments of a liquid plug propagation in a pipe flow were conducted using a cylindrical glass tube with inner diameter of 𝑑 𝑖 = 2𝑎 = 4 mm corresponding to the diameter of bronchioles at 6th-to -7th generation of an adult human lung. In order to reproduce the dynamics of a liquid plug, a controlled volume of the synthetic mucus is injected, at first, through the glass tube, and then pushed along using two different gases as described in detail below.

The mucus and gas injections are performed at a constant flow rate by using a flowmeter (Flow Unit, Fluigent), as shown in Fig. 3(a). A systematic pre-wetting surface treatment of the glass tube is done by cleaning with acetone, isopropanol and dichloromethane before each experiment. The impact of the rheological properties of synthetic mucus is studied under atmospheric conditions, using filtered air directly coming from the room.

The dynamic rupture of liquid plugs was recorded using a CCD camera (DMK 23U618, Monochrome, up to 120 images/s, The Imaging Source, Bremen, Germany) connected to a PC. The color difference between the liquid plug (white) and the background (dark gray) is exploited in our in-house Matlab ® code to post-process the dynamics of the liquid plug (see Fig. 3(b)). Each snapshot is processed by extracting a stripe of the width of one pixel along the projection of the pipe axis. The corresponding space-time diagram is thereafter obtained by gathering all the one-pixel stripes and sorting them in chronological order. As the background of the experimental setup is dark-gray and the liquid plug is white, the dynamics of the liquid plug is readily identified as illustrated in Fig. 3(b). A final image polishing is performed by using an in-house front-identification algorithm that takes advantage of the sharpness of the space-time diagram, and it characterizes the brightness contour between the white pixels of the plug and the darkgray pixels of the background. Within this final post-processing step, we also remove the blind region of the space-time diagram (see vertical gray stripe in Fig. 3(b) due to a metallic plate employed in our experiment to hold the capillary pipe in place). The image-processing protocol readily provides the parameters of interest in our statistics, i.e. the initial plug length, 𝐿 i = 𝐿 p (𝑡 = 𝑡 0 ), the final plug length, 𝐿 f = 𝐿 p (𝑡 = 𝑡 r ), the rupture time, 𝑡 r , and the distance covered by the plug before the rupture event takes place (rupture length 𝐿 r = 𝑧(𝑡 = 𝑡 r )).

A schematic of such quantities is presented in Fig. 3(b). Finally, the average plug velocity is estimated by the ratio of the rupture length to the rupture time, i.e. 𝑈 p = 𝐿 r ∕𝑡 r .

Dynamics of plug ruptures

The dynamics of the liquid plug in the capillary tube is characterized for several flow conditions using the four samples of synthetic mucus and varying the initial plug length 𝐿 i . Fig. 4 depicts the space-time diagram of the liquid plug rupture for two different concentrations, i.e. MUC-0.5% and MUC-1.5%, and three different initial plug lengths, 𝐿 i = 5, 13 and 22 mm. The markers denote the data resulting from our image post-processing protocol, and the dashed lines are cubic polynomials obtained by least-squares fit for tracking the rear and the front menisci. Cubic polynomials are employed because they result in a smooth and seamless interpolation while avoiding oscillations known as Runge phenomena [START_REF] Vetterling | Numerical Recipes in Fortran[END_REF]. Each subfigure includes two insets representing two experimental snapshots where the plug moves from left to right. On the top-right corner we depict the liquid plug as soon as both front and rear menisci entered the observation window of our experimental setup. This condition defines 𝐿 i = 𝐿 p (𝑡 = 0). On the bottom-left corner we depict the liquid plug right before rupture and we identify the final plug length 𝐿 f = 𝐿 p (𝑡 = 𝑡 r ). Hence, the time taken to pass from 𝐿 p = 𝐿 i to 𝐿 p = 𝐿 f is the rupture time 𝑡 r . The white stripe at about two thirds from the pipe inlet (see Fig. 3(b)) is a fixed metal plate employed to hold the capillary tube in place. Due to it, the plug cannot be detected when passing below such a metal plate (gray stripe in Fig. 4).

The four panels of Fig. 4 report four representative experimental space-time diagrams. Fig. 4(a) shows a typical space-time diagram observed for short plugs (𝐿 i = (𝑎)). Regardless of the Actigum concentration, the short plugs of our experiments always rupture following an axi-symmetric dynamics. The progressive drainage of the plug up to 𝐿 p (𝑡 r ) = 𝐿 r ≈ 0 is also observed for low Actigum concentrations, i.e. MUC-0.25% and MUC-0.5%, regardless of the initial plug length. This is demonstrated in Fig. 4(b). Upon an increase of the Actigum concentration, the liquid plug becomes more prone to asymmetric ruptures (see bottom-left inset of Fig. 4(c)). A detailed dynamics of the asymmetric plug rupture for MUC-1.0% and 𝐿 i = 13 mm is depicted in Fig. 5. By increasing the viscoplastic character of the liquid plug, the deposition of the liquid film becomes wavy and loses the axi-symmetry observed for low Actigum concentrations. This same phenomenon has been reported in figs. 14 to 16 by Caliman et al. [START_REF] Caliman | An experimental investigation on the Newtonian-Newtonian and viscoplastic-Newtonian displacement in a capillary tube[END_REF] and it results in a non-zero centerline plug length at rupture (see also space-time diagram of Fig. 4(c)). Another scenario is finally depicted in Fig. 4(d): A long liquid plug undergoes an initial drainage and then gets steadily advected in the capillary pipe without rupturing. Fig. 6 depicts the experimental findings for our four samples of synthetic mucus upon a change of the average plug speed 𝑈 𝑝 , the effective viscosity 𝜇, the relaxation time 𝜆, and the yield stress 𝜏 0 . Two quantities are employed for characterizing the dynamics of the liquid plug: (i) the rupture length non-dimensionalized by the difference between initial and final length 𝐿 𝑟 ∕(𝐿 𝑖 -𝐿 𝑓 ) and (ii) the rupture time. The first nondimensional group is employed to compare the distance traveled by a liquid plug with the total plug-length reduction until the rupture. Fig. 6(a,b) show increasing the average plug length does not correlate with a change of non-dimensional plug length, but it rather impacts on the rupture time. This is a combined effect of (i) the faster propagating plug that experiences higher capillary numbers, hence thicker coatings, and (ii) the higher-strained plug with therefore a lower effective viscosity (shear thinning) and lower relaxation times (see Figs. 1 and2). Also the yield stress plays a role in delaying the rupture, even if the component of the delay due to higher yield stresses seem to be a higherorder correction when compared to the capillary-number effect. The change of the plug rheology is rather more impactful on the normalized rupture length, that grows upon an increase of 𝜇, 𝜆 and 𝜏 0 . This is consistent with the physiological observations of non-Newtonian mucus plugs that develop a higher resistance to rupture when the polymeric concentration is increased. We further stress that our trend for 𝑡 𝑟 is in agreement with the results of Hu et al. [START_REF] Hu | A microfluidic model to study fluid dynamics of mucus plug rupture in small lung airways[END_REF], who showed that the liquid Fig. 7. Schematic of the reduced-order model for the liquid plug rupture in an airway of constant radius 𝑎. The rear air finger (dark blue) has a film thickness ℎ r and moves with tip velocity 𝑈 r , while the front air finger (dark red) has a film thickness ℎ f and recedes with tip velocity 𝑈 f . The core fluid (dotted box) is described by the recirculating flow 𝑢 p superposed to the averaged velocity 𝑈 p = (𝑈 r + 𝑈 f )∕2.

plug takes a longer time to deform and rupture at higher Bingham numbers. Therefore, they attributed the longer rupture time to the slowed plug deformation induced by high yield stresses. Moreover, a larger scatter of the experimental data is observed when c𝑎𝑐𝑡 ↑, which is due to a higher probability to run into an asymmetric rupture of the liquid plug.

Reduced-order model of the liquid plug rupture

A reduced-order model of our experiment is formulated as a onedimensional model based the concept flow across a liquid plug. Referring to Fig. 7, we assume that the pressure in the advancing air finger at the trailing meniscus (rear film, subscript r) of the liquid plug, 𝑃 r , and the pressure in the receding air finger at the leading meniscus (front film, subscript f) of the liquid plug, 𝑃 f , are constant. The liquid plug moves therefore under the effect of a given pressure difference across the plug, 𝛥𝑃 = 𝑃 r -𝑃 f . This imposed pressure difference is opposed by the pressure drop 𝛥𝑃 p due to the capillary and viscous forces. In the absence of inertia, the pressure drop is proportional to the liquid flow rate 𝑄 and can be written as 𝛥𝑃 p = 𝑅 p 𝑄, where 𝑅 𝑝 is the liquid plug resistance. Assuming that the front and the rear menisci do not interact, the resistance of the plug can be expressed as a system of three resistances in series, i.e. 𝑅 p = 𝑅 r +𝑅 c +𝑅 f , where 𝑅 r and 𝑅 f are the resistances due to the interaction of surface tension with the viscous rear and front thin film, respectively, and 𝑅 c is the viscous resistance of the core. Similar modeling approaches have been proposed in literature by Howell et al. [START_REF] Howell | The propagation of a liquid bolus along a liquid-lined flexible tube[END_REF] and Waters and Grotberg [START_REF] Waters | The propagation of a surfactant laden liquid plug in a capillary tube[END_REF], who treated the plug core inviscid, and relegated the resistance to the thin films ahead and behind the liquid plug.

Following the model of Fujioka et al. [START_REF] Fujioka | Reduced-dimension model of liquid plug propagation in tubes[END_REF] for a clean Newtonian liquid plug, the following system of ODEs can be used to represent the plug propagation.

𝑑 𝑑𝑡 ⎡ ⎢ ⎢ ⎣ 𝑍 p 𝑈 p 𝐿 p ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ ( 𝑈 r + 𝑈 f ) ∕2 ( 𝛥𝑃 -𝑅 p 𝑈 p 𝜋𝑎 2 ) ∕𝜌𝐿 p -(𝑈 p ∕𝐿 p )𝑑 𝑡 𝐿 p 𝜀 f ( 2 -𝜀 f ) 𝑈 f -𝜀 r ( 2 -𝜀 r ) 𝑈 r ⎤ ⎥ ⎥ ⎦ , ( 3
)
where 𝜌 is the density of the liquid, 𝑡 is the time, 𝑈 p , 𝑈 r and 𝑈 f are the plug, rear-bubble and front-bubble velocities, respectively, 𝐿 p is the plug length, and 𝜀 r = ℎ r ∕𝑎 and 𝜀 f = ℎ f ∕𝑎 are the non-dimensional rear and front liquid film thicknesses, respectively (see Fig. 7). The first equation describes the location of the center of the liquid plug, 𝑍 p , with 𝑑𝑍 p ∕𝑑𝑡 = 𝑈 p , the second equation is derived from the onedimensional Navier-Stokes equation, and the third one is the continuity equation. A detailed characterization of the Newtonian resistances is reported in the following section, partially modifying the model of Fujioka et al. [START_REF] Fujioka | Reduced-dimension model of liquid plug propagation in tubes[END_REF] to take into account the recirculation resistance in the liquid plug (see 𝑢 p in Fig. 7). The derivation of the last two equations is briefly summarized in Appendix. We further generalize this model by including non-Newtonian effects. Such extension of the model of Fujioka et al. represents an original contribution of our study.

Newtonian plug resistance

A Newtonian liquid plug travels in a capillary tube of radius 𝑎 and is characterized by trailing and a leading film thicknesses ℎ r and ℎ f , respectively, a length 𝐿 p and a mean traveling velocity 𝑈 p . Under the hypothesis that the rear and the front menisci do not interact, the rear meniscus can be characterized by the Bretherton lubrication model of a steady Taylor bubble advancing with velocity 𝑈 r . The pioneering work of Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] is concerned with asymptotically small capillary numbers of the rear meniscus, i.e. 𝐶𝑎 r = 𝜇𝑈 r ∕𝜎 → 0, where 𝜎 is the surface tension at the liquid-gas interface and 𝜇 is the dynamic viscosity of the liquid phase. The results of Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] have successively been extended by several other authors that corrected the Bretherton's correlation for finite 𝐶𝑎 r . The most refined correction has recently been proposed by Fujioka et al. [START_REF] Fujioka | Reduced-dimension model of liquid plug propagation in tubes[END_REF], who confirmed the formula given by Aussillous and Quéré [START_REF] Aussillous | Quick deposition of a fluid on the wall of a tube[END_REF] for the steady film thickness ℎ r and improved the result for the resistance 𝑅 r by fitting Stokesian numerical simulations of an advancing, steady, semi-infinite Bretherton bubble. The formulae by Fujioka et al. [START_REF] Fujioka | Reduced-dimension model of liquid plug propagation in tubes[END_REF] read

𝜀 r = ℎ r 𝑎 = 1.34𝐶𝑎 2∕3 r 1 + 1.34 × 2.5𝐶𝑎 2∕3 r , 𝑅 r = 𝜇 𝜋𝑎 3 ( 1 -𝜀 r ) 2 ( 2 × 1.79 × 2 2∕3 𝐶𝑎 -1∕3 r 1 + 1.41𝐶𝑎 1∕3 r + 1.1
) .

(

) 4
The film resistance of the front meniscus is considered by the correction proposed by Fujioka et al. [START_REF] Fujioka | Reduced-dimension model of liquid plug propagation in tubes[END_REF] that takes into account a finite capillary number and does not require ℎ f ≪ 𝑎. The formula derived by fitting the Stokes flow for a rear front of a semi-infinite Bretherton bubble is

𝑅 f = 𝜇 𝜋𝑎 3 ( 1 -𝜀 f ) 2 ( -2 × 3 2∕3 𝑓 𝐶𝑎 -1∕3 f ) , 𝑓 = 5 ∑ 𝑘=0 𝐶 𝑘 ( log 10 𝛼 ) 𝑘 + 1.02𝜀 0.348 f 𝛼 -0.594 1 + 0.0778𝛼 -0.594 (5)
where The last component of the liquid plug resistance is 𝑅 c that is obtained following the approach of Che et al. [START_REF] Che | An analytical model for plug flow in microcapillaries with circular cross section[END_REF], who considered the recirculating axisymmetric flow ⃗ 𝑢 𝑝 in the reference frame co-moving with the plug at speed 𝑈 p . Since a change of reference frame does not affect the flow resistance, we employ the formula reported in Che et al. [START_REF] Che | An analytical model for plug flow in microcapillaries with circular cross section[END_REF]

𝐶𝑎 f = 𝜇𝑈 f ∕𝜎, 𝑈 f is the velocity of the front bubble, 𝜀 f = ℎ f ∕𝑎, 𝛼 = 𝜀 f ∕ ( 3𝐶𝑎 f ) 2∕3 , 𝐶 0 = -0.149, 𝐶 1 = 2.
𝑅 c = - 2𝜇𝐿 p 𝜋 2 𝑎 4 ∞ ∑ 𝑛=1 1 -(-1) 𝑛 𝑛 𝛼 𝑛 𝛽 𝑛 , 𝛼 𝑛 = 2[(-1) 𝑛 -1] 𝑛𝜋𝜙 𝑛 , 𝛽 𝑛 = - 2𝑛𝜋𝑎𝐼 1 ( 𝑛 2 𝜋 2 𝑎 2 𝐿 2 p ) 𝐿 p 𝐼 2 ( 𝑛 2 𝜋 2 𝑎 2 𝐿 2 p )
𝜙 𝑛 = 𝑛𝜋𝑎𝐼 0 ( 𝑛𝜋𝑎 𝐿 p ) 𝐼 2 ( 𝑛𝜋𝑎 𝐿 p ) -𝑛𝜋𝑎𝐼 2 1 ( 𝑛𝜋𝑎 𝐿 p ) 𝐿 p 𝐼 2 1 ( 𝑛𝜋𝑎 𝐿 p ) (6)
where 𝐼 0 , 𝐼 1 and 𝐼 2 are the modified Bessel function of the first kind of order 0, 1 and 2, respectively. The core resistance reported in (6) assumes straight menisci, dry plug transport and tends to the viscous resistance of a fully-developed Hagen-Poiseuille profile for long liquid plugs, i.e. 𝐿 p ≫ 𝑎. Owing to the slow convergence of ( 6), it is computationally expensive to evaluate infinite hence we prefer to employ a least-square fit built on the following power law

𝑅 c = 𝜇𝐿 p 2𝜋𝑎 4
[ 16 + 77.9725 ( 𝐿 p 𝑎 + 0.5546

) -1.7670 ] , ( 7
)
where the Hagen-Poiseuille asymptote enters enforcing the first coefficient of the fit, i.e. 16. A comparison between ( 6) and ( 7) is reported in Fig. 8.

Corrections for non-Newtonian effects

The effect of the non-Newtonian behavior of the liquid phase modifies all three flow resistances of our reduced-order model. To the best of the authors' knowledge, no empirical, numerical or theoretical correlations are present in the literature for taking into account the viscous, elastic and plastic behavior of the mucus simultaneously. In our model, therefore, we assume that the viscoelastic and viscoplastic corrections can be superposed.

Following Gauri and Koelling [START_REF] Gauri | The motion of long bubbles through viscoelastic fluids in capillary tubes[END_REF], the viscoelasticity of a liquid strongly affects the film thickness. Based on their experimental measurements, they propose the following correction to pass from the Newtonian to the viscoelastic (VE) case

𝜀 * VE (2 -𝜀 * VE ) = 𝜀 * (2 -𝜀 * )[0.1217 ln(𝐷𝑒 * ) + 0.8461] (8)
where 𝐷𝑒 * = 2𝑈 * 𝜆∕𝑎 is the Deborah number and 𝜆 is the relaxation time of the liquid. This correlation is valid for 𝐷𝑒 * ≥ 3, and higherorder corrections required for 𝐷𝑒 * < 3 are neglected in our model. No viscoelastic correction is taken into account for the core resistance. This is consistent with the limit of the Hagen-Poiseuille profile for long plugs.

Several studies have focused on the correction of the rear-and frontbubble resistances and the corresponding non-dimensional thicknesses for Herschel-Bulkley liquids (see e.g. Laborie et al. [START_REF] Laborie | Yield-stress fluid deposition in circular channels[END_REF] and Jalaal and Balmforth [START_REF] Jalaal | Long bubbles in tubes filled with viscoplastic fluid[END_REF]). However, a comprehensive explicit correlation for 𝑅 𝑝 and 𝜀 * between Newtonian and Herschel-Bulkley liquids (hence as a function of 𝐾, 𝜏 0 , and 𝑛) is still missing. Hence, we limit the viscoplastic correction of our model to Bingham fluids (𝑛 = 1) and correct the film thickness and the plug resistance by fitting the results of Zamankhan et al. [START_REF] Zamankhan | Steady motion of bingham liquid plugs in two-dimensional channels[END_REF][START_REF] Zamankhan | Steady displacement of long gas bubbles in channels and tubes filled by a bingham fluid[END_REF]. In this case 𝐶𝑎 * reduces to the classic definition 𝐶𝑎 * = 𝜇𝑈 * ∕𝜎. The resulting corrections are

𝛥𝑃 VP p = 𝛥𝑃 p + 2.3𝐵𝑖 p 𝜇𝑈 p 𝑎 , ( 9a
)
𝜀 𝑉 𝑃 * (2 -𝜀 𝑉 𝑃 * ) 𝜀 * (2 -𝜀 * ) = 1 + 𝐵𝑖 * 𝐶𝑎 * { 0.06 [ tanh ( 20 -𝐵𝑖 * ∕𝐶𝑎 * 𝜋 ) + 1 ] + 0.03 [ tanh ( 𝐵𝑖 * ∕𝐶𝑎 * -20 𝜋 ) + 1 ]} ( 9b
)
where 𝐵𝑖 is the Bingham number 𝐵𝑖 * = 𝜏 0 𝑎∕𝜇𝑈 * , 𝜏 0 is the yield stress, and the fits are reported in Fig. 9. In our fitting model, the control parameter 𝐵𝑖∕𝐶𝑎 plays a relevant role in describing the viscoplastic correction of the thin film thickness. To better understand it, we consider the dynamics of a thin liquid film, i.e. 𝜀 = ℎ∕𝑎 ≪ 1. Within the framework of the thin-film approximation, the pressure inside the liquid film is simply projected in direction normal to the wall. Hence, to leading order, the pressure inside the liquid (𝑝 relative to the air core pressure) is 𝑝 𝜎 = 𝜎∕(𝑎 -ℎ) ≈ 𝜎∕𝑎. For slow viscoplastic liquid plugs, and especially in the thin film coating the capillary pipe, inertial effects are negligible and the film thickness results out of the balance between viscoplasticity (hence yield stress 𝜏 0 ), capillary pressure (𝑝 𝜎 ) and viscous stresses (𝜏 𝑣 = 𝜇𝑈 * ∕𝑎). Hence, the ratio between the Bingham number 𝐵𝑖 = 𝜏 0 𝑎∕𝜇𝑈 * and the capillary number 𝐶𝑎 = 𝜇𝑈 * ∕𝜎 represents a dimensionless group that quantifies the relative importance of yield and capillary-induced normal stresses inside the liquid film against the viscous stresses, i.e. 𝐵𝑖∕𝐶𝑎 = 𝜏 0 𝑝 𝜎 ∕𝜏 2 𝑣 . In this sense, the non-dimensional group 𝐵𝑖∕𝐶𝑎 represents a parameter that involves all the stresses relevant for the dynamics of a viscoplastic thin film. The choice of 𝐵𝑖∕𝐶𝑎 as a control parameter is confirmed by the data of the independent study of Zamankhan et al. where data almost collapse on the 𝐵𝑖∕𝐶𝑎-dependent fitting curve employed by our model. Considering that three very different capillary numbers have been simulated for producing the data used in Fig. 9, the remarkable data collapse is a strong indication that 𝐵𝑖∕𝐶𝑎 is a significant control parameter. We stress that the ranges of capillary and Bingham numbers used in the original study of Zamankhan et al. and fitted in (8) cover the range of physiological 𝐶𝑎 and 𝐵𝑖 of interest for the reopening of a distal airway. In this sense, the fit is expected to be general for the purpose of its application.

Since the elastoviscoplastic corrections include a change of the non-dimensional film thicknesses, two variants M1 and M2 of the reduced-order model should be taken into account. When the thickness of the liquid films is corrected, either we consider the bubble-tip velocity constant (model M1) or we assume that the flow rate does not change (model M2). Further details are provided in the following subsection.

Algorithm of the reduced-order model

We model the liquid plug rupture of an elastoviscoplastic liquid plug assuming: The following algorithm defines our reduced-order model for an elastoviscoplastic liquid plug:

1. Enforce the front-bubble film thickness 𝜀 f = 𝜀 EVP f ;

2. Enforce the Bingham-to-capillary ratio and the Deborah number of the front-bubble, i.e. 𝐵𝑖 f ∕𝐶𝑎 f and 𝐷𝑒 f , respectively; 3. Remove the effect of the fluid elastoviscoplasticity by inverting Eqs. ( 9)(right) and ( 8) to compute the non-dimensional film thickness of a bubble retracting in a Newtonian liquid 𝜀 f ;

Note: We test two variants for the non-Newtonian film thickness corrections (superscript corr) of the reduced-order model:

M1: velocity-matched model 𝑈 corr f = 𝑈 f with 𝛥𝑃 f = -2𝜎𝑓 (3𝐶𝑎 f ) 2∕3 ∕𝑎; M2: flow-rate-matched model 𝑈 corr f = 𝑈 f [(1 -𝜀 f )∕(1 - 𝜀 surf f )] 2 with 𝛥𝑃 f = -2𝜎𝑓 (3𝐶𝑎 f ) 2∕3 ∕𝑎;
4. Compute the capillary number and the film thickness of the front bubble by inverting (4)(left), which holds for retracting bubbles, too; 5. Compute the front-bubble velocity 𝑈 f from the capillary number of a Newtonian bubble; 6. Compute the pressure drop of a Newtonian bubble at the frontbubble tip using 𝛥𝑃 f ; 7. Based on the selected 𝐵𝑖 f ∕𝐶𝑎 f and 𝐷𝑒 f , compute the corresponding 𝜏 0 and 𝜆;

Note: 𝑈 f and 𝛥𝑃 f are derived by the enforced front-bubble film thickness 𝜀 f , and by the rheological properties of mucus. Hence, none among 𝑈 f , 𝛥𝑃 f and 𝜀 f is affected by the liquid plug dynamics;

8. Compute the variation in time of the liquid plug velocity by using the following algorithm: combining the momentum and continuity equations, the time derivative of the plug velocity is expressed as

𝑑𝑈 p 𝑑𝑡 = 𝛥𝑃 -𝛥𝑃 𝑘 p 𝜌𝐿 𝑘 p - 𝑈 𝑘 p [ 𝜀 f ( 2 -𝜀 f ) 𝑈 f -𝜀 𝑘 r ( 2 -𝜀 𝑘 r ) 𝑈 𝑘 r ] 𝐿 𝑘 p ; (10) 9.
Compute the variation in time of the liquid plug length by using the continuity equation:

𝑑𝐿 p 𝑑𝑡 = 𝜀 f ( 2 -𝜀 f ) 𝑈 f -𝜀 𝑘 r ( 2 -𝜀 𝑘 r ) 𝑈 𝑘 r ( 11
)
10. Compute the time variation of the axial displacement by integrating the liquid plug velocity:

𝑑𝑍 p 𝑑𝑡 = 𝑈 𝑘 p ( 12
)
11. Integrate the three ODEs and re-initialize 𝑈 p , 𝐿 p and 𝑍 p before iterating from (viii.a) on.

The algorithm is initialized by setting 𝑈 p (𝑡 = 𝑡 0 ), 𝐿 p (𝑡 = 𝑡 0 ) and assuming 𝑍 p (𝑡 = 𝑡 0 ) = 0.

Validation against experimental measurements and numerical simulations

The experimental validation of our reduced-order model is presented for two test fluids (0.25% and 1% of Actigum concentration) and three plug driving conditions per each test fluid. Owing to the limitations of the experimental apparatus in use, our equipment does not allow for an accurate control of the driving pressure 𝛥𝑃 and for a precise measurement of the front film thickness 𝜀 f . These are three of the ten parameters of our reduced-order model, namely 𝐷𝑒, 𝐵𝑖, 𝐶𝑎, 𝜇, 𝜌, 𝐿 0 , 𝑈 𝑝 (𝑡 = 0), 𝑎, 𝛥𝑃 , and 𝜀 f .

All the parameters measured experimentally, i.e. 𝐷𝑒, 𝐵𝑖, 𝐶𝑎, 𝜇, 𝜌, 𝐿 0 , 𝑈 p (𝑡 = 0), and 𝑎, are used in the model in order to compare the modeled and the experimental centerline plug dynamics. In our experiments, 𝛥𝑃 and 𝜀 f are not measured, and they are therefore assumed as fitting parameters of our reduced-order model to match the experimental average plug speed 𝑈 p . Fig. 10. Comparison between the space-time diagrams resulting from the reduced-order model (dashed line) and the experimental measurements (light-color markers) for Act. 0.25% (left) and Act. 1% (right). Three cases are compared for each concentration: blue (Case 1), green (Case 2), and red (Case 3).

Table 2

Parameters of the comparison between the experiments and the reduced-order model. Two concentrations are considered, i.e. Act. 0.25% and Act. 1%, for three driving conditions per test liquid (Cases 1 to 3). The parameters of our comparison are summarized in Table 2 for the six experiments considered. As shown in Fig. 10, the comparison between the dynamics predicted by the model and the experimental space-time diagram is satisfactory. We stress that only a narrow neighborhood of the selected parameter tuple (𝛥𝑃 , 𝜀 f ) allows to reproduce the average plug speed of the experiments. Hence the degree of freedom of our parameter matching is very narrow, which shows that our quantitative matching is not made arbitrarily. Still, the model results are robust to small parameters variations. We further stress that, based on the results of Fig. 10 the non-linear interplay between viscoelasticity and viscoplasticity seems to play a significant role only in the final phase of the rupture, when the interface curvature (partially neglected in our model) must also play a significant role. We also expect that local interesting features are determined by the complex rheology of the fluid, but their investigation would require detailed simulations and experimental techniques which are out of the scope of our study.

A second validation of our reduced-order model is carried out by comparing the conditions of steady propagation for a viscoplastic liquid plug. The results of Zamankhan et al. [START_REF] Zamankhan | Steady motion of bingham liquid plugs in two-dimensional channels[END_REF] are assumed as reference, as they result form the numerical simulation of a two-dimensional Stokes flow for a viscoplastic liquid plug. Fig. 11 depicts the comparison between our reduced order model (markers), and the full numerical simulations of Zamankhan et al. [START_REF] Zamankhan | Steady motion of bingham liquid plugs in two-dimensional channels[END_REF] (lines) for 𝐶𝑎 * = 0.1, 0.05, and 0.025. The equilibrium film thickness to achieve a steady plug propagation is well predicted by our reduced-order model for the whole interval of Bingham numbers considered 𝐵𝑖 * ∈ [0, 1]. The differences between the model results and the numerical simulations are partially due to the difference in the geometrical set-up, as our reduced-order model assumes an axisymmetric flow, while Zamankhan et al. solve the steady liquid plug flow in a plane channel.

Parametric study on the standalone non-Newtonian non-dimensional groups

The reduced-order model formulated in Section 3 is here employed to investigate the effect of viscoelasticity and viscoplasticity. Since these two effects have been superposed in the model equations proposed in Section 3.3, we will consider them separately in the following subsections. In this parametric study, we will employ our reduced-order model to characterize the trends of 𝐿 𝑟 and 𝑡 𝑟 upon a change of 𝐶𝑎∕𝐵𝑖, and 𝐷𝑒 for physiologically-relevant regimes that exceed the parameters range investigated experimentally.

All the following results presented in Section 3.5.1 and 3.5.2 are obtained for the same baseline clean Newtonian set-up, indicated as the solid black line in Figs. [START_REF] Griese | Pulmonary surfactant, lung function, and endobronchial inflammation in cystic fibrosis[END_REF] simulations are carried out until plug rupture by integrating the ODEs of Section 3 by an explicit Euler method and employing 𝛥𝑡 = 0.001 s as time step.

Effect of viscoplasticity

The effect of viscoplasticity is taken into account by increasing the leading plug resistance and the film thickness as shown in Fig. 12. The increase of liquid plug resistance has a direct effect on the liquid plug by decelerating it. On the other hand, the increase of film thickness reduces the menisci curvatures, reducing the capillary resistance of the plug and inducing an acceleration of the plug. We stress that in this and in the following subsections we use our reduced-order model to isolate the effect of 𝐵𝑖 f ∕𝐶𝑎 f on the pressure drop and the film thickness, and we compare with Newtonian plugs keeping the dynamic viscosity constant. The goal of such an analysis is to better understand the net deceleration of the liquid plug observed for higher yield stresses (see Fig. 10). For 𝐵𝑖 f ∕𝐶𝑎 f ∈ [0, 100], the latter effect prevails and the liquid plug experiences a net acceleration. Keeping the leading film thickness constant, for a given 𝐵𝑖 f ∕𝐶𝑎 f the velocity of the leading meniscus tip is fixed and an increase of the plug velocity increases the rear capillary number. This implies that the rear film thickness increases because of both, the Newtonian (4) and the viscoplastic effects [START_REF] Veen | Recurrent exacerbations in severe asthma are associated with enhanced airway closure during stable episodes[END_REF]. For 𝐵𝑖 f ∕𝐶𝑎 f ∈ [0, 100], therefore, the viscoplastic liquid plug deposits more liquid than a Newtonian one (of same effective dynamic viscosity) due to the increased trailing film thickness and, as a consequence, the rupture occurs at a shorter rupture time 𝑡 r . Despite the fact that 𝑈 p for a viscoplastic liquid plug is higher than that for a Newtonian plug, both, the rupture time 𝑡 r and the traveled distance 𝑍 p decrease upon an increase of 𝐵𝑖 f ∕𝐶𝑎 f keeping the effective dynamic viscosity constant. This is demonstrated in Fig. 12, where the results obtained for several Bingham-to-capillary-number ratios are shown, i.e. 𝐵𝑖 f ∕𝐶𝑎 f = 0 (black, Newtonian), 𝐵𝑖 f ∕𝐶𝑎 f = 5 (red), 𝐵𝑖 f ∕𝐶𝑎 f = 10 (blue), 𝐵𝑖 f ∕𝐶𝑎 f = 50 (green), and 𝐵𝑖 f ∕𝐶𝑎 f = 100 (magenta). The two variants M1 (dashed lines) and M2 (solid lines) are also compared in Fig. 12 showing that our reduced-order model is insensitive to the two matching corrections tested.

Effect of viscoelasticity

The effect of viscoelasticity on the liquid plug rupture is taken into account by an increase of leading and trailing film thickness defined by [START_REF] Gunther | Surfactant alterations in severe pneumonia, acute respiratory distress syndrome, and cardiogenic lung edema[END_REF]. Keeping the leading film thickness constant, for a given 𝐷𝑒 f , the tip velocity of the retracting air bubble is fixed and the viscoplastic liquid plug deposits more liquid than the Newtonian one because of the increased trailing film thickness. Hence, the plug rupture occurs at shorter rupture times 𝑡 r , keeping constant the effective dynamic viscosity. Depositing more liquid on the wall, the plug reduces the curvature and, consequently, the capillary pressure resistance. This leads to increase the velocity of the plug upon an increase of 𝐷𝑒 f . Despite such an increase, the shortening of the rupture time plays a dominant role in the plug dynamics and leads to a reduction of the distance traveled by the liquid plug before rupture, 𝑍 p . This liquid plug dynamics is demonstrated in Fig. 13 for 𝐷𝑒 f ∈ [0, 100]. The results are obtained for several Deborah numbers including 𝐷𝑒 f = 0 (black, Newtonian), 𝐷𝑒 f = 5 (red), 𝐷𝑒 f = 10 (blue), 𝐷𝑒 f = 50 (green), and 𝐷𝑒 f = 100 (magenta). The liquid-gas interface is assumed to be clean, and the models M1 (dashed lines) and M2 (solid lines) are compared. The quantitative differences between the corresponding results are less than 20% and the qualitative trends of the two variants of our reduced-order model provide a consistent prediction. We stress that these results are not intended for a direct comparison with the experimental measurements since the experiments are not carried out at constant 𝛥𝑃 , 𝜇, 𝜀 f , 𝐿 p (𝑡 = 0), as we do for producing the model results in Section 3.5.

Summary and conclusions

The rupture of a non-Newtonian liquid plug in a capillary tube has been investigated for reproducing the airway reopening in the 6th-to -7th generation of adult human lungs. Four synthetic mucus samples are used to study the elastoviscoplastic effect of the liquid on the plug rupture ranging over several orders of magnitude for the Deborah, Bingham and capillary numbers. A detailed rheological characterization of our liquid samples have been provided, including the determination of the model parameters for the Herschel-Bulkley model (viscoplasticity) and the White-Metzner model (viscoelasticity). Shear relaxation time measurements were carried out, as they are commonly employed in the literature for characterizing the relaxation time of biological and synthetic mucus [START_REF] Lai | Micro-and macrorheology of mucus[END_REF][START_REF] Patarin | Rheological analysis of sputum from patients with chronic bronchial diseases[END_REF]. However, considering the extensional rheology would be also of interest as the polymeric stretching plays a significant role during the liquid plug rupture. Future studies will therefore focus on measuring the extensional viscosity and relaxation time of mucus in order to better relate the plug dynamics to the fluid rheology.

In support to our experiments, we derived an original reduced-order model built on the existing literature. Our model takes into account the effect of viscoplasticity and viscoelasticity, assuming they can be linearly superposed. By means of space-time diagrams, we quantitatively characterize the dynamics of the liquid plug along the centerline of the pipe. Upon a decrease of 𝑈 𝑝 we observe a correlated increase of the rupture time, while 𝜇, 𝜆, and 𝜏 0 are rather impactful on the dimensionless traveling length 𝐿 𝑟 ∕(𝐿 𝑖 -𝐿 𝑓 ), that increases as 𝜇 ↑, 𝜆 ↑, and 𝜏 0 ↑. These same trends are predicted by our reduced-order model, when validating against the experiments. Our reduced-order model has also been used to isolate the viscoplastic effect from the viscoelasticity of mucus, taking into account the increase of the film thickness and of the pressure drop due to viscoplasticity [START_REF] Veen | Recurrent exacerbations in severe asthma are associated with enhanced airway closure during stable episodes[END_REF].

Comparing the synthetic mucus samples at different Actigum concentrations, our experiments uncovered a remarkable axi-symmetry breaking for the liquid plugs at high 𝐵𝑖 and 𝐷𝑒. As a result, the liquid plug does not rupture along the centerline, on the contrary of what has always been observed for almost-Newtonian liquids (MUC-0.25% and MUC-0.5%) in our experiment. This non-Newtonian effect is in agreement with the experimental results of Caliman et al. [START_REF] Caliman | An experimental investigation on the Newtonian-Newtonian and viscoplastic-Newtonian displacement in a capillary tube[END_REF], who unravel the various regimes of symmetry breaking for viscoplastic fluids in capillary pipes.

Our model prediction, corroborated by the experimental evidence, shows that a standalone increase of 𝐵𝑖∕𝐶𝑎 and of 𝐷𝑒 on the liquid film thickness and plug resistance can facilitate liquid plug rupture, compared to Newtonian plugs of equal effective viscosity. This conclusions help interpreting the experimental findings by identifying the contribution of single non-dimensional groups. Upon an increase of Actigum concentration, for a given 𝑈 p , 𝐶𝑎 ↑, 𝐵𝑖 ↑, and 𝐷𝑒 ↓ (where the overline denotes the non-dimensional groups computed using 𝑈 p , see Table 2). For comparable 𝑈 p , (the same case labeling in Table 2), the dimensional rupture time is not significantly affected (the same colors in Fig. 10). This, in combination with the one-order-of-magnitude increase of 𝜇 upon an increase of concentration from MUC-0.5% and MUC-1.5% (see Table 2), explains why the non-dimensional rupture time in capillary time units (𝑡𝜎∕𝜇𝑎) decreases upon an increase of 𝐷𝑒 and a decrease of 𝐶𝑎 and 𝐵𝑖, as observed experimentally. On the other hand, the dimensional rupture time increases upon an increase of non-Newtonian effects due to a higher effective viscosity. This same trend has been reproduced by our reduced-order model (see Table 2) and is consistent with the medical literature, that reports a significant enhancement of the dimensional rupture time and plug resistance in patients affected by asthma [START_REF] Dunican | Mucus plugs in patients with asthma linked to eosinophilia and airflow obstruction[END_REF]. In fact, increasing the polymeric concentration leads to a net increased difficulty of reopening the airways and a higher pressure driving must be applied to rupture the mucus liquid plugs in patients affected by asthma in order to allow air exchange at distal airways.

These conclusive considerations suggest the future studies should focus on including the non-linear interaction between complex rheological effects, as well as on validating the predictions of our reduced-order model, either by numerical simulations or by experiments.

Fig. 1 .

 1 Fig. 1. Flow curves of the synthetic mucus solutions obtained from rheological measurement at 21 • C. (a) Viscosity 𝑣𝑠 shear rate. The dashed line represents the rheometer low-torque limit which is multiplied by 3 (𝜇 < 3𝑇 𝑚𝑖𝑛 ∕(2𝜋𝑅 3 1 ̇𝛾), with 𝑅 1 = 50 mm and 𝑇 𝑚𝑖𝑛 = 0.5 μN m [57]). (b) Shear stress 𝑣𝑠 shear rate. The colored continuous lines represents fits based on the Herschel-Bulkley model applied to flow curves. The discrepancy between fits and the experimental results is attributed to slip of solution in low torques.

Fig. 2 .

 2 Fig. 2. (a) First normal-stress difference as a function of shear rate. Data are fitted with a power law model. Black solid line of slope 𝛼 = 2 represents the Oldroyd-B model (𝑁 1 = 2𝜂 𝑝 𝜆 ̇𝛾2 ), considering an identical relaxation time at polymer contribution to the zero shear viscosity 𝜂 𝑝 . (b) Relaxation time as a function of shear rate for different concentration of synthetic mucus using the White-Metzner model. Colored continuous lines are fitted to the experimental data by a power law model. The corresponding model constants can be found in Table1.

Fig. 3 .

 3 Fig. 3. Sketch of the experimental apparatus. (a) Experimental set-up. (b) Tracking of the rear and the front menisci, and space-time diagram.

Fig. 4 .

 4 Fig. 4. Space-time diagram of the liquid plug rupture for (a, b) MUC-0.5% and (c, d) MUC-1.5%, and three initial plug lengths, i.e. (a) 𝐿 𝑝 (𝑡 0 ) = 𝐿 i = 5, (c) 13 and (b,d) 22 mm. The markers denote the post-processed data along the pipe centerline, the dashed lines are least-squares cubic fits. Top-right inset: snapshot of the liquid plug as soon as the front and the rear menisci entered the observation window. Bottom-left inset: snapshot of the liquid plug rupture (a-c) or exit of the liquid plug from the observation window. The identified contour of the liquid plug is highlighted in red.

Fig. 5 .

 5 Fig. 5. Development of an asymmetric plug propagation for MUC-1.0% and 𝐿 𝑝 (𝑡 0 ) = 𝐿 𝑖 = 13 mm.

Fig. 6 .

 6 Fig. 6. Rupture time 𝑡 𝑅 (right) and non-dimensional rupture length 𝐿 𝑟 ∕(𝐿 𝑖 -𝐿 𝑓 ) (left) plotted against: the average plug speed 𝑈 p (a, b), the effective viscosity 𝜇 (c, d), the relaxation time 𝜆 (e, f), and the yield stress 𝜏 0 (g, h). Same color coding of Fig. 2 is employed.

  04, 𝐶 2 = 0.570, 𝐶 3 = 0.233, 𝐶 4 = 0.064, and 𝐶 5 = 0.00689.

Fig. 8 .

 8 Fig. 8. The markers are derived from the resistance formula in form of infinite series (6), data extracted from fig. 8 by Che et al. [65], while the solid line is obtained from the power-law fit (7). The dashed line denotes the Hagen-Poiseuille asymptote lim 𝐿 p ∕𝑎→∞ 2𝜋𝑎 4 𝑅 c ∕𝜇𝐿 p = 16. The inset depicts a cutout of the interval 𝐿 p ∕𝑎 ∈ [0, 10].

  (a) A fully-developed flow profile in the liquid plug core; (b) Superposition of viscoelastic-and viscoplastic-induced corrections for the non-dimensional film thickness of the advancing and the retracting air fingers; (c) Correction of the capillary resistances of the liquid plug menisci determined by the viscoplastic dynamics of a Bingham fluid.

Fig. 9 .

 9 Fig. 9. Effects of viscoplasticity. The solid lines in the left panel denote the two asymptotes observed for 𝐵𝑖 * ∕𝐶𝑎 * < 20 and 𝐵𝑖 * ∕𝐶𝑎 * > 20, the dashed line includes the smoothed transition between them and depicts (9)(b), the solid lines in the right panel depict the fit (9)(a),while the markers are the data of Zamankhan et al. [47,67]. Left panel: Non-dimensional film thickness for 𝐶𝑎 = 0.25 (triangles), 𝐶𝑎 = 0.1 (squares), and 𝐶𝑎 = 0.05 (circles) taken from Zamankhan et al. [67]. Right panel: Non-dimensional pressure drop for 𝐶𝑎 = 0.025 (triangles), 𝐶𝑎 = 0.05 (squares), and 𝐶𝑎 = 0.025 (circles) taken from Zamankhan et al. [47].

(

  viii.a) Compute the core pressure drop 𝛥𝑃 𝑘 c of a Newtonian bubble employing 𝑈 𝑘 p = 𝑈 p (𝑡 = 𝑡 𝑘 ) and 𝐿 𝑘 p = 𝐿 p (𝑡 = 𝑡 𝑘 ) resulting from the previous time step; (viii.b) Assuming that the plug moves at the average velocity between the two bubble tip velocities, i.e. 𝑈 p = (𝑈 r +𝑈 f )∕2, the rear-bubble tip velocity at 𝑡 = 𝑡 𝑘 is computed as 𝑈 𝑘 r = 2𝑈 𝑘 p -𝑈 f , that leads to the rear-bubble capillary number; (viii.c) Compute the capillary number for the bubble; (viii.d) Compute the non-dimensional film thickness for a bubble advancing in a Newtonian fluid; (viii.e) Compute the pressure drop at the rear-bubble tip 𝛥𝑃 𝑘 r for a Newtonian fluid; (viii.f) Compute the Bingham and Deborah number of the rear bubble and use them to calculate the elastoviscoplastic correction 𝜀 𝑘 r = 𝜀 EVP r ; (viii.g) Compute the viscoplastic correction of 𝛥𝑃 𝑘 p by summing up all the Newtonian resistances (selecting either M1 or M2) and applying (9)(left); (viii.h) Enforce a driving pressure difference across the plug, 𝛥𝑃 ;

Fig. 11 .

 11 Fig. 11.Comparison between the equilibrium non-dimensional film thickness from the reduced-order model (markers) and the numerical simulations (lines) for 𝐶𝑎 * = 0.1 (light-blue), 0.05 (violet), and 0.025 (orange). The Bingham number is varied to determine the equilibrium conditions for the propagation of a viscoplastic liquid plug.

  and 13. The baseline parameters are: 𝜀 f = 0.05, 𝜇 = 0.005 Pa s, 𝜌 = 1000 kg/m 3 , 𝜎 = 0.070 N/m, 𝑎 = 0.005 m, 𝐿 p (𝑡 = 0) = 0.1 m, 𝑈 p (𝑡 = 0) = 0.01 m/s and 𝛥𝑃 = 150 Pa. All

Fig. 12 .

 12 Fig. 12. Effect of the viscoplastic rheological properties on the liquid plug rupture. The Newtonian case (𝐵𝑖 f ∕𝐶𝑎 f = 0) is depicted in black solid line, while the viscoplastic cases are shown in red (𝐵𝑖 f ∕𝐶𝑎 f = 5), blue (𝐵𝑖 f ∕𝐶𝑎 f = 10), green (𝐵𝑖 f ∕𝐶𝑎 f = 50) and magenta (𝐵𝑖 f ∕𝐶𝑎 f = 100). In all simulations, we set 𝜀 f = 0.05, 𝜇 = 0.005 Pa s, 𝜌 = 1000 kg/m 3 , 𝜎 = 0.070 N/m, 𝑎 = 0.005 m, 𝐿 p (𝑡 = 0) = 0.1 m, 𝑈 p (𝑡 = 0) = 0.01 m/s, 𝛥𝑃 = 150 Pa, 𝐷𝑒 f = 0 and 𝛽 = 0. The two versions of the reduced-order model for clean plugs are tested: M1 (dashed lines), M2 (solid lines). Left panels: Liquid plug dynamics in terms of plug velocity 𝑈 p , traveling distance 𝑍 p , and length 𝐿 p . Right panel: spatio-temporal diagram of the plug rear and front meniscus tips.

Fig. 13 .

 13 Fig. 13. Effect of the viscoelastic rheological properties on the liquid plug rupture. The Newtonian case (𝐷𝑒 f = 0) is depicted in black solid line, while the viscoelastic cases are shown in red (𝐷𝑒 f = 5), blue (𝐷𝑒 f = 10), green (𝐷𝑒 f = 50) and magenta (𝐷𝑒 f = 100). In all simulations, we set 𝜀 f = 0.05, 𝜇 = 0.005 Pa s, 𝜌 = 1000 kg/m 3 , 𝜎 = 0.070 N/m, 𝑎 = 0.005 m, 𝐿 p (𝑡 = 0) = 0.1 m, 𝑈 p (𝑡 = 0) = 0.01 m/s, 𝛥𝑃 = 150 Pa, 𝐵𝑖 f ∕𝐶𝑎 f = 0 and 𝛽 = 0. The two versions of the reduced-order model for clean plugs are tested: M1 (dashed lines), M2 (solid lines). Left panels: Liquid plug dynamics in terms of plug velocity 𝑈 p , traveling distance 𝑍 p , and length 𝐿 p . Right panel: spatio-temporal diagram of the plug rear and front meniscus tips.

  Fig.2shows the results of the normal stress difference and relaxation time as a function of shear rate. Note that 𝑁 1 , 𝜎 and 𝜂( ̇𝛾)

		𝜂( ̇𝛾) and 𝜆( ̇𝛾) are the shear rate dependent total viscosity and
	relaxation time, respectively, and 𝜂 𝑠 is the solvent viscosity. The dif-
	ference, 𝜂( ̇𝛾) -𝜂 𝑠 , gives the polymer viscosity contribution 𝜂 𝑝 . Finally,
	we can calculate the relaxation time by rewriting (1) to isolate 𝜆( ̇𝛾) by
	introducing the ratio of the first normal-stress difference 𝑁 1 to shear
	stress 𝜎 = 𝜂( ̇𝛾)	̇𝛾
	𝜆( ̇𝛾) =	1 2 ̇𝛾	𝜂( ̇𝛾) 𝜂( ̇𝛾) -𝜂 s (	)	𝑁 1 𝜎	(2)

Acknowledgments

Support from National Institutes of Health (NIH), USA, grant number R01-HL136141, and the Scientific and Technical Research Council of Turkey, (TUBITAK), grant number 119M513, is kindly acknowledged.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix. Conservation laws

The second and third equations of (3) are here derived starting from first principles. The cross-section-averaged conservation of momentum in the liquid plug reads

where 𝑉 p is the plug volume approximated as 𝑉 p = 𝜋𝑎 2 𝐿 p . Simplifying the constant factor common to all terms of the equation, i.e. 𝜋𝑎 2 , and dividing by the liquid density, it yields

Finally, re-arranging the l.h.s. and substituting the liquid plug pressure drop with the corresponding plug resistance times the flow rate (𝛥𝑃 p = 𝜋𝑎 2 𝑅 p 𝑈 p ), the one-dimensional Navier-Stokes equation reads

which is equivalent to the second equation of ( 3), where the first term of the l.h.s. has been moved to the r.h.s. and we divided all the term of the equation by 𝐿 p .

The last equation of ( 3) is derived by the cross-section-averaged conservation of mass for the liquid phase. Considering a control volume between the trailing and the leading film, the one-dimensional conservation of mass reads

Substituting the plug volume with 𝑉 p = 𝜋𝑎 2 𝐿 p and dividing by 𝜌𝜋𝑎 2 , the conservation of mass reads

which is equivalent to the last equation of (3).


Valentin Noël

Thomas Rodet

Identification and characterization of damaged fiber-reinforced laminates in a Bayesian framework

Valentin NOËL a,b , Thomas RODET a , and Dominique LESSELIER b,1 a Systèmes et Applications des Technologies de l'Information et de l'Energie, Ecole Normale Supérieure Paris-Saclay, Gif-sur-Yvette, France b Laboratoire des Signaux et Systèmes, Université Paris-Saclay, CNRS, CentraleSupélec,Gif-sur-Yvette, France

Abstract. Non-destructive testing of damaged composite laminates modeled from the homogenization of fiber-reinforced polymers has always been and still is a challenge, both by its underlying complexity and the difficulties encountered in the quantification of uncertainties related to the identification and characterization of defects. In order to provide a rigorous framework that accepts data from different modalities and allows data fusion as well, we propose here a Bayesian neural network (BNN) [START_REF] Kononenko | Bayesian neural networks[END_REF] with two input streams, with a focus on local inter-layer delaminations identification and characterization.

Keywords. thermography, composites, Bayesian neural network, data fusion

Composite materials are widely used in numerous fields, ranging from the automotive industry to green engineering, due to their low cost and lightness. Additionally, their superior strength/ and stiffness/weight ratios vs. traditional materials make them appealing. Yet their non-destructive testing may be highly challenging, with the need for more modeling and computational imaging, before envisaging sound laboratory-controlled experimentations.

Active thermography as a contactless procedure is more and more used for the inspection of isotropic workpieces, refer to [START_REF] Ratsakou | Shape reconstruction of delamination defects using thermographic infrared signals based on an enhanced Canny approach[END_REF] and references therein, but an increasing volume of contributions appears on fibered laminates, with a strong emphasis on carbonfiber-based ones, those exhibiting high anisotropic behavior due to their thermal conductivity function of orientation, e.g., [START_REF] Bard | Influence of fiber volume content on thermal conductivity in transverse and fiber direction of carbon fiber-reinforced epoxy laminates[END_REF], as illustrated in the recent review [START_REF] Ciampa | Recent advances in active infrared thermography for non-destructive testing of aerospace components[END_REF].

The above stated, focused on thermographic modalities of inspection, one needs (i) proper mathematical and physical models of damaged laminates, (ii) versatile enough codes that can provide a good amount of accurately simulated data when those assumed to be inspected by both modalities, and last but not least, (iii) efficient computational imaging tools that can use in an optimized fashion the information provided by the data. The three-fold challenge is strong and from the start one needs simplifications.

The first simplification is that intact laminates are homogeneous per horizontal ply and exhibit within each ply a uniaxial anisotropic behavior, which is manifested by a transverse and a longitudinal thermal conductivity.

The second simplification is that damages generally can be seen as thin flat delaminations occurring between plies, which are amenable in thermography (in the case of a flash lamp as a heat source at least) to an extension of the Truncated Region Eigenfunction Expansion (TREE).

The third simplification is that one has available imaging tools that can cope with data provided by the modality(ies). Data fusion can be approached in the context of an active thermography procedure with the use of two (or more) IR cameras capturing the temperature evolution on the top and bottom surfaces of the part, which can be considered as two fusionable data streams, making data fusion more readily available. So, if data simulated on a host of damages, double-stream convolutional neural networks (CNN) appear as the way forward, also having insurance that brute force codes (FEM) could always provide good data, if the semi-analytical algorithms were to fail in some cases. This approach was shown highly effective in [START_REF] Qin | Breast imaging by convolutional neural networks from joint microwave and ultrasonic data[END_REF].

However, because of the large number of parameters included in a deep learning model, many questions arise about their initialization, while other questions also arise about the quantification of uncertainties as well as the existence of a physical framework in what may look like a black box. To tackle this problem, BNNs have become highly attractive due to their ability to quantify uncertainties [START_REF] Moloud | A review of uncertainty quantification in deep learning: Techniques, applications and challenges[END_REF]. Thus, the combination of knowledge related to the data and the Bayesian model should allow a rigorous and precise identification, while ensuring a detailed characterization, reinforced by quantification of the uncertainties.

To conclude, this investigation aims at pointing out that semi-analytical solutions are indeed capable at least in simple enough cases (thin inter-layer delaminations) to yield proper data to mimic thermographic inspections, once accounted for homogenization of the intact structures. Then, data can be used and/or fused like was done by the authors and co-workers in breast imaging [START_REF] Qin | Breast imaging by convolutional neural networks from joint microwave and ultrasonic data[END_REF], [START_REF] Qin | Fused microwave and ultrasonic breast imaging within the framework of a joint variational Bayesian approximation[END_REF] but in a Bayesian realm enabling uncertainty quantification.

Corresponding Author: dominique.lesselier@centralesupelec.fr.


Valentin Noël

Thomas Rodet

Dominique Lesselier

Data fusion and non-destructive testing of damaged fiber-reinforced laminates

Keywords: Electromagnetics, thermography, composites, data fusion

Data fusion and non-destructive testing of damaged fiber-reinforced laminates

Valentin Noël, Thomas Rodet, Dominique Lesselier

Composite materials are widely used in automotive industry, aeronautics, and green engineering, since low cost, light weight, and of large strength/ and stiffness/weight ratios vs. traditional materials. Yet their non-destructive testing may be highly challenging, with need for more modeling and computational imaging, before envisaging sound laboratorycontrolled experimentations.

Active thermography as a contactless procedure is more and more used for inspection of isotropic work pieces, refer to [START_REF] Ratsakou | Shape reconstruction of delamination defects using thermographic infrared signals based on an enhanced Canny approach[END_REF] and references therein, but an increasing volume of contributions appears on fibered laminates, with strong emphasis on carbon-fiberbased ones -those exhibiting strong anisotropic behavior due to their thermal conductivity function of orientation, e.g., [START_REF] Bard | Influence of fiber volume content on thermal conductivity in transverse and fiber direction of carbon fiber-reinforced epoxy laminates[END_REF]-as illustrated in the recent review [START_REF] Ciampa | Recent advances in active infrared thermography for non-destructive testing of aerospace components[END_REF], with deep learning coming also to the fore [START_REF] Tong | Automated defect detection in CFRP using flash infrared thermography and deep learning method, 2nd European NDT & CM Days[END_REF]. In rather sharp contrast, electromagnetic probing of such composite structures appears much less studied, possibly since standard eddycurrent probes tend not to work so well for carbon-fiber-based laminates, and may call for tailor-made sensors [START_REF] Heuer | Review on quality assurance along the CFRP value chain -Nondestructive testing of fabrics, preforms and CFRP by HF radio wave techniques[END_REF] in the MHz range and beyond, so as usual eddy-current hypotheses fail.

The above stated, focused on thermographic and electromagnetic modalities of inspection, one needs (i) proper mathematical and physical models of damaged laminates, (ii) versatile enough codes that can provide a good amount of accurately simulated data when those assumed to be inspected by both modalities, and last but not least, (iii) computational imaging tools that can fuse in a proper way the information which those data provide independenty, with the benefit of absence of registration since the structure is unchanged during testing. The three-fold challenge is strong and simplifications are needed.

The first simplification is that intact laminates are homogeneous per horizontal ply and exhibit within each ply an uniaxial anisotropic behavior, which is manifested by a transverse and longitudinal dielectric permittivity or thermal conductivity, being underlined that one has the joint benefit that first-order electromagnetic homogenization [START_REF] Li | On homogenization of fiber-reinforced polymers[END_REF] holds close to the same in thermography.

The second simplification is that damages generally can be seen as thin flat delaminations occuring between plies, which is not too demanding for the electromagnetic Green-dyadic-based semi-analytical approaches in [START_REF] Zhong | Fast calculation of scattering by 3-D inhomogeneities in uniaxial anisotropic multilayers[END_REF], while they are amenable in thermography (in the case of a flash lamp as heat source at least) to extension of the Truncated Region Eigenfunction Expansion (TREE), so far applied to isotropic laminates [START_REF] Ratsakou | Shape reconstruction of delamination defects using thermographic infrared signals based on an enhanced Canny approach[END_REF], since only the heat equation is impacted by the uniaxial behavior.

The third simplification is that one has available imaging tools that can cope with data provided by the two modalities. So, if data simulated on a host of damages, doublestream convolutional neural networks (CNN) appear the way forward, also having insurance that brute force codes (FEM) could always provide good data, if the semi-analytical algorithms were to fail in some cases. This approach was shown highly effective in [START_REF] Qin | Breast imaging by convolutional neural networks from joint microwave and ultrasonic data[END_REF] for fused ultrasound and microwave imaging, if enough physics within it to guide to sound solutions, with the additional advantage that a stand-alone solution with a single modality can always be developed. With now the advantage of far less data needed than with CNN, the Bayesian framework as also devised for breast imaging in [START_REF] Qin | Fused microwave and ultrasonic breast imaging within the framework of a joint variational Bayesian approximation[END_REF] enables in principle to bypass the need of many priors on regularisation factors (those are produced in effect during the construction of the maps of the damages) while there should be good profit in the Bayesian formalism to draw from the assumption of thin delaminations.

To conclude, this preliminary investigation aims at pointing out that data fusion should enable testing of carbon-fibered-based laminates. Semi-analytical solutions are indeed capable at least in simple enough cases (thin inter-layer delaminations) to yield proper data to mimick electromagnetic and thermographic inspections, once the intact structures properly homogenized. Then, data from both modalities could be fused, e.g., in deep learning and Bayesian realms, as inspired, e.g., from [START_REF] Qin | Breast imaging by convolutional neural networks from joint microwave and ultrasonic data[END_REF][START_REF] Qin | Fused microwave and ultrasonic breast imaging within the framework of a joint variational Bayesian approximation[END_REF].

Corresponding Author: dominique.lesselier@centralesupelec.fr.


Juan Forero-Saboya

Parth Desai

Roman Healy-Corominas

Encarnacion Raymundo-Piñero

Aurélien Canizarès

Dominique Foix

Jean-Marie Tarascon

Sathiya Mariyappan

à la diffusion de documents scientifiques 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

7

  Figure S1. Recorded temperature inside 18650 cells during the three formation cycles, depending on the temperature employed during formation. A) 0 °C, B) 25 °C, and C) 55°C. The peaks associated to heat release events in the first charge are marked with dashed lines, and the resting regions between charge and discharge are highlighted in light blue.

Figure S2 .Figure S3 .

 S2S3 Figure S2. Fitting of the obtained PEIS spectrum right after formation (measured at 100 % SOC and 10 °C). Points represent the obtained experimental data and the black line is the fitted spectrum using the equivalent circuit shown as inset.

Figure S4 .

 S4 Figure S4. Performance metrics of the NVPF|HC cells as depending on the formation protocol used. A) Coulombic efficiency evolution upon cycling, and B) Evolution of polarization (calculated as the difference between the average voltage during charge and the average voltage during discharge)

Figure S6 .

 S6 Figure S6. Effect of the temperature on a Raman band of EC. No significant displacement is observed in a solvent upon changes in temperature.

Figure S7 .

 S7 Figure S7. XPS spectra of NVPF electrodes' surface after formation.

Figure S8 .

 S8 Figure S8. Formation of 18650 NVPF|HC cells. a) Cycling profiles during formation for the 55°C x 3 formation protocol, b) Cycling profiles during the 55-25-0 °C formation protocol, c) Nyquist plot of the cell impedance measured right after formation, at 25 °C and 0 % state of charge (SOC), and d) First cycle of both cells right after formation.

Figure S9 .

 S9 Figure S9. Evolution of the cycling curves of the cells during cycling at different temperatures.

Table S2 . Integrated intensity of the deconvoluted XPS spectra showed on Figure 5 and Supplementary Figure S6

 S2

				C 1s			F 1s	Na 1s		P 2p		V 2p
	HC electrodes	HC 51 5 9,3 55°C Pt1 0,2 Pristine 0°C 25°C	PVdF Deposited C PVdF + POxFy 24 -20 6,7 5,1 11,3 10 8,4 13,6 0,6 8,2 8,3	NaF 3,4 25,1 20,3 25,8	Na -37,6 29,6 28	134,0 eV -0,7 0,6 4,3	135,3 eV -0,5 0,3 -	137,8 eV -0,8 0,6 0,9	V ---4,7
		55°C Pt2 1,1	1,4	5,5	9,5	25	37,9	0,9	1,7	1,2	0,8
			C	PVdF Deposited C PVdF +	NaF	NVPF 133,8	135,3	136,9	V NVPF
	NVPF electrodes	Pristine 0°C 25°C	black 38 36,4 37,1	22 19 19,1	---	POxFy 23 23,2 23,1	and NVPF 1,9 2,7 2,5	2,7 2,6 2,5	eV NVPF 2,9 2,7 2,7	eV -0,4 0,3	eV -0,5 0,4	1,5 1,4 1,6
		55°C	38,1	18,6	-	21,5	2,4	2,3	2,8	0,7	-	1,6

Supporting Information

Influence of formation temperature on cycling stability of sodium-ion cells: A case study of Na 3 V 2 (PO 4 ) 2 F 3 | HC cells.