I Charfeddine

J C Majesté

C Carrot

O Lhost

A MODEL FOR THE PREDICTION OF THE MORPHOLOGY OF IMMISCIBLE BLENDS OF POLYMERS

Keywords: Polyolefin blends, Morphology prediction, Co-continuity

This paper presents a micro-rheological model to predict the morphology of immiscible polymer blends, established during melt processing, on the whole composition range. This model based on elementary mechanisms of dispersion described in the literature, uses two different descriptions of the topology in the case of dispersed and co-continuous phases. By using a small number of material and process parameters (shear rate, viscosity of the phases and interfacial tension) and by analysing the stability of the droplets or threads of one phase in the matrix, it enables to calculate diagrams of continuity, sizes of phases and interfacial areas.

The model was tested in comparison with experimental results on selected polyolefin blends, regarding the type of morphology, the interfacial area and the excess of elasticity at low frequency in oscillatory shear.

INTRODUCTION

Polyolefins represent the largest commodity thermoplastics used today because of their low cost, high-speed production, high mechanical performance and good barrier properties. The polyethylene (PE) and the polypropylene (PP) are the most produced polymers. They are used in a wide variety of applications such as packaging, automotive and medical fields. For the PE, the control of the molecular chain structures using the catalyst system and the polymerization process impacts the branching density LCB or SCB (long or short chain branching), the weigh average molecular weight M w and the molecular weight distribution.

The choice of the catalysis system allows the control of the tacticity of the PP (isotactic or syndiotactic) and its crystallinity. Hence, the development of improved catalyst systems such as single site metallocene catalysts allows the production of a resulting polymer with a more uniform molecular structure [START_REF] Qiao | Recent advances in polyolefin technology[END_REF]. In general, isotactic polypropylene is used for application requiring high modulus. Polyethylene is used for application demanding higher resistance to brittle fracture. Though, olefinic homopolymers or copolymers can hardly fulfill all the requirements and blends of PE and PP may be a useful solution to extend their application domain. A large number of scientific works have been dealing with the more general topics of blends of olefinic polymers [START_REF]Polyolefin blends[END_REF].

As for any type of immiscible polymer blend, depending on the concentration [START_REF] Utracki | On the viscosity concentration dependence of immiscible polymer blends[END_REF], two main types of morphologies can be observed: dispersed and co-continuous morphology. They are responsible for different functional properties: mechanical, barrier or electrical. For example, the addition of PE or olefinic elastomers dispersed in a PP matrix enhances the impact resistance. It was shown that PE/PP blends in co-continuous morphology may have a better Young modulus than that of homopolymers. The prediction of the morphology development is therefore of great interest and it is the subject of micro-rheology. Lots of theoretical aspects of this issue have been dealt with in the literature during the last decades. These studies especially focus on the morphology control during melt processing in single and twin-screw extruders or internal mixers and are reported in reviews [START_REF]Mixing and compounding of polymers: theory and practice[END_REF][START_REF] Utracki | Development of polymer blend morphology during compounding in a twin-screw extruder. Part I: Droplet dispersion and coalescence-a review[END_REF][START_REF] Tokita | Analysis of Morphology Formation in Elastomer Blends[END_REF]. They generally start from fundamental knowledge developed on blends of diluted Newtonian liquids for which the coalescence mechanism is neglected or for viscoelastic systems for which the volume fraction of the dispersed phase is below the percolation threshold. Under these assumptions, it is possible to figure out the fundamental deformation and breakup mechanisms during the flow in order to understand the evolution of the morphology during processing.

In the literature, many scientific papers [START_REF] Willemse | Co-continuous morphologies in polymer blends: the influence of the interfacial tension[END_REF][START_REF] Janssen | Droplet breakup mechanisms: Stepwise equilibrium versus transient dispersion[END_REF][START_REF] Elemans | Transient phenomena in dispersive mixing[END_REF] illustrate that the blend morphology depends mainly on a limited number of parameters: the interfacial tension, the kinematics of the flow, the rheological properties of the components. Among the large number of experimental and theoretical studies, a few milestones of work on this understanding are worth recalling. For dilute Newtonian systems, Taylor [START_REF] Taylor | The Viscosity of a Fluid Containing Small Drops of Another Fluid[END_REF][START_REF] Taylor | The Formation of Emulsions in Definable Fields of Flow[END_REF] showed that the evolution of the morphology is governed by the viscosity ratio and the capillary number that defines the equilibrium drop size through a balance between the interfacial and viscous forces. The comparison of the capillary number () to the critical capillary number ( cr ) is useful to predict whether a deformed drop might disintegrate into small droplets or might keep an equilibrium shape and to predict the smallest drop size that can be reached under flow if the breakup mechanism is based on the Taylor mechanism. The variation of the critical capillary number as a function of the viscosity ratio was examined experimentally by Grace [START_REF] Grace | Dispersion phenomena in high viscosity immiscible fluid systems and application of static mixers as dispersion devices in such systems[END_REF] in shear or hyperbolic (extensional) flow.

On this basis and from its own work, De Bruijn proposed empirical relationships for this variation in shear and extension [START_REF] De Bruijn | Deformation and breakup of drops in simple shear flows[END_REF]. Rumscheidt and Mason [START_REF] Rumscheidt | Particle motions in sheared suspensions XII. Deformation and burst of fluid drops in shear and hyperbolic flow[END_REF] sorted the break-up mechanisms under shear and extension for couple of fluids having different viscosity ratios and a wide range of interfacial tensions. In shear flow, if the capillary number is much larger than the critical capillary number, droplets might extend until micron size and then break into a series of a small identical droplets. This breakup mechanism initially reported by Rayleigh for gases, was analyzed by Tomotika [START_REF] Tomotika | On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous fluid[END_REF][START_REF] Tomotika | Breaking up of a drop of viscous liquid immersed in another viscous fluid which is extending at a uniform rate[END_REF] for viscous liquids and it occurs during flow but also while returning at rest. Sheeting in shear flows was also described as possible mechanism to get dispersed phase morphology with final sizes independent on the capillary number [START_REF] Sundararaj | Sheet formation in immiscible polymer blends: model experiments on initial blend morphology[END_REF][START_REF] Willemse | Morphology development in immiscible polymer blends: initial blend morphology and phase dimensions[END_REF]. Coalescence and its competition with particle break-up have also been studied on model systems especially [START_REF] Elmendorp | A study on polymer blending microrheology: Part IV. The influence of coalescence on blend morphology origination[END_REF][START_REF] Minale | Study of the morphological hysteresis in immiscible polymer blends[END_REF][START_REF] Ziegler | On the role of drop break up for coalescence processes and morphology development in polymer blends under shear[END_REF][START_REF] Vinckier | Droplet size evolution during coalescence in semiconcentrated model blends[END_REF].

For a real-life system, the parameters influencing the final morphology are more numerous. In addition to the viscosity of the matrix phase, viscosity ratio, interfacial tension and deformation rate, other parameters such as the composition especially in concentrated systems, the flow history [START_REF] Minale | Effect of shear history on the morphology of immiscible polymer blends[END_REF] during the mixing and the elasticity of the components strongly interfere [START_REF] Wu | Formation of dispersed phase in incompatible polymer blends: Interfacial and rheological effects[END_REF]. Therefore, the prediction of the morphology of polymer blends in the melt state has interested many research groups in the 1980s and 1990s. The goal was to combine fluid mechanics with the previous micro-rheological models including the balance between coalescence and the break-up mechanism. Tokita [START_REF] Tokita | Analysis of Morphology Formation in Elastomer Blends[END_REF] developed a model to predict the equilibrium particle size of a natural rubber NR dispersed in a matrix of ethylene-propylenediene (EPDM) during mixing in an internal mixer. A balance between breaking mechanism and coalescence in a simple shear flow is supposed to be the origin of the equilibrium particle size. Obtaining a small particle size requires an increase of the shear rate, a decrease of the volume fraction of the dispersed phase and a low interfacial tension. On the basis of Tokita's work, Castro [START_REF] Castro | Experimental and theoretical description of low frequency viscoelastic behaviour in immiscible polymer blends[END_REF] derived a model to predict the interfacial area on a larger composition range including those giving rise to a co-continuous morphology. This enables to calculate the boundaries of the dispersed or co-continuous domains. Utracki and Shi [START_REF] Utracki | Development of polymer blend morphology during compounding in a twin-screw extruder. Part I: Droplet dispersion and coalescence-a review[END_REF][START_REF] Utracki | Development of polymer blend morphology during compounding in a twin-screw extruder. Part II: Theoretical derivations[END_REF] proposed a model to predict the morphology evolution of binary polymer blends during extrusion in twin screw extruders. A shear flow with  > 4  cr is assumed, giving rise to large deformation, for which the drops are extended into threads and the model considers two break-up mechanisms: either the threads disintegrate into small droplets when the flow is stopped or the break-up occurs under shear during the extrusion. The flow conditions in the various zones of the extruder (screw design, speed of rotation, temperature, filling) and the properties of the materials (viscosity, interfacial tension) are used to calculate the local pressure and deformation profile to define the average diameter of a dispersed phase. Huneault et al. [START_REF] Huneault | Development of polymer blend morphology during compounding in a twin-screw extruder. Part IV: A new computational model with coalescence[END_REF] extended the model of Utracki and Shi for concentrated dispersed phases taking into account a balance with coalescence with an approach similar to that of Tomotika. Additionally, breakup occurs either when the thread diameter decreases below the sub-micron range (critical diameter), or in the screw region where the stress is low. This makes the development of the morphology more sensitive to the screw profile than to the speed of rotation in the extruder.

Such results were also used in the "2-zones" model by Janssen et Meijer [START_REF] Janssen | Dynamics of liquid-liquid mixing: A 2-zone model[END_REF].

Most of the models were established in the case of dispersed phase morphology. The cocontinuous morphology was rarely taken into consideration in the models of the literature. In the composition diagram proposed by Castro, the intermediate compositions are considered as an extension of the domain of the dispersed phases. In the Tokita model, the probability of coalescence and the energy of break are poorly defined. Co-continuity is often described by a phase inversion criterion, comparing the viscosity of emulsions A in B or B in A and taking the lowest viscosity under flow for selection [START_REF] Utracki | On the viscosity concentration dependence of immiscible polymer blends[END_REF][START_REF] Epinat | Predicting phase inversion based on the rheological behavior in Polyamide 6/Polyethylene blends[END_REF]. Indeed, if the descriptors of the dispersed phase morphology are trivial, this is not the case for the co-continuous morphology. Among the attempts in the literature, it is worth mentioning the use of the interfacial area and of a tensor of surface orientation by Lee and Park [START_REF] Lee | Rheology and dynamics of immiscible polymer blends[END_REF], the use of gyroids by Sakurai et al. [START_REF] Sakurai | Gyroid structures and morphological control in binary blends of polystyrene-block-polyisoprene diblock copolymers[END_REF] and the geometrical models of Yu [START_REF] Yu | Linear viscoelasticity of polymer blends with cocontinuous morphology[END_REF]. Willemse et al. [START_REF] Willemse | Co-continuous morphologies in polymer blends: a new model[END_REF] considered the minor phase of a fully co-continuous blend as a dispersion of randomly oriented rods at their maximum packing density in the major phase. Combination with micro-rheological requirements for the existence of such rods, results in an equation giving a range of compositions within which fully co-continuous structures can exist.

In this paper, a micro-rheological model that allows to predict the morphology development of immiscible polymer blends during extrusion is proposed. The main goal is to use morphology descriptors for the dispersed phase and co-continuous morphologies in the whole composition range, to build a continuity diagram and to predict the type of morphology and its characteristic size at any composition after melt blending. The model presented hereafter proposes a scenario including some elements from the previous theories and includes a rational determination of the limits of morphological domains to predict the characteristic size of the morphology. At last, a comparison between the model prediction and experimental results will be provided.

THEORETICAL MODEL

Descriptors for the topology

The dispersed phase morphology is assimilated to spherical droplets inserted in a continuous matrix. For a volume fraction , the interfacial area per unit volume is given by Q d :

Q d = 3Φ R v (1)
Where: R v is the volume average radius of particles.

As mentioned before, the descriptors of the co-continuous topology were borrowed from an idea of Yu [START_REF] Yu | Linear viscoelasticity of polymer blends with cocontinuous morphology[END_REF] to describe the viscoelastic properties of polymer blends with co-continuous morphology. The co-continuous morphology is described by the duplication of an elementary cell with characteristic dimension l c , made of 3 parallelepipeds with orthogonal faces (Figure 1). l c is also the length of the parallelepipeds and (a) is the width of their square section. (a) and l c are related by the volume fraction  through the equation:

Φ=3 ( a l c ) 2 -2 ( a l c ) 3 (2)
The interfacial area per unit volume of co-continuous structure Q c is calculated as: In the case of a dispersed phase morphology under flow, if breakup is driven by a stepwise mechanism, the diameter of stability (D 0 ) of diluted droplets can be calculated from the Taylor analysis [START_REF] Taylor | The Viscosity of a Fluid Containing Small Drops of Another Fluid[END_REF][START_REF] Taylor | The Formation of Emulsions in Definable Fields of Flow[END_REF] in comparison to the critical capillary number [START_REF] Grace | Dispersion phenomena in high viscosity immiscible fluid systems and application of static mixers as dispersion devices in such systems[END_REF] by:

Q c = 12 a ( l c -a ) l c 3 (3)

Scenarios for the morphology development

D 0 =2 σκ cr ɳ m γ ( 4
)
Where γ is the shear rate,  is the interfacial tension,  m is the matrix viscosity, and  cr in shear may be calculated by [START_REF] De Bruijn | Deformation and breakup of drops in simple shear flows[END_REF]:

log κ cr =-0.506-0.0994logp+0.124(logp) 2 - 0.115 logp-0.611 ( 5
)
where p is the droplet to matrix viscosity ratio.

When the concentration increases, coalescence happens and the equilibrium diameter resulting from the balance between breakup and coalescence increases with the concentration as [START_REF] Huneault | Development of polymer blend morphology during compounding in a twin-screw extruder. Part IV: A new computational model with coalescence[END_REF]:

D=D 0 ( 1+ ( Φ Φ * ) 4 3
)

The value of the coalescence threshold of 1% agrees well with experimental results [START_REF] Sundararaj | Sheet formation in immiscible polymer blends: model experiments on initial blend morphology[END_REF][START_REF] Huneault | Development of polymer blend morphology during compounding in a twin-screw extruder. Part IV: A new computational model with coalescence[END_REF] or with adaptation of the coalescence theories to take into account the hydrodynamic interactions in concentrated systems [START_REF] Ziegler | On the role of drop break up for coalescence processes and morphology development in polymer blends under shear[END_REF][START_REF] Vinckier | Droplet size evolution during coalescence in semiconcentrated model blends[END_REF].

However, because the diameters may become large because of coalescence, and also because the capillary number may be considerably larger than the critical value, the particles may also be deformed into linked stable filaments of diameter (a) which give rise to a co-continuous morphology during mixing. By considering a primitive cubic cell with dimension l c , containing a droplet of diameter D, that evolves towards a pattern such that of Figure 1, volume conservation allows to write:

Φ= πD 3 6l c 3 =3 ( a l c ) 2 -2 ( a l c ) 3 (7)
where:

l c =D ( π 6Φ ) 1/ 3 (8)
and dimension (a) can be calculated from the third order equation. The analytical solution is:

a= l c 2 ( 1-cos ( θ 3 ) +√3sin | θ 3 |) (9)
where:

θ=arctan

( 2√Φ-Φ 2 1-2Φ ) (10)
The co-continuous morphology remains stable if (a) > a * with a * in the order of micron size (Path 3 in Figure 2). On the contrary, if (a) < a * , filaments disintegrate into small droplets by capillary instabilities and the resulting droplet diameter is in the order of a * (Path 2 in Figure 2). This is similar to the proposal of Willemse [START_REF] Willemse | Co-continuous morphologies in polymer blends: stability[END_REF] who suggested that the stability of the filament, and therefore of the co-continuous morphology, is obtained when the dominant wavelength of capillary instabilities is larger than l c .

In the intermediate zone, at concentrated fractions of both components, when D > a * and (a) > a * , the morphology is co-continuous with characteristic sizes l c and (a).

At low volume fraction of one component, when a * > D, the droplet diameter is more likely to be equal to a * (Path 1 in Figure 2). The occurrence of the capillary instabilities was strongly suggested by experimental observations as shown later in the paper.

(a * ) is an important parameter however it remains difficult to determine this critical dimension with precision. It can be theoretically calculated from the deformation analysis of Tomotika [START_REF] Tomotika | On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous fluid[END_REF][START_REF] Tomotika | Breaking up of a drop of viscous liquid immersed in another viscous fluid which is extending at a uniform rate[END_REF] as suggested by Janssen and Meijer [START_REF] Janssen | Droplet breakup mechanisms: Stepwise equilibrium versus transient dispersion[END_REF][START_REF] Janssen | Droplet breakup mechanisms: Stepwise equilibrium versus transient dispersion[END_REF] in the case of thread breakup in quiescent conditions or flows with very low shear rate for Newtonian fluids or in extensional flows with constant stretching rate [START_REF] Janssen | Droplet breakup mechanisms: Stepwise equilibrium versus transient dispersion[END_REF]. However, these authors also mention that (a * ) can hardly be calculated in simple shear flows for which the stretching rate is continuously changing. For these reasons, in this paper, the value of 1 μm often reported in the literature [START_REF]Mixing and compounding of polymers: theory and practice[END_REF][START_REF] Huneault | Development of polymer blend morphology during compounding in a twin-screw extruder. Part IV: A new computational model with coalescence[END_REF][START_REF] Sundararaj | Morphology development in polymer blends[END_REF] was used for all blends. Actually, at this critical value, It was further assumed that whenever the condition for capillary instabilities (a < a*) was fulfilled, there was enough time for breakup of the filaments in the extruder.

Inputs, calculations and outputs for the micro-rheological model

The model requires the knowledge of a limited number of data. The mean shear rate in the processing machines γ is used to describe the flow in the twin-screw extruder. The model does not describe the flow locally contrary to other models in the literature [START_REF] Utracki | Development of polymer blend morphology during compounding in a twin-screw extruder. Part II: Theoretical derivations[END_REF][START_REF] Huneault | Development of polymer blend morphology during compounding in a twin-screw extruder. Part IV: A new computational model with coalescence[END_REF][START_REF] Janssen | Dynamics of liquid-liquid mixing: A 2-zone model[END_REF]. The viscosity of the matrix  m and that of the dispersed phase  d are measured at the mean shear rate which is also necessary to define the viscosity ratio p. The interfacial tension  is an important parameter for the calculation of the capillary number.

All calculations were performed with Scilab. The various quantities in Equations 4-10 were first calculated on both sides of the composition diagram. Firstly, dispersed phase morphology is assumed and the equilibrium breakup diameter including coalescence (D) is calculated by Equations 4-6. Secondly, assuming co-continuous morphology, the critical dimensions (l c and a) are calculated with Equation 7-10.

Then, the final dimensions of the phases (R v , l c , a) were obtained from the results of both calculations on each side of the phase inversion point at which they converge. For this, the most stable structure is selected after comparing the characteristic size to the critical diameter a * that was set at 1 µm. If D < a * , a dispersed phase morphology with size R v = a * /2 is expected, the level of continuity is zero for the minor phase and 1 for the major phase. If D > a * and a > a * , the level of continuity is one for both phases, a co-continuous morphology with characteristic sizes l c and (a) will be observed. Between these two well defined morphologies, on each side of the composition diagram, there is an intermediate zone where both structures, filaments and drops, can be stable and coexist. In this zone, an arbitrary rule of moments was applied to determine the amount of each structure as well as the continuity level d c :

d c = D -a * D -a (11)
This continuity level was used to calculate the total interfacial area according to the mixing rule:

Q=Q d ( 1-d c ) +Q c d c ( 12
)
Where Q d and Q c are given by Equations ( 1) and ( 3) respectively.

The outputs are the type of morphology (dispersed, co-continuous), the characteristic dimensions of the dispersed phase (R v ) or co-continuous morphologies (l c and a), the continuity diagram and the interfacial area Q as a function of the composition.

MATERIALS AND EXPERIMENTS

Materials and blends

Polyethylenes and polypropylenes used in this study were supplied by the Total company (Feluy, Belgium). The characteristics are given in Table 1

Rheological characterization

The rheological study was carried out on a rotating strain-controlled rheometer (ARES, TA Waters) with parallel plane geometry (diameter 25mm). All rheological measurements were performed under nitrogen atmosphere at 190°C which is the processing temperature.

Dynamic frequency sweep tests were carried out from 100 to 0.01 rad/s in the linear domain.

The stability of the polymers and blends was checked and the content of stabilizer mixture was adapted for the duration of the test for LDPE.

Interfacial tension measurements

The interfacial tension between LDPE/iPP and HDPE/sPP was measured on blends containing 20% of dispersed minor phase, using the Gramespacher and Meissner model [START_REF] Gramespacher | Interfacial tension between polymer melts measured by shear oscillations of their blends[END_REF] and looking at the relaxation time of the interface τ interface . The interfacial tension can be calculated from this time:

τ interface = ɳ m R v σ (19 p+16 )(2 p+3) 40 ( p+1) ( 1+Φ 5(19 p+16) 4( p+1)(2 p+3) ) (13)
Where: p = viscosity ratio between the dispersed phase and matrix viscosity, It was shown that R v /R n was smaller than 2 for all blends which is a necessary condition for the use of the volume average radius in Equation 13. Then the relaxation time of the interface was taken from a fit of the imaginary part of the complex viscosity by a combination of three Cole-Cole equation for the viscoelastic contributions from the matrix, dispersed phase and interface. It was finally checked that the values of interfacial tensions are consistent for the 80/20 and 20/80 compositions. The results agree well with calculation carried out alternatively by the Palierne method [START_REF] Palierne | Linear rheology of viscoelastic emulsions with interfacial tension[END_REF].

 m =

Scanning electron microscopy SEM

SEM observations were carried out using a FEI NovananoSEM 200 instrument equipped with an electron backscatter diffraction detector EBSD. The accelerating voltage is about 10 kV.

Ruthenium tetroxide RuO 4 was used as a staining agent to enhance the contrast of images of blends by increasing the electron density [START_REF] Sawyer | Polymer Microscopy[END_REF]. Samples were trimmed with a glass knife under nitrogen atmosphere below the glass transitions of the PE and PP at -125°C [START_REF] Trent | Ruthenium tetraoxide staining of polymers for electron microscopy[END_REF] to obtain small analysis surface using an ultramicrotome. Samples were stained during 3 to 4 h for LDPE/iPP blends and during 24 h for HDPE/sPP blends in the vapor phase. The staining solution was prepared by reacting hydrated ruthenium chloride (0.2 g RuCl.xH 2 O, Aldrich) with sodium hypochlorite (10 ml 5% active Cl, Aldrich). Then, the stained samples are washed with distilled water and dried for 48 h. After staining, samples become hard and this makes possible to cut the surface at room temperature to obtain a smooth surface using a diamond knife. Before observation, samples are coated with a thin layer of conductive carbon using a sputter coater. The stained phase appeared bright and the less stained phase appeared dark in SEM [START_REF] Brown | New method for the characterization of domain morphology of polymer blends using ruthenium tetroxide staining and low voltage scanning electron microscopy (LVSEM)[END_REF]. Finally, ImageJ software was used to measure the volume average radius of the dispersed morphology or the characteristic length l c for the co-continuous morphology.

The measurements of the characteristic length was made by drawing lines along the length and width of the SEM micrograph and by counting the number of intersections (N) between the two phases. The characteristic length (l c ) is then deduced from the line length divided by N/2. 200 particles or intersections are analyzed from SEM micrographs taken from different samples. A slight distortion of the co-continuous structure was sometimes observed because of the passage through the die. For this reason, the measurements of the l c dimension was made in two different orthogonal directions and averaged. The interfacial area was calculated using Equation ( 1) for dispersed phase morphology and Equation (3) for co-continuous structures. In the transition zone where both structures coexist, the calculation was performed by both equations defining the upper and lower limit of an error bar.

Shear rate in the twin-screw extruder

The Ludovic ® software (Sciences Computers Consultants) was used to simulate the flow in the co-rotating twin-screw extruder. It enables the calculation of pressure, residence times, mean temperature and shear rates along the screw during processing of a given material in given conditions. In this study, it was used to simulate the shear rate as a function of the screw profile and processing conditions. The calculated shear rate was an averaged value in the screw section where the polymer is molten. A weighting by local residence times was performed.

RESULTS AND DISCUSSIONS

Data used for the modeling

Figure 3 shows the viscosity at 190°C measured in oscillatory shear for the various polymers. given in Table 2. The interfacial tension at 190°C were measured by the method proposed by

Gramespacher and Meissner and found to be 0.7 mN/m between LDPE and iPP and 1.0 mN/m between HDPE and sPP. These results are in agreement with those found in the literature [START_REF] Samuel | Melt compatibility between polyolefins: Evaluation and reliability of interfacial/surface tensions obtained by various techniques[END_REF][START_REF] Veenstra | On the coarsening of co-continuous morphologies in polymer blends: effect of interfacial tension, viscosity and physical crosslinks[END_REF].

Experimental observation of the morphology

Figure 4 shows a superposition of the SEM pictures with the storage modulus at low frequency. The excess of elasticity in comparison with a simple blending law is due to the relaxation of shape of the phases in immiscible blends that induces long relaxation times [START_REF] Galloway | Comparison of methods for the detection of cocontinuity in poly (ethylene oxide)/polystyrene blends[END_REF].

The storage modulus shows a higher value in dispersed phase morphologies because of the large interfacial area for these morphologies. The occurrence of continuity of any phase is marked by a maximum of the storage modulus since the co-continuous morphologies lower the interfacial area. These maxima are observed at 30 ± 5 % and 70% ± 5 % of iPP for LDPE/iPP blend. For the HDPE/sPP blend, the maximum for low amount of sPP is not so clearly visible, probably because the frequency of 0.159 rad/s does not match with that of the shape relaxation in this case of very different viscoelastic behaviors between the two phases.

On the right side of the graph, when sPP is the major component, the maximum is observed at 80% ± 5 % of sPP. The comparison between the calculated and experimental interfacial area is given in Figure 5c. Large error bars for the experimental data are obtained in the intermediate zones due to the coexistence of the dispersed phase and co-continuous morphologies. The order of magnitude of the interfacial areas is generally captured by the model though slightly underestimated in comparison to that from SEM pictures particularly at 70% of iPP. Apart from any defect of the model itself, there might be some experimental causes for a finer structure, responsible for both the higher interfacial area and the narrowing of the co-continuous window. This cannot be attributed to the shear in the die because the shear rate in the die is lower than in the extruder (50 s -1 ) and the residence time is also short (< 5 s). Additional breakup of the filaments in the co-continuous phase at rest at the exit of the die could be a possible reason since the samples were cooled at room temperature.

Continuity diagrams, characteristic sizes and interfacial areas from the model

Figure 6 shows the same set of diagrams for the HDPE/sPP blend in the same processing conditions. For this blend, the diagram of continuity (Figure 6a) is strongly asymmetric because the viscosities of the pure components are very different. From the continuity 22 diagram, the change of morphology (d c =0.5), is observed at 25% and 54% of sPP. From 0 to 19%, sPP is dispersed in HDPE with an average diameter of 1.0 µm. From 35% to 40%, both polymers are continuous, above 63%, HDPE is dispersed in sPP with an average diameter of 1 µm. The co-continuous zone is narrower than for LDPE/iPP blend because the interfacial tension is higher. The phase inversion is calculated at slightly different compositions for the various dimensions (D, l c , a * ) because the dependence of these quantities with the composition is not the same. This is strongly visible when the viscosities of the components are very different. It was not observed in the case of the LDPE/iPP blend for which the viscosity ratio is close to one. Calculated and experimental interfacial area of Figure 6c shows that the order of magnitude and the variation with the composition is correct unless for the blend with equal amount of PE and PP. Except for this particular composition, the results of the calculation agree well with the storage modulus data of Figure 4b which showed that the maximum of the storage modulus on the left side of the diagram is not visible and this result is consistent with a weak maximum of the interfacial area.

Evidences for the capillary instabilities as the main breakup mechanism

The model proposed in this study assumes the capillary instabilities to be the main mechanism for the breakup of droplets of the dispersed phase. This was based on the high capillary number in comparison to the critical capillary number and on the high deformability of coarsened droplets following coalescence. Figure 7 compares the interfacial area for the LDPE/iPP blends calculated assuming a Taylor stepwise breakup and capillary breakup. The latter is clearly more suitable for modeling the experimental data, since the diameter is not very sensitive to the composition in the dispersed phase zone. Indeed, the theoretical equilibrium size by the stepwise mechanism is very small and leads to large interfacial areas that are not observed under practically relevant conditions. Actually, some SEM pictures taken at the edges of the continuity diagram for PE / PP blends show a series of align droplets with the same radius that can originate from a thread breakup suggesting the capillary instabilities of drawn filaments, though satellite droplets were never observed, probably because of coalescence.

Effect of changes of the flow conditions on the continuity diagram

Figure 8 shows the predicted effect of an increase of the shear rate from 70 s -1 to 200 s -1 on the interfacial area in the case of the LDPE/iPP blends. This graph contains all the relevant information and shows two major effects. Firstly, the co-continuous zone narrows, the lower limit increasing from 27% to 38% and the upper limit decreasing from 72% to 60%. Secondly, interfacial area is higher, meaning that the characteristic size of the morphology decreases in the co-continuous phase. Among possible improvements of the model, the effect of elasticity in conditions of large shear rates could be considered. Elasticity of the dispersed phase increases the resistance to deformation of droplets. It could be implemented easily in the model in the capillary number by subtracting the normal stress difference (calculated as two times the storage modulus) to the shear stress as proposed by Sundararaj et al. [START_REF] Sundararaj | Drop breakup and coalescence in polymer blends: the effects of concentration and compatibilization[END_REF]. Basically, the effect can be expected to be similar to that of decrease of the shear stress at constant viscosity and viscosity ratio. This results in an enlargement of the co-continuity window as shown on Figure 8. Practically, due to the shear shinning behavior of viscoelastic polymer melts, any increase of the shear rate may change both viscosity and elasticity of the phases simultaneously. Taking into account the viscosity and elasticity ratio in the model could help to predict and control the dispersedto-continuous transition (DCT) for blends with large asymmetric viscoelastic behavior as suggested by You et al. [START_REF] You | Control of the dispersed-to-continuous transition in polymer blends by viscoelastic asymmetry[END_REF].

CONCLUSION

Using some fundamental knowledge on the morphology development presented in the literature, a micro-rheological model that allows predicting the morphology of immiscible polymer blends was proposed. In particular, the model was based on a description of the topology of phases in both dispersed and co-continuous morphologies, which is non-trivial in the latter case. The model predicts the continuity, the characteristic sizes and the interfacial area diagrams for the entire composition range. It requires the knowledge of easily accessible input data, such as the mean shear rate in the mixing device, viscosity of each phase and viscosity ratio at this shear rate and the interfacial tension between the polymers. It was shown that for the investigated blends, morphology development is governed by capillary instabilities rather than Taylor break-up mechanism. The core of the model is a criterion of stability of elongated droplets forming threads of the minor phase that defines the type of morphology. The criterion relies on the comparison of the dimension of threads with a critical dimension that was taken to be 1 µm according to literature. The rational for this critical dimension is unfortunately not easy in simple shear flow and was not established. This remains a weakness of the analysis though the model is not very sensitive to the choice of this critical value: a 25% change only induced a 6% variation of the maximum interfacial area and a change of the limits of the continuity of 3 vol-% on each side of the continuity diagram. The model was tested by comparison of the prediction with experimental data on two polyolefin blends LDPE/iPP and HDPE/sPP processed by twin screw extrusion on the entire range of compositions. This opens the route to the prediction of functional properties of these blends, that are known to be strongly dependent on the morphology.

GLOSSARY

Figure 1 :

 1 Figure 1: Schematic of tree parallelepipeds with orthogonal faces inserted in a cubic mesh as

A

  model describing the evolution of the blend morphology was established on the basis of the balance between droplet or filament breakup and coalescence. The various events and paths of this evolution are summarized in Figure 2. The characteristic dimensions in this diagram are: very diluted droplets originating from Taylor's mechanism, a * the critical diameter of stability for filament break-up, (a) and l c the filament diameter and length of elongated filaments in the co-continuous structure.

Figure 2 :

 2 Figure 2 : Schematic flow diagram of the micro-rheological model: the different evolution

  matrix viscosity taken at the relaxation frequency,  d  dispersed phase viscosity taken at the relaxation frequency,  = volume fraction of the dispersed phase, = interfacial tension, R v = volume average radius of the dispersed phase from SEM pictures.

Figure 3 :

 3 Figure 3: Complex viscosity of the various polymers, measured in oscillatory shear flow at

Figure 4 :

 4 Superposition of the SEM pictures with the variation of the storage modulus at low frequency (190°C, 0.159 rad/s). a) LDPE/iPP, b) HDPE/sPP.

FigureFigure 5 :

 5 Figure 5a shows the diagram of continuity calculated for the LDPE/iPP blend at a shear rate of 70 s -1 . The corresponding characteristic sizes are shown on Figure 5b. For this blend, the

Figure 6 :

 6 Figure 6: Results of the calculation for the HDPE/sPP blend prepared at 190°C, 70s -1 . a)

Figure 7 :

 7 Figure 7: Comparison of the predicted interfacial area for capillary instabilities and Taylors's breakup mechanisms in the case of LDPE/iPP blend (190°C, 70s -1 ).

Figure 8 :

 8 Figure 8: Prediction of the effect of an increase of shear rate from 70 s -1 to 200 s -1 on the

Table 1 :

 1 Physico-chemical properties of the polyolefins (rad: radicalar, m: metallocene, ZN:

	including information on catalyst

Table 2 :

 2 Data of the blends for calculations (Viscosity at 190°C, 70s -1 ).

	Major phase	Minor phase	 m (Pa.s)	 d (Pa.s)	p	 (mN/m)
	LDPE	iPP	330	340	1.03	0.7 ± 0.1
	iPP	LDPE	340	330	0.97	0.7 ± 0.1
	HDPE	sPP	920	1630	1.77	1.0 ± 0.2

ACKNOWLEDGEMENTS

We wish to acknowledge Ms Stéphanie Reynaud, Research Engineer at the Hubert Curien Laboratory, UMR CNRS 5516, for her valuable help in SEM analysis.


Pierre Bochet

Francis Ramond

Renaud Touraine

Thierry Thomas

Hubert Marotte
email: hubert.marotte@chu-st-etienne.fr

Arthritis associated to cardio-facio-cutaneous syndrome related to a MAP2K1 mutation

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.

2

Image of a 17 years old female patient with cardio-facio-cutaneous (CFC) syndrome (Figure 1), who developed since six months a symmetric polyarthritis affecting hands (interphalangeal without metacarpophalangeal involvement), wrists, elbows, knees, and ankles without night awakening or morning stiffness associated with biological inflammation. No immunological abnormality or erosion (radiographs or ultrasonography) was observed. Clinical and biological remissions were obtained after three months of methotrexate at 10 mg/week in subcutaneous.

CFC is a rare genetic disease where activating mutations of the RAS-MAPK genes (Ras/mitogen activated protein kinase) [START_REF] Nava | Cardio-facio-cutaneous and Noonan syndromes due to mutations in the RAS/MAPK signalling pathway: genotype-phenotype relationships and overlap with Costello syndrome[END_REF]. Here, the CFC was confirmed with a heterozygous c.275T>G variant of MAP2K1 gene.


Laurent Chircop-Reyes

What It Takes To Be A Caravaneer?

Keywords: Human Mobility, Network, Security, Migration, Caravan Knowledge, China

The study of historical mobility invites reflection on adaptation to environmental and social change. The different ways in which humans interacted with their environment in the context of uncertainty, economic and social crises or hostile conditions, as well as the multiple ways in which various communities have practised space by moving. There are many excellent works on human mobility, both in Asia and in other parts of the world (for Asia see for example Guerassimoff et al., 2020;[START_REF] Jami | Mobilité humaine et circulation des savoirs techniques (XVII e -XIX e siècles)[END_REF]Li, 2010;[START_REF] Yi | Transforming Inner Mongolia : Commerce, Migration and Colonization on the Qing Frontier[END_REF] or in other regions see [START_REF] Noiriel | L'Immigration : naissance d'un 'problème[END_REF][START_REF] Ahmida | Bridges Across the Sahara: Social, Economic, and Cultural Impact of the Trans-Sahara Trade During the 19th and the 20th Centuries[END_REF]. They are too numerous to mention here. However, it seems to me that many questions remain unanswered concerning the plural functioning of Chinese human mobility across the land in the North (see fig. 3), and the different local knowledge used to move by caravaneers and migrant communities before and at the turn of the 20th century. In this article, we attempt to understand the plural functioning of circulations, particularly in spaces where the trajectories of different groups were deployed in hostile environments and on the periphery of the scope of political action. In particular, we discuss what it took to move through these regions before the great cycles of modernisation at the beginning of the 20th century.

Constrained and Risky Travel

The cost of crossing the Loess Plateau and the Sino-Mongolian steppe was sometimes high. Before the first half of the 20th century, these remote areas were indeed risky, for natural and human reasons. Windy and icy in winter, subject to intense drought in summer, these regions were also characterised by a relative absence of government. Various forms of brigandage were practised on caravans travelling north and west, especially from Shahukou (Northern Shanxi) (see figs. 5 and 7), also known as the "gateway to the West" (Xikou 西口). The abuses were particularly violent against groups moving with high-value goods, such as merchants from Shanxi Province (Jinshang 晉商), some of whose families went to trade far away and migrated as far as Kiakhta, a "trade town" (maimaicheng 買賣城) in Russia. But also migrants with hopes of a better life in the vastness of the steppe "beyond the fortifications" (saiwai 塞外). Coming from diverse professions (farmers, artisans and also traders), they were forced to leave Shanxi and the surrounding provinces (Shaanxi, Hebei) from the 15th century onwards, among other reasons because of environmental problems (drought, aridity of the soil) and overpopulation.

Merchants and migrants struggled to provide their own security. However, the powerful and corporately organised merchants could still seek the private protection services of escortcaravaneers (baobiao 保鏢). Between the eighteenth and early twentieth centuries, escortmasters (biaoshi 鏢師), as interstitial actors, bridged the gap between merchants and brigands, on the one hand, and participated in the maintenance of trans-regional trade as well as in the marking of trails for commercial and migratory journeys, on the other. Some previous work (including my own) on the phenomenon of escort companies (biaoju 鏢局) has enabled us to learn a great deal about the role of knowledge for moving in a context of insecurity. These studies concern especially communicative knowledge (verbal and non-verbal codes, jargon, vehicular language) and defensive knowledge (martial arts practices).

However, these escorted human circulation are only the tip of the iceberg that was Chinese human mobility in the North.

Walking to the West

Trade routes and migrations were intensely connected over several generations in North China. Some have shaped the phenomenon of "zou Xikou 走西口" (walking [towards] the West gate). This term is derived from local folklore. It refers to the local songs, dances, and operas still practised today by the people of Shanxi, whose ancestors once migrated to the northwest. These stories tell, among other things, of the nostalgia for the native land, but also of the arduous and dangerous nature of the journeys during which the caravans were regularly victims of brigandage. Zou Xikou thus has a strong local historical and cultural dimension. Its presentday heritage is a valuable source of field data. This mobility also participated to shape the sino-russian "tea trade routes" (chaye shanglu 茶 葉商路) in the early 18th century. The settlement of Shanxi people along these routes was supposed to be temporary. Partly for reasons related to respect for ancestor cults. But Cui [START_REF] Sipeng | Qing Daoguang shiqi Shanxiren 'zou Xikou' cheng yin tanxi » 清道光時期山西人 '走西口' 成因探析 [The causes of the 'westward walk' of the Shanxi people during the Daoguang reign, Qing Dynasty[END_REF] point out to us that, in order not to "dishonour the families left behind", only those groups who were "successful in trade could harbour the hope of returning to their homeland".

No study to date can indicate a precise number of Shanxi people who left the province to trade and migrate to the north-west. In comparison, more accurate data are available for the northeastern migration, known as the "Manchuria rush" (chuang Guandong 闖關東) -lit. "East gate rush". A report in 1933 (Zhang, 2001: 322) gives the (underestimated) figure of 170,000 having migrated to the provinces of Jilin, Liaoning, Heilongjiang, as well as to the Siberian and current Inner Mongolia regions. Representatives of Shanxi merchant companies, who migrated to Russia and then returned to their province following the October Revolution in 1917, estimated that at the beginning of the 20th century more than 10,000 Shanxi merchants were living in Russia. Migration and trade were, according to historian Zhang [START_REF] Zhengming | Jinshang xingshuai shi 晉商興衰史 [History of the rise and fall of Shanxi Merchants[END_REF], the two most common activities in the daily lives of Shanxi people.

Desert Ship and Professional Escort

What did long-distance caravan trade consist of? First, animals were indispensable. These were animals exploited for their carrying power (camels, mules and horses), but also for security (guard dogs) or communication (domesticated pigeons). While the mule accompanied (and still does) the caravans from Yunnan to the Tibetan, Burmese and Vietnamese regions (Ma Jianxiong, 2016), the camel, a true "desert ship", renowned for its endurance and resistance, played a central role in these expeditions. The caravan size varied between twenty and eighty individuals, and a hundred animals per expedition. This is of course incomparable with the colossal number of sub-Saharan caravan families between the 8th and 16th centuries, who could lead caravans of thousands of dromedaries to protect themselves from brigand attacks (Bernard Nantet, 2008). In order to overcome the problem of brigandage, the Shanxi Merchants entrusted their goods to private armed militias, some of which maintained tacit "friendly" relationships with the brigands to "rent the passage" (jielu 借路) of the caravans. Note that this passageway trade worked in two ways: the caravaneers bought their way from the brigands, and the brigands could also buy their traffic from the regular soldiers posted at the interregional border posts (guankou 關口). The Shanxi caravans (also those of Hebei and Shandong) were thus unique in that they were composed of professional escorts (biaohang 鏢行). That is to say, specialists mastering, among other things, human power transportation techniques (driving handcarts), brigands' jargon and local languages (such as 'dark language', heihua 黑話), as well as martial arts techniques and the use of the rifle. Both the escorts and their companies belonged to family or clan lineages of martial arts schools, and recruitment into a company was done internally and according to the rules of the particular lineage. The members of an escort company had a different status depending on their seniority and skills, and on the knowledge they possessed by natural disposition or transmission. For example, there are those who are recruited for their knowledge of the environment and ways of communicating with different local groups, or for their mastery of a particular skill (healer, craftsman). There are also those who are recruited for their good eyesight, others for their keen hearing, and still others for their olfactory sensitivity -being able to recognise which animal a particular smell belongs to (via faeces) in order to detect the possible presence of brigands, who usually travel with horses. Long-distance travel required a keen sense of direction and a fine "reading" of the environment. These parameters were vital to avoid getting lost and to be able to locate water points. An expedition could take between one and three months from Shanxi to Hebei or Inner Mongolia, or even six months to the Russian border. The caravaneers were also affected by the deleterious natural conditions. A simple badly healed chapped skin or a superficial lesion of the skin could cause cracks, followed by a skin infection, fever, etc., and could compromise the progress of the whole group. Among the many "prescriptions" to be followed was that of not washing during the entire journey. This protected them from the dry, cold climate by allowing the fat secreted by the pores of the skin to accumulate and cover the face until it formed an effective protection against the wind and sand. Moving through these spaces meant mastering the techniques of fire production, cave dwelling, etc. For fear of poisoning, the master escorts were strictly forbidden to eat food from the inns. The latter, which resembled modest caravanserais, were rest stops on the road, but not places of trust. According to local oral tradition, it was not uncommon for innkeepers to be associated with local brigands

The studies carried out on Chinese human mobility across the land in the North (see figs. 2 and 3) prior to the 20th century indicate that caravaneers could not do without individuals who mastered local knowledge and the most elementary techniques of daily life. They highlight the relevance of investigating a hitherto little-investigated everyday life, whose apparent "banality" reveals situations that are rich in information about know-how (knowing how to repair and make shoes, for example) and social functioning.

  figure 2: Provinces in Qing China and Great Wall after 1550. Source: William T. Rowe, China's Last Empire. The Great Qing (2009).

  figure 4 : A caravan of goods on camels between Shanxi and Hebei provinces. Source: T.C Chamberlin, 1909.

figure 7 :

 7 figure 7: An ancient route in Shahukou 杀虎口, Shanxi/Inner Mongolia "border". Source: photographs by Laurent Chircop-Reyes, June 2018.

  What It Takes To Be A Caravaneer? Laurent Chircop-Reyes 2021 7 GIS Asie -French Academic Network on Asian Studies Articles of the month figure 5: Routes recognised and taken by escort companies during the Qing period. Source: Liu Yinghai, Qiao Zengguang, Biaoxing sihai 鏢行四海 (2014).

	What It Takes To Be A Caravaneer?
	Laurent Chircop-Reyes 2021
	figure 6: A handcart passing through a loessic valley in northern Shanxi.
	Source: Yadong yinhua ji, 1924.
	8
	GIS Asie -French Academic Network on Asian Studies
	Articles of the month


Zicheng Liu
email: zicheng.liu@telecom-paristech.fr

Changyou Li

Yu Zhong

Dominique Lesselier

Electromagnetic modeling of damaged fiber-reinforced laminates

Keywords: nondestructive testing, array scanning method, multipole method, damaged periodic structure, fibered laminates, Green's function, equivalence theory, plane-wave expansion, cylindrical-wave expansion

As a prerequisite to nondestructive testing of damaged fibered laminates, the Green's function, corresponding with an undamaged structure, and the electromagnetic fields associated with the damaged one are investigated herein. For the undamaged fibered laminate, benefiting from the periodicity of the fibers within each layer, the field solution follows the scattering-matrix-based method using the Floquet theorem. Yet, the periodicity is destroyed by the analytical source (for the Green's function) or by damages, and the Floquet theorem cannot be directly applied to compute the associated scattering matrices. The array scanning method is introduced to that effect. Inserting fictitious sources to get a quasi-periodic source array, the modeling approach for undamaged laminates can be used to compute the field with the source array, the integration of which cancels the effects of the fictitious sources and yields the Green's function. With the multipole method, field disturbances by damages, which include missing, displaced, shrunk, and expanded fibers and circular inclusions inside fibers, are accurately modeled by setting equivalent sources inside sound fibers, and the array scanning method applies. Modeling accuracy and efficiency of the approaches are illustrated by numerical simulations.

Introduction

Fiber-reinforced laminates are widely used in aerospace and automotive industries due to advantages in stiffness and weight. Each layer of the laminate is composed of a planar slab and an embedded periodic fiber array. Fiber material might be glass or graphite, and the slab itself is usually made of epoxy or polyester. Stacking the reinforced slabs yields the concerned structure and, to provide strength along various directions, fiber orientations in different layers usually differ. Nondestructive testing of potential damages inside is necessary in order to ensure reliability and safety. Electromagnetic (EM) means can be used to that effect because of broad applicability on both conductive [START_REF] Heuer | Review on quality assurance along the CFRP value chain -Non-destructive testing of fabrics, preforms and CFRP by HF radio wave techniques[END_REF] and dielectric [START_REF] Zoughi | Microwave Non-Destructive Testing and Evaluation Principles[END_REF] composites. Correspondingly, accurate modeling and high-resolution imaging methods are requested to guarantee the detection performances, the former rarely reported in literature and being focused onto in the present contribution.

As indicated, modeling of damaged and undamaged laminates is required to better study impacts of damages and configure the testing systems, but may not be easily handled by commercial softwares (e.g., COMSOL) as the boundary conditions are not clear due to the periodicity at different directions in each layer. Efforts have been made by (some of) present authors. Without damages, the modeling approach refers to works in [START_REF] Li | Full-wave computational model of electromagnetic scattering by arbitrarily rotated 1-D periodic multilayer structure[END_REF] where fields near slab boundaries are expanded into plane-wave expansions (PWEs) in order to facilitate modes matching, and the multipole method [START_REF] Yasumoto | Electromagnetic Theory and Applications for Photonic Crystals[END_REF] is adopted to analyze the fiber array. However, the applicability of such methods to the disorganized structure is prevented by the fact that the periodicities are destroyed. Approximated field solutions can be provided by brute-force codes running FDTD [START_REF] Wang | Efficiency-improved LOD-FDTD method for solving periodic structures at oblique incidence[END_REF] or FEM [START_REF] Tyrus | A local finite element implementation for imposing periodic boundary conditions on composite micromechanical models[END_REF] methods. Since regions near the damages are those of main interest, a truncated model (a cell) of the laminate centered on the damages is often investigated. To guarantee accuracy, testing the field convergence is needed by enlarging the cell, i.e., increasing the truncation size. A similar methodology is followed in our previous work [START_REF] Liu | Electromagnetic modeling of damaged single-layer fiber-reinforced laminates[END_REF], save periodic conditions being set at the end sides of the cell to fabricate a new, fictitious periodic structure; then, quasi-periodic theories apply. But the unknown approximation errors (due to lack of the true solution) and possibly heavy computational cost in case of low-loss materials show the necessity of developing a more sophisticate modeling.

Contributions have been made by authors in [START_REF] Liu | Fast full-wave analysis of damaged periodic fiber-reinforced laminates[END_REF] on efficiently and accurately modeling damaged single-layer laminates with missing, displaced, shrunk, and expanded fibers and circular inclusions inside fibers. Defects are electromagnetically equivalent to sources within the initially organized fibers by reformulating the changed fiber boundary conditions or via application of Graf's addition theorem. Then the field solution is a summation of responses due to the exterior emitting source and to the equivalent ones. The former is obtained with the method dealing with undamaged structures [START_REF] Li | Full-wave computational model of electromagnetic scattering by arbitrarily rotated 1-D periodic multilayer structure[END_REF]. The latter is obtained via the array scanning method (ASM) [START_REF] Capolino | Fundamental properties of the field at the interface between air and a periodic artificial material excited by a line source[END_REF] where a phased source array is assumed so that quasi-periodic theories become applicable and the final solution results from integration over the first Brillouin zone.

In this contribution, the equivalence theory is applied on the modeling of damaged multilayer laminates. Treating the concerned damages as equivalent sources, the solution to the illumination of undamaged multilayer laminates by exterior and equivalent sources is required and obtained with the help of scattering matrices [START_REF] Li | Recursive matrix schemes for composite laminates under plane-wave and Gaussian beam illumination[END_REF], which describe mutual scattering between two layers. Playing the key role in both modeling and imaging, the solution to the Green's function, i.e., the response to a line source in 2-D cases, is also explored following the same computational approach as for equivalent sources.

Representing the field in the whole region by PWEs or CWEs (short for cylindrical-wave expansions), the above solution is achieved through the retrieval of expansion coefficients, for which PWEs and CWEs are truncated and the truncation size should be properly set. A numerical strategy, which makes use of the knowledge that PWEs and CWEs should yield the same field solution in the commonly supported region, is proposed, while the possible different structures (in terms of material and of geometry) of layers are considered to yield the correspondingly different truncation sizes.

The contribution is structured as follows. Notations and configurations are introduced in Section 2, the undamaged structure being of first concern. The exterior illumination is studied in Section 3. The solution to the Green's function is presented in Section 4. With quantities obtained above, different types of damages are treated as equivalent sources in Section 5 and how to reach the true solution is discussed as well. The truncation number of PWEs and CWEs is the key parameter in the modeling. A numerical strategy is presented in Section 6 to adaptively select the optimal value. The numerical analysis is considered in Section 7 and conclusions are proposed in Section 8. Appendices provide complementary materials. The description of undamaged laminates goes first, and then the one of concerned damages. The structure of the multilayer laminate is sketched in Fig. 1, where it is sandwiched by two half spaces at top and bottom, each layer of the laminate being made by periodically inserting similarly-orientated fibers inside a planar slab. The composition of different layers can differ w.r.t. fiber radius, slab thickness and materials. The fiber periods are to be the same here.

Configurations and formulation

Assume L layers are stacked on one another. Labeling them with integers from 1 to L, the upper and lower half spaces are indexed by 0 and L + 1, respectively. Setting the origin of the global coordinate system at the center of a fiber in the 1st layer, the local coordinate system for the l-th layer is defined by setting the origin at the fiber center o (l,0) with global coordinates (o

(l,0) x , o (l,0) z ), t l+1 < o (l,0) z < t l , -d/2 < o (l,0) x ≤ d/2, l = 1, . . . , L, the plane z = t l
denoting the lower boundary of the l-th layer and d the fiber period. Then, the center of fibers in the l-th layer has local coordinates (x l , z l ) = (vd, 0), v ∈ Z, based on which fibers in the l-th layer are labeled by integer couples (l, v).

Remark that the subscript of x and z indicates the coordinate system, i.e., "l" stands for the local coordinate system of the l-th layer.

All materials, including the one in the inclusions, are linear, isotropic and homogeneous. Fiber and slab materials are layer-dependent, characterized by relative permittivities ε l j (complex value for lossy materials) and relative permeabilities µ l j , j = " f " for fiber and "s" for slab. Consider now the damaged laminate. Indexes of damaged fibers compose the set D, i.e., the v-th fiber in the l-th layer is damaged if (l, v) ∈ D. A tilde is put above the corresponding notations to denote the damaged status. One respectively denotes the center and radius of the (l, v)-th fiber by o (l,v) and c l for the undamaged laminate, the local coordinates of õ(l,v) = (vd, 0) for displaced fibers and c(l,v) = c l for shrunk/expanded ones. A circular inclusion is configured by its center ō(l,v) and radius c(l,v) , with a bar above corresponding notations. Without loss of generality, half spaces and inclusions are assumed filled by air and missing fibers are "fibers" sharing the same material with the background slab.

With time-dependence exp(-iωt), i imaginary unit, ω angular frequency, the incident waves considered are either time-harmonic plane waves of (x, z) plane of incidence or line sources orientated along the y axis, with either TM or TE field polarization. This configuration leads to a two-dimensional scalar scattering problem, which means only the y-component of electric field (for TM) or magnetic field (for TE) needs to be tackled.

Electromagnetic modeling of undamaged laminates

The modeling of damaged laminates is based on the analysis of the corresponding undamaged ones. For completeness, this Section revisits the modeling approach [START_REF] Li | Recursive matrix schemes for composite laminates under plane-wave and Gaussian beam illumination[END_REF] for undamaged laminates illuminated by a planar wave (considering that other types of sources might be analyzed by exploiting their PWEs [START_REF] Yang | Two-dimensional scattering of a Gaussian beam by a periodic array of circular cylinders[END_REF][START_REF] Watanabe | Accurate analysis of electromagnetic scattering from periodic circular cylinder array with defects[END_REF]).

Analyzing wave propagations within each layer through the combination of multipole method and field matching on boundaries, scattering behaviors in different layers are linked by scattering matrices (rather than possibly illconditioned transfer matrices [START_REF] Li | Recursive matrix schemes for composite laminates under plane-wave and Gaussian beam illumination[END_REF]), which are computed according to a recursive scheme. As shown in Fig. 2, for the l-th layer, the fiber array divides the layer into three slices, region above the array, region below it, and the array itself, denoted by Γ l+ , Γ l-and Γ l f , respectively. Due to the periodicity of the structure, the Floquet theorem applies and the total field follows a pseudo-periodic distribution, i.e., V (r + d x) = V (r)e iα 0 d , where V stands for the y-component of the E-or H-field depending upon the wave polarization, r denotes the observation position, x is the unit vector along the x axis, and α 0 = k 0 s sin θ inc , k 0 s being the wavenumber of the upper half space and θ inc being the incident angle of the plane wave measured counter-clockwise from the z axis.

The pseudo-periodic property of V yields its representation in PWEs [START_REF] Botten | Formulation for electromagnetic scattering and propagation through grating stacks of metallic and dielectric cylinders for photonic crystal calculations. Part I. Method[END_REF]. For the region of Γ l+ and Γ l-, the expansion is with coefficients f l+ p , g l+ p and f l- p , g l- p , respectively. The subscript "p" denotes the expansion order and "g", "f" stand for field components propagating along the z-axis and the opposite direction, respectively. Thus, in the local coordinate system,

V (r) = ∑ p∈Z f l± p e -iβ l p z l + g l± p e iβ l p z l e iα p x l , r ∈ Γ l± , (1)
where

α p = α 0 + p2π/d, β l p = (k l s ) 2 -α 2 p .
Matching fields at the slab boundary z = t l , expansion coefficients at the two sides of this boundary are related as

f l+ g (l-1)-= E 1 l E 2 l E 3 l E 4 l f (l-1)- g l+ . ( 2
)
where

f l± = [ f l± p ], g l± = [g l± p ]
, and the definition of E u l , u = 1, 2, 3, 4, is given in Appendix A. Since f l+ , g (l-1)- stand for waves scattered by the boundary and f (l-1)-, g l+ for incident ones, the matrix composed by E u l indicates the scattering property of the slab boundary at z = t l . Such relation w.r.t. the fiber array is also obtained with the multipole method [START_REF] Yasumoto | Electromagnetic Theory and Applications for Photonic Crystals[END_REF]. Expressing the total field in the vicinity of and inside the (l, v)-th fiber in cylindrical-wave expansions (CWEs),

V (r) = ∑ m∈Z a (l,v) m ϕ J m k l s , r (l,v) + b (l,v) m ϕ H (1) m k l s , r (l,v) , (3a)
V (r) = ∑ m∈Z c (l,v) m ϕ J m k l j , r (l,v) + q (l,v) m ϕ H (1) m k l j , r (l,v) , (3b)
where r (l,v) is the relative position of the observation point r w.r.t. the fiber center o (l,v) , i.e., r (l,v) = r-o (l,v) . Denoting

(r (l,v) , θ (l,v) ) as the polar coordinates of r (l,v) , ϕ Z m (k l j , r (l,v) ) = Z m (k j , r (l,v) ) exp(imθ (l,v) )
, where Z = J stands for the Bessel function and Z = H (1) for the first-kind Hankel function. The coefficients of the CWEs are sketched in Fig. 3.

While the coefficient of c

m e iα 0 vd ϕ H (1) m (k l s , r (l,0) -vd x) = ∑ p∈Z ∑ m∈Z K l± m,p b (l,0) m e i(α p x l ±β l p z l ) , (4)
where K l± m,p = 2(-i) m exp ±imθ l p /(dβ l p ), cos θ l p = α p /k l s , and the Floquet theorem has been applied with b

(l,v) m = b (l,0)
m e iα 0 vd . Thus, the total field in the l-th layer is expressed as a PWE as

V (r) = ∑ p∈Z e iα p x l f l+ p e -iβ l p z l + g l- p e iβ l p z l + ∑ m∈Z K l± m,p b (l,0) m e ±iβ l p z l , (5)
the matrix form of which is

V (r) = Ξ l-f l+ + Ξ l+ g l-+ Ξ l± K l± b (l,0) , (6)
where Ξ l± = [e i(α p x l ±β l p z l ) ], K l± = [K l± m,p ]. Signs + andare applicable in the region of z l ≥ 0 and z l < 0, respectively, for Ξ l± and K l± . The relationship between b (l,0) m , f l+ p , g l- p is exploited by analyzing the scattering behaviors of the fiber array.

With the field representation as CWEs in Eq. ( 3), the fiber boundary conditions yield the linear relations of CWE coefficients,

b (l,v) m = R l m a (l,v) m + T l m q (l,v) m , (7a)
c (l,v) m = T l, m a (l,v) m + R l, m q (l,v) m , (7b)
where R l m , R l, m , T l m , T l, m are reflection or transmission coefficients defined in Appendix A. Remark that q (l,v) m equals 0 since no inner source exists. Together with the so-called Rayleigh identities [START_REF] Wilcox | Modeling of defect modes in photonic crystals using the fictitious source superposition method[END_REF],

a (l,0) m = ∑ n∈Z S l m-n b (l,0) n + ∑ p∈Z J l+ m,p g l-+ J l- m,p f l+ , (8)
which is obtained with the Jacobi-Anger expansion [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] and application of the Graf's addition theorem [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF], the expression of b (l,0) m as a function of f l+ p and g l- p is written as

b (l,0) = I -R l S l -1 R l J l+ g l-+ J l-f l+ , ( 9
)
where

I is the identity matrix, b (l,0) = [b (l,0) m ], R l = diag R l m , S l = [S l m-n ], J l± = [J l± m,p ], J l± m,p = (i) m exp(∓imθ l p
) and S l m-n is the lattice sum defined as

S l m = ∑ v∈Z\{0} H (1) m (k l s |v| d)e i(α 0 ld-m arg(v)) , (10)
which enjoys a fast convergence about v for lossy slabs whereas the computation approaches in [START_REF] Yasumoto | Efficient calculation of lattice sums for free-space periodic Green's function[END_REF][START_REF] Denlinger | A fast summation method for oscillatory lattice sums[END_REF] apply for lossless cases.

The substitution of Eq. ( 9) into Eq. ( 6) leads to V as a function of PWE coefficients f l+ and g l-. Together with Eq. ( 1), the scattering matrix of the fiber array is obtained and described as

f l- g l+ = F 1 l F 2 l F 3 l F 4 l f l+ g l-, (11)
where

F 1 l = I + Ω l-R l J l-, F 2 l = Ω l-R l J l+ , F 3 l = Ω l+ R l J l-, F 4 l = I + Ω l+ R l J l+ , Ω l± = K l± I -R l S l -1 .
Then, the scattering effects of the fiber array are equivalent to a slab boundary with the scattering matrix defined by F u l , u = 1, 2, 3, 4. Representing the equivalent boundary by the slashed line, the sketch in Fig. 2 is reformulated as Fig. 4. Following the recursive scheme in [START_REF] Li | Recursive matrix schemes for composite laminates under plane-wave and Gaussian beam illumination[END_REF], scattering matrices linking PWE coefficients in any two slices are computed. Based on the computed scattering matrices, the PWE and CWE coefficients can be solved for.

The scattering matrix linking the slice Γ 0 and Γ l± is defined by

f l± g 0 = A l± 0 B l± 0 C l± 0 D l± 0 f 0 g l± , (12)
and the one linking Γ l± and Γ L+1 by

f L+1 g l± = A L+1 l± B L+1 l± C L+1 l± D L+1 l± f l± g L+1 . (13)
Considering the coefficient vector f 0 (PWE coefficients of the incident wave) is known and g L+1 = 0 due to no incident field from the bottom, one has equations

g l± = C L+1 l± f l± , (14a)
f l± = A l± 0 f 0 + B l± 0 g l± , (14b)
which yield the solution of f l± and g l± as

f l± = I -B l± 0 C L+1 l± -1 A l± 0 f 0 , (15a)
g l± = C L+1 l± f l± . ( 15b
) When l = L + 1, one has f L+1 = A L+1 0 f 0 and g 0 = C L+1 0 f 0 .
CWE coefficients are computed by substituting f l± and g l± into Eq. ( 8) and [START_REF] Capolino | Fundamental properties of the field at the interface between air and a periodic artificial material excited by a line source[END_REF]. With PWE and CWE coefficients, the full-field solution is achieved.

The Green's function

The Green's function, denoted by G(r, r s ), is defined as the field solution at r when the line source locates at r s and emits field H

(1) 0 (k l j r s )/(4i), j being "s" or " f " depending on the source location. Here, due to the role in the modeling of damaged laminates in Section 5, only the case of the line source inside the fiber is considered. Remark that the computation of G(r, r s ) w.r.t. the single-layer laminate has been presented in the authors' earlier works [START_REF] Liu | Electromagnetic modeling of damaged single-layer fiber-reinforced laminates[END_REF].

This Section gives the similar approach but with the combination of the scattering matrix to reach the solution of G(r, r s ) for the multilayer laminate.

Since G(rvd x, r svd x) = G(r, r s ), without losing the generality, the line source is assumed located inside the (l, 0)-th fiber. Section 3 concerns the exterior illumination on the undamaged structure. Thus, the presented modeling approach in Section 3 cannot be applied directly since the source is inside the laminate and destroys the quasi-periodic property of the field distribution. The array scanning method (ASM) overcomes this difficulty by fabricating a phased source array, as shown in Fig. 5. All sources share the same current except the one inside the (l, v)-th fiber is with phase exp (iα 0 vd), i.e., the true source is with phase 1. With the phased array, the field distribution becomes quasi-periodic and scattering due to the phased array can then be analyzed. Denote the field solution w.r.t. the source array by G(r, r s , α 0 , ∞). Following the identity d 2π

π/d -π/d e iα 0 vd dα 0 = δ v,0 , (16)
where δ v,0 equals 1 when v = 0 and 0 otherwise, the impacts of fictitious sources are expected to be canceled after integration over the first Brillouin zone, i.e., the true field solution is

G(r, r s ) = d 2π π/d -π/d G(r, r s , α 0 , ∞)dα 0 . (17)
Numerical methods [START_REF] Gonnet | Increasing the reliability of adaptive quadrature using explicit interpolants[END_REF] are used to perform this integration, while the computation of the integrand is introduced next.

Due to the presence of the source array, q (l,v) m in Eq. ( 7) is non-zero. Since the wave radiated by the line source is with expression H

(1)

0 (k l f r s )/(4i), r s = r s -o (l,0) , one has q (l,0) m = J m (k l f r s -o (l,0
) ) exp -im arg r s -o (l,0) /4i with the application of the Graf's addition theorem. Thus, Eq. ( 11) is updated as

f l- g l+ = F 1 l F 2 l F 3 l F 4 l f l+ g l-+ Ω l-T l q (l,0) Ω l+ T l q (l,0) , (18)
where q (l,0) = [q (l,0) m ], T l = diag T l m . Together with two equations, Eq. (14a) and (14b) (f 0 = 0 since no exterior source exists), the solution to f l± and g l± is obtained as

g l+ = Θ l+ C L+1 l-Ψ -1 l Ω l-+ Ω l+ Ξ l -Θ l+ C L+1 l-Ψ -1 l Θ l-B l+ 0 T l q (l,0) , f l-= Θ l-B l+ 0 Ξ -1 l Ω l+ + Ω l- Ψ l -Θ l-B l+ 0 Ξ -1 l Θ l+ C L+1 l-
T l q (l,0) ,

f l+ = B l+ 0 g l+ , g l-= C L+1 l-f l-, (19)
where

Θ l± = I + Ω l± R l J ± , Ξ l = I -Ω l+ R l J l-B l+ 0 , Ψ l = I -Ω l-R l J l+ C L+1 l-.
With the solution to f l± and g l± , as shown in Fig. 4, the laminate can be treated as two parts, layers below and above the l-th layer. The field solution for the former part is achieved by treating the field component with PWE coefficients f l- p as the incident wave and following the computational approach in Section 3. The field solution for the later part is similarly obtained with PWE coefficients g l+ p .

Damages modeled as equivalent sources

Through the analysis of the scattering behaviors, the authors found that the field disturbance due to damages, which as already indicated include missing, displaced, expanded, and shrunk fibers and circular inclusions, is equivalent to setting line sources inside the initially organized fibers, on condition that the CWE coefficients of the emitted fields satisfy specific relations. A brief introduction of equivalence theory in [START_REF] Liu | Fast full-wave analysis of damaged periodic fiber-reinforced laminates[END_REF] is given here for completeness.

Fig. 6. CWE coefficients inside the fiber with an circular inclusion.

For a circular inclusion, express fields around the inclusion with CWE coefficients shown in Fig. 6. The application of the Graf's addition theorem [START_REF] Liu | Electromagnetic modeling of damaged single-layer fiber-reinforced laminates[END_REF] leads to

q (l,v) m = ∑ n∈Z ϕ J n-m (k l f , o (l,v) -ō(l,v) ) b(l,v) n , (20a) ā
(l,v) m = ∑ n∈Z ϕ J n-m (k l f , ō(l,v) -o (l,v) )c (l,v) n . (20b)
Considering the linear relations of CWE coefficients in Eq. ( 7) and the following one from the boundary condition of

the inclusion, b(l,v) m = R(l,v) m ā(l,v) m , (21)
where R(l,v)

m is similarly defined as R (l,v) m , q (l,v) m
can be expressed as a function of a

(l,v) m , q (l,v) = Π (l,v) a (l,v) , (22)
where

a (l,v) = [a (l,v) m ], Π (l,v) = I -Φ (l,v) R l, -1 Φ (l,v) T l, , R l, = diag R l, m , T l, = diag T l, m , Φ (l,v) = [Φ (l,v) m,n ], and
Φ (l,v) m,n = ∑ u∈Z R(l,v) u ϕ J u-m (k l f , o (l,v) -ō(l,v) )ϕ J n-u (k l f , ō(l,v) -o (l,v) ). (23)
From the above analysis, one can conclude that a circular inclusion can be equivalently modeled as a line source at the fiber center on condition that its CWE coefficient q (l,v) m satisfies Eq. ( 22).

Like equivalences are valid for the other concerned damages. For the displaced fiber, following an analysis similar with the one for the circular inclusion, the relation in Eq. ( 22) is obtained but with a different definition of Π (l,v) ,

Π (l,v) = (T l ) -1 ϕ (l,v) 1 R l ϕ (l,v) 2 -R l , (24)
where ϕ

(l,v) 1 = [ϕ J n-m (k l s , o (l,v) -õ(l,v) ], ϕ (l,v) 2 = [ϕ J n-m (k l s , õ(l,v) -o (l,v
) )], õ(l,v) being the center of the displaced fiber. For the missing, expanded, shrunk fiber, one gets the corresponding definition of Π (l,v) by reformulating the linear relation of Eq. (7a). Since no fields are scattered by the missing fiber, b (l,v) m = 0. Then Eq. (7a) yields

q (l,v) m = - R l m T l m a (l,v) m . ( 25
)
For the expanded and shrunk fiber, the CWE coefficients are related by

b (l,v) m = R(l,v) m a (l,v) m , (26)
where R(l,v) m is similarly defined with R l m but with a different fiber radius. The rewritten form of Eq. ( 26),

b (l,v) m = R l m a (l,v) m + T l m 1 T l m R(l,v) m -R l m a (l,v) m , (27)
leads to

q (l,v) m = 1 T l m R(l,v) m -R l m a (l,v) m . (28)
Thus, Π (l,v) is defined as diag -R l m /T l m for the missing fiber and diag R(l,v) m -R l m /T l m for the expanded and shrunk fiber. The coefficient of q (l,v) m (or a (l,v) m ) needs to be solved to determine the equivalent source and get the field solution.

According to the scattering linearity about multiple sources, the field solution is decomposed as

V (r) = V inc (r) + ∑ (l,v)∈D ∑ n∈Z q (l,v) n G n (r, o (l,v) ). ( 29
)
where V inc is the field solution w.r.t. the exterior source and G n (r, o (l,v) ) as the field response due to the cylindrical wave ϕ H (1) n (k l j , r (l,v) ), i.e., the field solution by setting q (l,v) m = δ m,n in Section 4. The second summed term in Eq. ( 29) is based on the CWE representation of the field scattered by the equivalent source ∑ n∈Z q

(l,v) n ϕ H (1) n (k l f , r (l,v)
). The associated CWE coefficients of Eq. ( 29) have the similar relation,

a (l,v) m = a (l,v) m,inc + ∑ (l,v)∈D ∑ n∈Z q (l,v) n a (l,v) m,n . (30)
The solution to a

(l,v)
m,inc and a

(l,v)
m,n is obtained during the computation of V inc and G n (r, o (l,v) ). The combination of Eq. ( 22) and Eq. (30) yields the solution of q (l,v) m as

q D = I -Π D Λ D -1 Π D a D inc , (31)
where the vector

q D = [q (l,v) m ], a D inc = [a (l,v) m,inc ] and matrix Λ D = diag Λ (l,v) , Λ (l,v) = [a (l,v) m,n ], Π D = diag Π (l,v) , (l, v) ∈ D.
Finally, with the solution of q (l,v) m , the field solution is obtained by following Eq. ( 29).

Remark that, since only the coefficients of damaged fibers are concerned in Eq. ( 31) and the number of damaged fibers is usually small, the size of matrices in (31) is not large and the computational cost for the involved matrix inversion is low. Moreover, since V inc , G n , a D inc , Λ D are only related with the structure of the undamaged laminate and the incident wave, which are usually known, these quantities can be computed in advance and stored. 1. Initialize the set of layer indexes L = {1, 2, . . . , L}, M l = M max l , l ∈ L, P = 0, threshold τ, and position of the observation points (x l , z l ), x l = x 0 l , z l = (z 0 l + t l )/2. 2. Compute field values at (x l , z l ), l ∈ L, by CWE and PWE and denote the field solution by V l c and V l p . 3. If max{ V l p -V l c , l ∈ L} < τ, go to step 4; otherwise, update P l = P l + 1 and go to step 2. 4. Update M l = M l -1 for l ∈ L and compute corresponding V l p and

Truncation of expansions

V l c . 5. If V l p -V l c ≥ τ, l ∈ L, update M l = M l + 1 and L = L \ {l} . If L is empty, stop; otherwise, go to step 4.
PWEs and CWEs are infinite series and approximated by truncations ∑ P p=-P and ∑ M p=-M in the computation. There is a trade-off in the selection of the truncation number P and M. Larger P and M lead to a higher computational cost (especially considering hundreds even thousands of computation of the integrand for numerical integrations, e.g. in Eq. ( 17)), while the modeling accuracy cannot be guaranteed with small values of P and M. Moreover, too large M may lead to ill-conditioned matrices, which are inverted during the computation of scattering matrices.

In the authors' previous work [START_REF] Liu | Electromagnetic modeling of damaged single-layer fiber-reinforced laminates[END_REF], P and M for the single-layer laminate are estimated by a numerical strategy which benefits from the knowledge that field values computed by PWE and CWE in the vicinity of fibers are the same.

For multilayer structures, since the geometry and composing material might differ, optimal P and M might be changed accordingly from one layer to another. Thus, the numerical strategy is adjusted to fit the modeling of multilayer laminates. It is summarized in Table 1.

Step 2 and 3 are to select P, the value for different layers being chosen to be the same to facilitate mode matching on slab boundaries. Step 4 and 5 aim at reducing M l , which denotes the truncation number of CWE for the l layer, from its maximum value, which is set as

M max l = Int ℜ 4.05 × k l s c l 1/3 + k l s c l +7 [ 19
], 7 being a security number. Parameter τ, which determines the modeling accuracy, is a small positive constant.

Remark that τ cannot be too small in order to avoid too large P and M. Values from 10 -4 to 10 -3 seem to be robust and generate accurate enough field solution.

Numerical results

Truncation numbers P and M are first discussed with a 2-layer undamaged laminate. The material of the slab is with relative permittivity ε 1 s , ε 2 s = 3.6 + i0.072 and the fiber materials are with

ε 1 f = 6, ℜ(ε 2 f ) = 12, conductivity σ 2 f = 330 S/m.
The fiber period is d = 0.1mm and the thickness of each slab is d. Fibers are placed in the middle of each layer (i.e., z (l,0) = (t l + t l+1 )/2) and no shift appears between the fiber arrays (i.e., x (1,0) = x (2,0) ). The laminate is illuminated by a plane wave E inc = e -iα 0 (z-d/2) , α 0 = k 0 s sin θ inc , θ inc = π/6, with the TE polarization. Setting Validation of the proposed computational approach is performed with a 3-layer laminate sandwiched by air. Let the fiber period d = 0.1mm, radius c l = (l + 1)d/10, l = 1, 2, 3, and the origin of the local coordinates systems (x (l,0) , z (l,0) ) as z (l,0) = (t l + t l+1 )/2, t l = 3d/2ld. Fiber arrays in 2nd and 3rd layers are shifted with x (2,0) = 0.1d,

x (3,0) = -0.1d. Fiber materials are with ε 1 f , ε 3 f = 6, ℜ(ε 2 f ) = 12, conductivity σ 2 f = 330 S/m, while the slabs are slightly lossy (loss tangents of 0.02) with ε 1 s , ε 3 s = 2.8 + i0.056, ε 2 s = 3.6 + i0.072. Permeabilities µ l j = 1, j = s, f . The parameter τ in Table 1 is set as 10 -4 to ensure the validation accuracy and the truncation number selected by the truncation scheme is P = 5, M 1 = 10, M 2 = 11, M 3 = 13, M l the optimal M for the l-th layer. Set the line source at o (1,0) and run COMSOL with default mesh settings [START_REF]COMSOL Multiphysics R v5.3[END_REF] (physics-controlled mesh with "Extremely fine" element size). Since scattered fields are often collected above the laminate, fields are sampled along the line -5d ≤ x ≤ 5d, z = d, which is close to the laminate boundary so that evanescent components are counted in the validation. The number of fibers, 43, is decided by observing the convergence of the field solution. As shown in Fig. 9, with 41 fibers, adding one more fiber at each end of the laminate, the total-field solution generated by COMSOL has almost no differences. In contrast, with 15 fibers, since the the fibers at the end of the laminate locate close to the emitting source and the sampling region, their scattering effects cannot be neglected and, as a consequence, significant differences are observed at tails when putting more fibers. The field solution with TE polarization has a similar convergence. The solution to the Green's function computed by the presented approach is validated. As observed from Fig. 10, Challenges in imaging are numerous. Generally, sources and receivers can only be put above the laminate in the testing system to collect the information. Yet, for the damages inside deep layers, their impacts on the scattered field might be reduced significantly due to losses and, therefore, these damages should not be easily identified. Investigations here are led on a 4-layer laminate by moving the missing fibers from layer to layer. Each layer of the laminate enjoys the same structure and composition. Maps of anomalous field are found in Fig. 14, where only amplitudes smaller than 0.02 are displayed for better appraisal. With graphite fibers and epoxy slabs, due to lossy fibers, the amplitude of the anomalous field in the upper half space significantly reduces as the missing fibers move deeper and deeper. When the missing fibers show up at the bottom layer, one has almost no hint of damages when above the laminate. However, changing the fiber material to lossless glass, as shown in the left figures of Fig. 14, the field decays much slower in the propagation. As a result, the impacts of damages can be observed above the laminate even if the missing fibers are located in the bottom layer.

•10 -2 1 •10 -2 1 •10 -2 1 •10 -2 1 •10 -2 1 •10 -2 1 •10 -2 1 •10 -2

Glass Graphite

Conclusion

The modeling of damaged periodic fibered laminates and the solution to the associated Green's function are the topics of concern herein. Modeling the undamaged laminates by computing the scattering matrices with a recursive scheme, the solution to the Green's function is obtained with the array scanning method, where the hypothesized phased array is first analyzed with the quasi-periodic theory and the final solution follows from integration over the first Brillouin zone. The field disturbance by missing, displaced, shrunk, and expanded fibers and circular inclusions is modeled as equivalent sources. According to the scattering linearity, the field solution to the damaged laminate is a summation of responses due to the exterior illuminating source and equivalent ones. The proposed modeling approach is validated by comparing results from COMSOL or the Born approximation.

Circular damages associated with fibers are those concerned in this contribution. For damages within a slab but outside fibers, the equivalence theory is to applied and solving the corresponding Green's function is the first step. For defects with arbitrary shapes, a method of moments within a proper field integral formulation can apply. The direct modeling of twisted (instead of straight) fibers, which are randomly distributed in the background slab, is another strong challenge ahead of the authors.

Behaviors of photonic crystals may be analyzed with the equivalence theory due to similar structures. For instance, the so-called defect mode may be observed by selecting frequencies for a singular matrix in Eq. (31) [START_REF] Faggiani | Lower bound for the spatial extent of localized modes in photonic-crystal waveguides with small random imperfections[END_REF].

Appendix A. Some variable definitions

In ( 2), E u l , u = 1, 2, 3, 4, are defined as the diagonal matrix of E u,p l and E 1 l,p = 2w l-1 p e iα p (x l -x l-1 ) e i w l-1 p + w l p E 4 l,p = 2w l p e -iα p (x l -x l-1 ) e i[β l p (t l-1 -z l )-β l-1 p (t l-1 -z l-1 )] w l-1 p + w l p where w l p = β l p /(µ 0 µ l s ) for TM and w l p = β l p /(ε 0 ε l s ) for TE, µ 0 , ε 0 the permeability and permittivity of air.

Fig. 1 .

 1 Fig. 1. Sketch of the damaged laminate with defects of concern, L = 3.

Fig. 2 .

 2 Fig. 2. Coefficients of plane-wave expansions in the l-th layer for upgoing and down-going waves.

Fig. 3 .

 3 Fig. 3. Field representation by cylindrical-wave expansions in the (l, v)-th fiber.

m

  stands for the wave propagating towards the fiber center, b (l,v) m and q (l,v) m denote the field scattered by the fiber and the inner source (e.g., the analytical line source inside the fiber for the scalar Green's function), respectively. The total field outside fibers of each layer is a summation of fields scattered by the corresponding upper and lower slab boundaries with PWE coefficients f l+ p and g l- p and fields scattered by all fibers with CWE coefficients b (l,v) m . The CWE representation of the latter is with the identity ∑

Fig. 4 .

 4 Fig. 4. Equivalence of fiber arrays as planar boundaries with scattering matrices E u l and F u l , u = 1, 2, 3, 4, defined by Eq. (2) and (11).

Fig. 5 .

 5 Fig. 5. Phased line-source array, solid dot standing for the true source and hollow ones for fictitious sources.

Fig. 7 .

 7 Fig. 7. Line graphs of CWE coefficients b m , M = 12, two-layer laminate with the plane-wave illumination, θ inc = π/6, λ inc = d.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Geometry of COMSOL module, 43 fibers in each layer of the laminate, fiber radius c 1 = 0.2d, c 1 = 0.3d, c 1 = 0.4d.

Fig. 10 .

 10 Fig.10. Validation of the computation of the Green's function by comparisons with COMSOL, r s = (0, 0), λ inc = d, z = d, "COMSOL D" denoting COMSOL solution with the default mesh settings and "COMSOL M" with manual settings.

2 )Fig. 12 .

 212 Fig.[START_REF] Watanabe | Accurate analysis of electromagnetic scattering from periodic circular cylinder array with defects[END_REF]. Validation of modeling the damaged laminate through comparisons with COMSOL, r s = (0, d), λ inc = d, "COMSOL D" denoting COMSOL solution with the default mesh settings and "COMSOL M" with manual settings.

Fig. 13 .

 13 Fig. 13. Anomalous field computed by the proposed approach and the Born approximation, the (1, -3)-th and (3, 3)-th fibers are missing, r s = (0, d), λ inc = 100d, TM, z = 1.5d.

Fig. 10 (

 10 Fig. 10(a) Fig. 10(b) Fig. 12(a) Fig. 12(b) COMSOL D 20 18 16 18 COMSOL M 53 130 86 96 Proposed 11 15 0.26 0.23

Fig. 14 .

 14 Fig. 14. Maps of anomalous field with the (l, 2)-th and (l, 2)-th fibers missing, l = 1, 2, 3, 4, line source at (0, d), λ inc = d, TM.

-i2β l- 1 p

 1 [β l p (t l-1 -z l )-β l-1 p (t l-1 -z l-1 )] (t l-1 -z l-1 )

Table 1 .

 1 Numerical strategy to truncate PWE and CWE.

Table 2 .

 2 Comparison of computation time (in seconds)

with default mesh settings, although the field solution of the proposed approach fits COMSOL well with TM polarization, differences appear with TE polarization. The "maximum element growth rate", the parameter in COMSOL determining the maximum rate at which the element size can grow [START_REF]COMSOL Multiphysics R v5.3[END_REF], is then reduced from the defaulted 1.1 to 1.01 to have a finer mesh to more accurately represent the geometry. The field solution associated with this manual setting is denoted by "COMSOL M" in Fig. 10 11, where the amplitude is limited at 0.3 for easier appraisal.

In the TM case, most of the energy is reflected into the upper half space and little penetrates into the laminate, while more energy propagates inside in the TE case.

To validate the computational approach for modeling defects, 500 values are respectively sampled along the line z = 0, -5d ≤ x ≤ 5d and the line z = -d, -5d ≤ x ≤ 5d. The comparisons with results from COMSOL are shown in Fig. 12. Two solutions again are generated by COMSOL with default mesh settings and manual settings (maximum element growth rate changed as 1.01), respectively. Slight differences appear among the three solutions. However, the fields computed by the proposed approach is between the two COMSOL solutions for most points, despite strong fluctuations due to the size of component, i.e., slab thickness and fiber radius, close to the wavelength.

When the wavelength is much larger than the component size, the field inside the fiber should only be slightly impacted by the damages and the Born approximation can be used as described in Appendix B. With the same composite but setting λ inc = 100d and assuming only missing fibers exist, the comparison of computed anomalous field is shown in Fig. 13. The peak value appears at the position corresponding to the (3, 3)-th missing fiber due to its larger radius (c 3 = 0.4d while c 1 = 0.2d).

The modeling efficiency is indicated by Table 2, which presents the computational time of reaching line graphs

Appendix B. Application of the Born approximation

With missing fibers, the field disturbance can be valued by the Fredholm integral equation [START_REF] Groby | Localization and characterization of simple defects in finite-sized photonic crystals[END_REF]

where D (l,v) is the cross-section of the (l, v)-th fiber, V is the field solution for the sound laminate and Ṽ for the damaged one. With large wavelengths, Ṽ (r s ) ≈ V (r s ), thus

Denote g 0 p as the PWE coefficients for G when the observation point is above the laminate and C

(l,v) m as the CWE coefficients for V (r s ). Substituting the expression of G and V into the integrand, Eq. (B.2) is rewritten as

is the m-th column of Φ (l,v) ,

where Γ (l,v) = e i(α p (x-vd)+β 0 p (z-t 1 )) and J = diag k l f c l J m+1 (k l f c l )J m (k l s c l )k l s c l J m (k l f c l )J m+1 (k l s c l ) . The integration in (B.3) is computed by numerical methods.