domain: spi math phys

text:
Luiz F T Fernandez

Murat Bronz
email: murat.bronz@enac.fr

Thierry Lefebvre
email: thierry.lefebvre@onera.fr

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

FROM MDO TO MANUFACTURING: APPLICATION CASE FOR UNMANNED AERIAL VEHICLES

Keywords: Additive Manufacturing, MDO of UAVs, UAV design and manufacturing

This paper explores the concept of coupling an MDO problem directly to a manufacturing process. The design parameters required for various disciplinary analyses are used to automatically generate geometries that comply with the selected manufacturing processes. Such geometries are then manufactured to validate our proposal. This strategy reduces the manufacturing time of unmanned aerial vehicles (UAVs) as the CAD modeling becomes automatic and parametric, also allowing for easier comparison of different geometries. It can also potentially improve the MDO process by closing the information loop with manufacturing constraints. We employ the Engineering Sketch Pad (ESP) and leverage from additive manufacturing techniques usually applied in UAVs. The choice of ESP facilitates reproducibility, as it is an operating system agnostic open source CAD, and gradient based optimization, due to its capability of computing gradients of the geometric outputs with respect to the design parameters. The manufacturing processes of wings and propellers are addressed. For the same representative wing geometry, we present different modeling strategies for mass and inertia prediction and for manufacturing using 3D printing. ESP is also employed to predict wing mass and inertia. The maximum difference in weight between the manufactured wing and its predicted value using ESP is roughly 7%. The propeller is 3D printed in resin with stereolithography (SLA) technique, and is defined by means of its airfoils, radius, and chord and twist angle distributions. We conclude about the applicability of the presented strategy, as well as its potential and limitations. All the scripts are shared with the community so researchers can apply it to validate their own MDO problem.

INTRODUCTION

With the establishment of unmanned aerial vehicles (UAVs), a significant increase in the interest in such kind of vehicle has been observed in recent years for several applications. Ranging from widespread ideas such as delivery of goods or search and rescue to more specialized ones, such as fire fighting [START_REF] Ausonio | Drone swarms in fire suppression activities: A conceptual framework[END_REF][START_REF] Aydin | Use of fire-extinguishing balls for a conceptual system of drone-assisted wildfire fighting[END_REF][START_REF] Akhloufi | Unmanned aerial vehicles for wildland fires: Sensing, perception, cooperation and assistance[END_REF], debris mapping [START_REF] Taddia | Uav approach for detecting plastic marine debris on the beach: A case study in the po river delta (italy)[END_REF], covid mitigation [START_REF] Restás | Drone applications fighting covid-19 pandemic;towards good practices[END_REF], and even nutrition delivery in the form of edible vehicles [START_REF] Kwak | Towards edible drones for rescue missions: design and flight of nutritional wings[END_REF]. Within different configurations of UAVs, winged designs stand out when higher flight range and endurance are needed. Several manufacturing techniques have been applied for winged UAVs fabrication. In terms of wing manufacturing, hobbyists predominantly utilize balsa wood, which was primarily employed during the early stages of the UAV era. This conventional approach has subsequently been enhanced by combining balsa with composite materials, as in Chung et al. [START_REF] Chung | Design, manufacturing, and flight testing of an experimental flying wing uav[END_REF], and even full composite wings as in Nugroho et al. [START_REF] Nugroho | Manufacturing process and flight testing of an unmanned aerial vehicle (uav) with composite material[END_REF]. Full composite wings tend to be light and of high quality, but the manufacturing process is complicated and time consuming. Recently, there has been attention towards additive manufacturing techniques due to their ability to save time and material, fabricate complex geometries, and achieve higher levels of automation. Such techniques have been successfully applied in several cases. Bronz et al. [START_REF]Mission-Oriented Additive Manufacturing of Modular Mini-UAVs[END_REF] presented a 3D printing strategy for UAVs using Onyx material, a blend of carbon fiber and nylon, which provided the desired structural properties and surface quality while keeping a reasonable weight. Taylor et al. [START_REF]Design Optimization, Fabrication, and Testing of a 3D Printed Aircraft Structure Using Fused Deposition Modeling[END_REF] and Muir et al. [START_REF]The Use of MDO and Advanced Manufacturing to Demonstrate Rapid, Agile Construction of a Mission Optimized UAV[END_REF] coupled optimization processes to additive manufacturing techniques in order to obtain an aircraft structure and a mission optimized UAV, respectively. Laliberté et al. [START_REF] Laliberté | Design and manufacturing of biologically inspired micro aerial vehicle wings using rapid prototyping[END_REF] used 3D printing to build a biologically inspired micro aerial vehicle. In 2015, Aurora Flight Sciences also presented the first jet-powered 3D-printed UAV [START_REF] Atlas | World's first jet-powered, 3d-printed uav debuts at dubai airshow[END_REF]. While several other applications of additive manufacturing for UAVs can be found in the survey by Goh et al. [START_REF] Goh | Additive manufacturing in unmanned aerial vehicles (uavs): Challenges and potential[END_REF], this technique has also been successfully employed for teaching activities [START_REF]Aero-engine model design and assembly practice teaching based on 3D printing[END_REF] and even in experiments in hypersonic conditions [START_REF] Rêgo | Ground experimentation with 3d printed scramjet inlet models at hypervelocities[END_REF]. For propellers, there are several off-the-shelf options commercially available. Their good quality combined with small prices is usually enough for the UAV community, which has not yet explored propeller design and manufacturing with the exception of a few studies, as in Rutkay and Laliberté [START_REF] Rutkay | Design and manufacture of propellers for small unmanned aerial vehicles[END_REF]. We already employed additive manufacturing for propellers in [START_REF] Fernandez | Assessment of methods for propeller performance calculation at high incidence angles[END_REF], where a specific geometry was needed in order to compare different aerodynamic methodologies. The higher level of automation from 3D printing also facilitates the integration between design and fabrication processes. Such integration might be beneficial in two ways: enriching design and optimization process with manufacturing constraints, anticipating possible problems or even unfeasible designs; and reducing the time needed to complete the full cycle between design and flight. Regarding the insertion of manufacturing information within a multidisciplinary design and optimization (MDO) process, Laan et al. [START_REF]Bringing Manufacturing into the MDO domain using MBSE[END_REF] used model-based system engineering techniques to account for the manufacturing cost of high lift devices in the design of experiments performed for the initial sizing of an aircraft. Doneli et al.. [START_REF] Donelli | A model-based approach to trade-space evaluation coupling design-manufacturing-supply chain in the early stages of aircraft development[END_REF][START_REF] Donelli | Value-driven model-based optimization coupling designmanufacturing-supply chain in the early stages of aircraft development: Strategy and preliminary results[END_REF] also addressed, within the MDO, the choice of raw materials, manufacturing and assembly processes as well as the supply chain, ensuring a more realistic and comprehensive optimization process. In this paper, our focus is to present a viable way to automatically connect the output of a given MDO process to the fabrication of a UAV. The optimization loop is represented by means of its output design parameters or geometry, and the fabrication is performed with 3D printing. We leverage from the size and affordable characteristics of UAVs to fully test our proposal, presenting the entire process, from modeling to real manufacturing, for wings and propellers. This paper is organized as follows: in Section 2 we present the context in which this study takes place. Section 3 explains how we share the data. The methodology and results are presented in Sections 4 and 5, respectively. Our conclusions and perspectives are then presented in Section 6.

CONTEXT

This paper is part of a research line dedicated to the application of MDO strategies to small vertical takeoff and landing (VTOL) vehicles. Apart from being useful for different applications, such vehicles also offer the possibility of experimentally validating design strategies and decisions due to their low cost and easier manufacturing and flight testing (compared to bigger vehicles) procedures. On the other hand, manufacturing the vehicle is usually a time consuming task that includes not only the fabrication itself but also the generation of 3D models, needed to guide the integration of all the systems. Building a good model often requires experience to avoid several iterations and minimize cost. In [START_REF] Fernandez | Development of a missiontailored tail-sitter mav[END_REF], we developed a simplified MDO problem to design a vehicle for the 2022 International Micro Aerial Vehicle Competition. The outcome of this problem was a set of parameters, shown in Tab.1. In order to "transform" such numbers into a real vehicle, we manually built the 3D model shown in Fig. 1. Such process consumed significantly more time than running the MDO. Similarly, in [START_REF] Fernandez | Assessment of methods for propeller performance calculation at high incidence angles[END_REF] we fabricated propellers to compare different aerodynamic techniques for evaluating their performance. After designing one by defining its airfoils, twist and chord distributions, a manual 3D model was also needed to design the hub and integrate it with the blades in order to have a usable propeller. This papers addresses this gap between the output of a given MDO problem and the generation of a manufacturing compliant 3D model, allowing for time and cost reduction and permitting less experienced designers to build models more easily.

DATA SHARING

In the subsequent sections, we outline the methodology utilized and present the results obtained. All files, ESP geometries and Python scripts, are shared so the community can use and improve the methods. They are available in ENAC drones group GitHub at https://github.com/enac-drones/AeroBest23_MDO_UAV/tree/paper_codes/. We provide two examples for wings manufacturing, one simpler without twist and using a NACA airfoil (under Wo twist folder) and other with twist and Kulfan airfoils (under With twist folder). The scripts using ESP's Python interface for wing weight estimation are also provided. The Python scripts for blade generation and ESP files for propellers are available under the folder propellers.

METHODOLOGY

Instead of defining an MDO process, we only represent the output of such kind of problem, ensuring that the same strategy can be used for any process with a similar parametrization. As explained by Hajdik et al. [START_REF]Aerodynamic Shape Optimization with CAD-Based Geometric Parameterization[END_REF], despite the wide application of free form deformation (FFD) in MDO, as in [START_REF] Kenway | A CAD-free approach to high-fidelity aerostructural optimization[END_REF][START_REF] Shi | Natural laminar flow wing optimization using a discrete adjoint approach[END_REF][START_REF] Bons | Aerostructural Wing Optimization of a Regional Jet Considering Mission Fuel Burn[END_REF], coupling FFD geometries in general design and manufacturing workflows is not straightforward, as information might be lost throughout this process. So, a CAD tool is employed, specifically the Engineering Sketch Pad (ESP) [START_REF] Dannenhoffer | The engineering sketch pad[END_REF], an open source and operating system agnostic software. ESP incorporates drivers for aerodynamic and structural tools. Additionally, it is compatible with OpenM-DAO [START_REF] Gray | OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization[END_REF] and natively allows for gradient calculation. In order to test our proposal for wings, we select a representative set of design parameters and present the different modeling strategies needed for manufacturing and for mass/inertia estimation. For propellers, the manufacturing process is simpler and relies on BladeX by Gadalla et al. [START_REF] Gadalla | BladeX: Python Blade Morphing[END_REF].

Wing and propeller parametrization

Each wing section is defined by the parameters shown in Tab.2, where the number of sections can also change. ESP has primitives for different airfoil families, but we employ NACA family for simplicity. The structural parameters, shown in Tab.3, are defined for the whole wing. Spars are assumed to be cylindrical as off-the-shelf carbon tubes are usually employed in this form for small UAVs. The rib offset is the angle between the ribs and the chord line. For propeller parametrization we employ three parameters: airfoil(s), chords and twist distributions along the radius.

Manufacturing setup

The wings are manufactured with a Raise3D N2 Pro by Raise 3D 1 . We use the lightweight PLA from colorFabb 2 for all tests. Once the geometry has been defined, as an STL file for our case, it needs to be "sliced" in order to be printed. There are several open source slicers. We use the IdeaMaker 3 due to its compatibility with our printers. For the propeller, the 3D printer Form 3+, from FORMLABS 4 is employed. The material used is the "Tough" resin, provided by the same company. To prepare the printing process we use the PreForm 5 , also open source.

Modeling wings for manufacturing

The 3D printing process is based on the one presented in [START_REF] Station | How i designed a 3d printed wing[END_REF] and has already been used in [START_REF] Fernandez | Development of a missiontailored tail-sitter mav[END_REF]. This fabrication strategy ensures a more continuous printing, minimizing the need to remove the nozzle throughout the process, which tends to increase the weight of the piece and reduce its quality. It roughly consists in modeling the wing as a solid body while keeping the ribs as empty cavities. Figure 2 shows the obtained model. ESP is capable of directly outputting the resulting geometry as an STL file, which is then used to obtain the G code needed to print the part. Manufacturing wings with high twist angles and cambered airfoils might require some adaptations to the original method in order to correctly generate the ribs. An example for such case is also available for comparison.

Modeling wings for mass and inertia estimation

Understanding the modelization strategy for mass and inertia estimation is straightforward. The CAD model needs to represent the final manufactured body as accurately as possible, as seen in Fig. 3. With such model and the material definition, shown in Tab.4, we calculate mass and inertia with the Masstran feature present in PyCaps [START_REF] Durscher | pycaps: A python interface to the computational aircraft prototype syntheses[END_REF]. The membrane thickness was obtained after averaging different printed pieces.

Modeling propellers

The propeller process starts with the definition of one blade using BladeX [START_REF] Gadalla | BladeX: Python Blade Morphing[END_REF]. The blade is defined by means of its airfoils and chord and twist distributions. Once the blade is generated, it is exported as an IGES file. The ESP script then imports this file to generate both blades. It also models the hub and exports the propeller as an STL file to be used for printing.

RESULTS

In the following sections we present the results obtained with the manufacturing processes as well as some aspects regarding their validation. We present the comparison between predicted and measured weight for wings and discuss the time needed to obtain the printed pieces.

Wings

Figure 4 shows two printed wings, different from each other because of the airfoils, geometric twist and number of ribs, and also a zoom in the airfoil and wing internal structure. The geometry of the wings can for now be validated by means of measuring span and chord. A further analysis with 3D scanners is left for future work. However, any mismatch between the printed wing and the one defined during the optimization is most likely caused by the printing itself, as the geometry generated with ESP matches the one intended. We also compare the expected and observed mass using Masstran from PyCaps [START_REF] Durscher | pycaps: A python interface to the computational aircraft prototype syntheses[END_REF]. Each wing mass was measured and the results are shown in Fig. 5. Table 5 shows the difference between predicted and measured mass. The accuracy of the outcome is highly dependent on the quality of the mass models. For more complex wings with, for example, command surfaces, servos, and hinges, the prediction might loose precision, depending on how such geometry is modeled. For our case, the simple wing geometry is based on our vehicle Falcon [START_REF] Fernandez | Development of a missiontailored tail-sitter mav[END_REF] shown in Fig. 1, a UAV with four motors and no command surfaces. As Masstran also is capable of calculating inertia, this tool can be very useful for flight dynamics analysis. However, any validation with this regard is left for future work. The CAD model and "manufacturing view" are automatically generated once the geometry is known. So the total time to obtain the piece is driven by the manual slicing process, which is straightforward, and the printing time. The printing time depends highly on the geometry, printer, and material. For both wings manufactured for this paper, the process took roughly 4 hours each considering the setup presented in Section 4.2. This result shows a potential of reducing the manufacturing time in days where the final goal is to have only the printing time as constraint. However, the real reduction is hard to estimate because it highly depends on the designer's ability to generate the model and on the complexity of the geometry.

Propellers

Figure 6 shows the obtained propellers. The printer generates a supporting structure to ensure quality, as shown in Fig. 6a. This structure needs to be removed and the propeller cleaned, as shown in Fig. 6b. We can also validate propellers by measuring the chords and twist angles. The 3D scan technology is even more important for them as it is difficult to measure the angles properly. However, the similar performance in experiments and analytical methods observed in [START_REF] Fernandez | Assessment of methods for propeller performance calculation at high incidence angles[END_REF] indicates that the propellers are close to the original geometry definition.

CONCLUSION

In this paper, we present a preliminary way of integrating the manufacturing process of wings and propellers to an MDO framework for UAVs. We employ ESP to automatically generate the fabrication and mass prediction geometric "views". Such geometries are then used to manufacture two sample wings and assess their predicted weight. The weight prediction shows good accuracy for our case, mainly because of the simplicity of the wings analyzed. We also present a similar strategy for propeller manufacturing. Such strategy has also been used in [START_REF] Fernandez | Assessment of methods for propeller performance calculation at high incidence angles[END_REF] and the performance obtained in wind tunnel tests was close to the one expected. It is expected that this solution can improve MDO processes by providing new manufacturing constraints and anticipating potential issues. It can also reduce the time needed to build test vehicles to perform flight and/or wind tunnel test campaigns, thus accelerating the design process. Every script used for this paper is open source so the community can use the methods to build UAVs and subscale vehicles to improve MDO processes.

Figure 1 :

 1 Figure1: 3D model and real vehicle from[START_REF] Fernandez | Development of a missiontailored tail-sitter mav[END_REF]

Figure 2 :

 2 Figure 2: Modeling the wing to comply with the manufacturing technique: solid body with empty internal ribs (in detail on the right) and spars.

Figure 3 :

 3 Figure 3: Modeling the real printed wing for mass and inertia estimation, with internal ribs being detailed on the right.

  (a) Printed wings with and without twist.(b) Airfoil and internal structure.

Figure 4 :

 4 Figure 4: Printed wings.

( a )

 a Twisted wing with 3 ribs per side. (b) Straight wing with 10 ribs per side.

Figure 5 :

 5 Figure 5: Wing mass measurements.

( a )

 a Printed propellers straight out of the printer. (b) After the cleaning and "sanding" process.

Figure 6 :

 6 Figure 6: Printed propellers.

Table 1 :

 1 Optimal design from[START_REF] Fernandez | Development of a missiontailored tail-sitter mav[END_REF]

	Design variables	Value
	Root chord	0.12 m
	Tip chord	0.095 m
	Wing span	0.7 m
	Cruise α	10 •
	Turn α	8 •
	Battery mass	0.2 kg

Table 2 :

 2 Wing section parameters.

	Parameter	Unit
	Chord	[m]
	Position in X, Y, and Z	[m]
	Twist angle	[deg]
	Airfoil location of max camber [percent chord]
	Airfoil maximum thickness	[percent chord]
	Airfoil maximum camber	[percent chord]

Table 3 :

 3 Wing structural parameters.

	Parameter	Unit
	Number of ribs	-
	Thickness of the ribs	[m]
	Ribs offset	[deg]
	Number of spars	-
	Diameter of each spar	[m]
	Location of each spar [percent chord]

Table 4 :

 4 Masstran material definition.

	Parameter	value
	Material	Lighweight PLA
	Material density Material type	390 [ kg m 3 ] Isotropic shell
	Membrane thickness	0.74[mm]

Table 5 :

 5 Mass prediction results.

	Number of ribs Masstran result [g] Measured mass [g] Error [%]
	3	25.30	26.1±0.1	3.16
	10	31.53	29.5±0.1	6.88

ACKNOWLEDGEMENTS

This work is part of the activities of ONERA-ISAE-ENAC joint research group.


Étude de la Projection des Forces Magnétiques pour le Calcul des Vibrations dans les Machines Électriques

Raphaël PILE 1,2,3 , Jean LE BESNERAIS 1 , Yvonnick LE MENACH 2 , Guillaume PARENT 3 , Thomas HENNERON 2 , and Jean-Philippe LECOINTE 3

1 EOMYS ENGINEERING, Lille-Hellemmes 59260, France 2 Univ. Lille, Arts et Metiers ParisTech, Centrale Lille, HEI, EA 2697 -L2EP -Laboratoire d'Électrotechnique et d'Électronique de Puissance, F-59000 Lille, France 3 Univ. Artois, UR 4025, Laboratoire Systèmes Électrotechniques et Environnement (LSEE), F-62400 Béthune, France De nombreux travaux récents traitent de calcul des forces magnétiques pour le calcul du bruit et des vibrations dans les machines électriques. La principale difficulté réside dans la méthodologie d'identification des ondes de forces rotatives les plus bruyantes. Les deux méthodes utilisées pour calculer les forces magnétiques sont le Principe des Travaux Virtuels (PTV) et le Tenseur de Maxwell (TM). Cette communication traite des vibrations obtenues sur une machine à aimants en surface grâce à une résolution électromagnétique Eléments Finis (EF) 2-D et une simulation mécanique harmonique EF 3-D. Les résultats de deux calculs distincts sont comparés : un cas avec le PTV, et l'autre cas basé sur le TM. En particulier, le but de ce travail est d'étudier l'effet de modulation mécanique. Cet effet a pour conséquence de filtrer les ondes rotatives de forces à cause de la nature discrète des dents du stator. Deux méthodes différentes pour tenir compte de l'effet de modulation sont comparées par rapport au PTV pour la machine étudiée. Ce travail montre qu'une des deux méthodes est plus pertinente que l'autre.

Vibration, Machines électriques, Force magnétique, Magnétomécanique, Tenseur de Maxwell, Principe des travaux virtuels.

INTRODUCTION

L'analyse du comportement vibro-acoustique des machines électriques est le sujet de nombreux travaux récents. Le problème, toujours pas clairement résolu, est même accru dans la mesure où l'accroissement des puissances massiques et volumiques impacte les niveaux vibratoires et sonores. Ces phénomènes doivent être pris en compte dès la conception, et l'un des moyens pour y parvenir est d'avoir recours à la simulation. Un point crucial de ces simulations est la modélisation des forces Electro-Magnétiques (EM).

Dans ce contexte, la méthode la plus couramment utilisée est le Tenseur de Maxwell (TM) : une force surfacique équivalente est estimée dans l'entrefer à partir de la connaissance du champ magnétique à cet endroit. Le TM a l'avantage d'être compatible avec les méthodes analytiques, semi-analytiques et Éléments Finis (EF). Grâce à l'analyse harmonique, le TM permet d'identifier les ondes de force rotatives les plus gênantes, et ainsi de les résoudre avec un minimum d'effet sur les performances. Pour réaliser cette analyse, le TM est décomposé dans la base de Fourier selon ses ordres angulaires (ou nombre d'onde angulaire, ou encore pulsation spatiale). Dans le cas d'un couplage avec une simulation mécanique EF, le TM est généralement intégré pour appliquer une force résultante par dent [START_REF] Saito | Vibration Prediction Method of Electric Machines by using Experimental Transfer Function and Magnetostatic Finite Element Analysis[END_REF].

Une autre méthode utilisée est le Principe des Travaux Virtuels (PTV) : appliqué à la méthode des EF, le PTV donne une résultante d'effort magnétique nodale [START_REF] Coulomb | Finite Element Implementation of Virtual Work Principle for Magnetic or Electric Force and Torque Computation[END_REF]. Le PTV est la méthode de référence pour le calcul des forces locales (notamment sur le bout des dents), mais il ne permet pas d'identifier simplement les sources de bruit EM comme dans le cas du TM. Pour ces raisons, cette étude propose de constuire un cas qui permet de simuler directement le déplacement issu d'un chargement PTV (non-modifié). Pour réaliser cela, on conserve le maillage magnétique pour la simulation EF mécanique. La simulation magnétique a été construite afin de limiter l'impact des irrégularités et des asymétries du maillage magnétique sur le résultat mécanique. Le résultat vibratoire est comparé à différents chargements TM avec un maillage mécanique.

Cette simulation PTV de référence permet d'étudier de manière purement numérique l'effet de modulation mécanique [START_REF] Fang | Modulation Effect of Slotted Structure on Vibration Response in Electrical Machines[END_REF][START_REF] Wang | Analysis of Zeroth-Mode Slot Frequency Vibration of Integer Slot Permanent-Magnet Synchronous Motors[END_REF] : des relevés expérimentaux ont mis en évidence un repliement des ordres élevés de forces sur les ordres plus faibles. Cette étude propose de vérifier numériquement l'ampleur de cette effet. Pour cela les modèles de modulation proposés sont comparés par rapport au calcul direct basé sur le PTV, en particulier pour un nombre d'onde de force égal au nombre d'encoches au stator Z s = 36 comme dans l'étude [START_REF] Wang | Analysis of Zeroth-Mode Slot Frequency Vibration of Integer Slot Permanent-Magnet Synchronous Motors[END_REF].

CALCULS DES FORCES D'ORIGINE MAGNÉTIQUE 2.1. Tenseur de Maxwell

Cette section présente le principe du calcul des forces magnétiques avec le tenseur de Maxwell, ainsi que la méthode d'application des efforts sur le modèle mécanique. La méthode classique consiste à prendre un contour circulaire dans l'entrefer sur lequel est défini l'effort surfacique à la position angulaire α ∈ [0, 2π] et l'instant t à partir de la perméabilité du vide µ 0 et de la densité de flux magnétique radial B r et tangentiel B t [START_REF] Pile | Comparison of Main Magnetic Force Computation Methods for Noise and Vibration Assessment in Electrical Machines[END_REF] :

σr(α, t) = Br(α, t) 2 -B θ (α, t) 2µ0 2 (1)
En appliquant une décomposition de Fourier selon le temps et la position angulaire, on obtient les harmoniques selon :

σr(α, t) = n ω σr(n, ω)e jnα (2)
avec j le nombre imaginaire, ω la pulsation, n le nombre d'onde. Ainsi, cette force surfacique est définie dans l'entrefer sans être rattaché à une surface. Une méthode consiste alors à intégrer le morceau de l'onde qui correspond à chaque dent de position angulaire

α i , i ∈ [1, Z s ] [1]
en ne considérant qu'un seul ordre spatial σ(n, ω). On obtient alors une onde de force radiale tournante intégrée par dent.

Effet de modulation

De récentes études ont mis en évidence la présence d'un effet dit de modulation : un nombre d'onde n élevé va être aussi perçu par la culasse comme un nombre d'onde m = n ± Zs plus petit. Ainsi le déplacement Y r,i (n, ω) peut être calculé grâce à Y r,i (m, ω) selon :

Yr,i(n, ω) = K mod (n, m)Yr,i(m, ω) (3)
L'amplitude K mod du repliement est différente selon la source. Cela dépend principalement de la méthode d'échantillonage de l'onde de force. Ainsi, dans la méthode de [START_REF] Fang | Modulation Effect of Slotted Structure on Vibration Response in Electrical Machines[END_REF] appelée TM-FANG, la force radiale tournante est intégrée sur le pas dentaire, de sorte que :

K mod (m, n) = m n sin(nπ/Zs) sin(mπ/Zs) (4)
Tandis que dans la méthode appelée TM-WANG [START_REF] Wang | Analysis of Zeroth-Mode Slot Frequency Vibration of Integer Slot Permanent-Magnet Synchronous Motors[END_REF], l'intégration se fait par rapport à la l'ouverture angulaire β du bout de la dent :

K mod (m, n) = m n sin(nβ/2) sin(mβ/2) (5)
Néanmoins aucune de ces deux études n'effectue de comparaison avec le PTV : c'est ce que la Section 3 propose de faire.

Principe des Travaux Virtuels

Le PTV utilisé est le même que dans [START_REF] Pile | Comparison of Main Magnetic Force Computation Methods for Noise and Vibration Assessment in Electrical Machines[END_REF]. Un exemple de résultat des forces radiales à un instant t obtenu pour cette étude est présenté dans la Fig. 1. Afin de réaliser une simulation mécanique harmonique avec le même chargement, le maillage de triangles 2-D est extrudé en 3-D sur 8 couches de même longueur, afin d'obtenir un maillage prismatique. Le chargement PTV est ensuite dupliqué sur les noeuds, excepté aux extrémités pour éviter les effets de bords.

APPLICATION À UNE MACHINE SYNCHRONE

Une comparaison numérique des méthodologies est réalisée sur une machine synchrone à aimants en surface. Avec le PTV, on conserve le maillage magnétique pour la simulation EF mécanique grâce à une extrusion, alors qu'avec le TM on utilise un maillage mécanique composé d'hexaèdres. La topologie choisie est constituée de 36 encoches au stator, et de 18 paires de pôles. Les 36 aimants sont régulièrement placés sans interstices. Cette topologie donne presque uniquement les ordres (0, 36, -36) dans le spectre des forces : cela permet une comparaison fiable avec le PTV en l'absence d'ordres qui parasitent le résultat. D'après l'effet de modulation de la Section 2.2, on voit que les ordres de force (36 , -36) sont susceptibles de se replier sur l'ordre 0. Grâce aux symétries, la simulation EF est réduite à une seule encoche : cela permet de symétriser le maillage afin de réduire le bruit numérique sur la simulation EF mécanique. Pour la même raison, on utilise 306000 éléments triangulaire au total. De plus l'utilisation d'une bande de glissement entre le rotor et le stator permet de conserver le même maillage à tous les pas de temps.

Le résultat de la comparaison est présenté sur la Fig. 2. L'accélération du TM est obtenue en simulant seulement le nombre d'onde 0, et en reconstruisant la contribution des ordres 36/-36 grâce aux formules (4) et [START_REF] Wang | Analysis of Zeroth-Mode Slot Frequency Vibration of Integer Slot Permanent-Magnet Synchronous Motors[END_REF]. On observe que la méthode TM-FANG donne un résultat beaucoup plus proche de la référence calculée avec le PTV par rapport à la méthode TM-WANG, en particulier pour la raie à 3564 Hz. La méthode de TM-FANG permet d'obtenir une résultat similaire au PTV avec une erreur inférieure à 5dB (Ref :1e-6 m/ 2 ) jusqu'à 6000 Hz.

CONCLUSIONS

Lors de cette étude, une méthodologie de validation pour le couplage magnéto-mécanique a été proposée dans le cadre de la simulation numérique du bruit et des vibrations dans les machines électriques. Les résultats de cette étude permettent de mieux expliquer l'intérêt du tenseur de Maxwell. Un cas particulier permettant une comparaison entre un chargement dentaire classique du tenseur de Maxwell et un chargement local calculé avec le principe des travaux virtuels a été réalisé. En particulier, l'effet de modulation et les différentes méthodes pour le prendre en compte furent étudiés. L'une des deux méthode étudiées est clairement plus précise pour l'application étudiée. La version étendue de ce papier détaillera les paramètres de simulation de la machines et proposera d'améliorer les différences pour prendre en compte les nombres d'ondes plus élevés à l'aide de l'effet de modulation.

RÉFÉRENCES

Fig. 1 .

 1 Fig. 1. Forces nodales calculées avec le PTV à 720 Hz.

Fig. 2 .

 2 Fig.2. Amplitude des déplacements (RMS) en fonction de la fréquence.

  D'après les travaux récents, la corrélation entre le TM et le PTV semble bonne pour les nombres d'onde et harmoniques les plus faibles [1, 3]. Cependant, peu d'études comparent directement les résultats vibratoires du PTV avec le TM pour les vibrations d'origine EM. En effet, le résultat du PTV donne une résultante calculée par rapport au maillage magnétique. Dans ces conditions la projection sur un maillage mécanique n'est pas immédiate, alors que les maillages magnétiques ne sont pas toujours adaptés pour calculer un résultat EF vibratoire précis. De plus, l'algorithme du PTV doit être adapté pour chaque type d'éléments et de résolution, contrairement au TM.


Rai -Amsterdam

Eddy Renaud

Alain Husson

Nicolas Rodriguez

Marie Collet

Vincent Couallier

PI to take into account water loss reduction in the SWDE IAM strategy

OUTLINE

METHOD -PIPE WATER LOSS PERFORMANCE INDICATOR: PIWL

 PIWL is designed to prioritize pipe sections that have a high probability of leaking, that are part of a District Metered Area (DMA) with high level of water losses, and for which Active Leakage Control (ALC) has limited effectiveness. This priority is also set higher if the DMA presents risks of water shortage. Its formula is the following:

 Where  PLP is the pipe leakage probability (between 0 and 1)  ∆LLI is a fonction of the difference between the yearly average Linear Losses Index (LLI) and its reference value for the DMA (m 3 /km/day)  An evaluation protocol remains to be put in place to assess the effectiveness of the approach and, if necessary, adjust the method.

 NFDI Night Flow Drift Index,

RESULTS -PLP: CALIBRATION

  performance indicator • SWDE (the Walloon Water Company), manages a distribution pipe network 27,000 km long divided into 2,000 DMAs. INRAE is a French public research organization • SWDE and INRAE have teamed up to conduct an applied research project entitled 'multi-scale asset management of drinking water networks' (GePaME, 2020-2023). • One of operational results of the project is a multi-criteria decision support tool for pipes renewal based on Performance Indicators



  is a function of the slope of the statistical distribution of daily LLIs of the DMA  WSI Water Shortage Index, is a function of the number of water shortages recorded for the DMA in the past Available information: Leakage state of a pipe (Z = 0/1) during a Leak Detection Campaign (LDC) within a DMA, knowing that some of the pipes are not inspected.  Z = R x Y  Y = 0/1 : leak detection at the time t with covariate ∈ℝ  R = 0/1 : pipe inspection at the time t with covariate ∈ℝ  As R and Y follow a logistic law, we can write : The municipalities of SWDE are classified according to 3 categories :  Cat1: Restrictions on use have occurred in the last three years  Cat2: Temporary measures have been put in place over the last three years (filling of tanks by trucks) but without restrictions on use  Cat3: No restrictions or temporary measures in the past  For each DMA, L1, L2 and L3 are the lengths of pipes located respectively in municipalities of categories 1, 2 and 3 WSI varies between 1 (DMA entirely located in a category 3 municipality) and 2 (DMA entirely located in a category 1 municipality)

PIWL COMPONENTS FOR THE TOP 1% OF LENGTH

  As part of the GePaME a water loss indicator calculated at the pipe scale was developed to feed a multi-criteria decision support tool for pipes renewal with the aim of limiting water losses and water shortages  A new model, giving promising results, has been developed to estimate the probability of pipe leaks  It was combined with indicators calculated at the DMA scale to try to maximize water savings and minimize water shortages while targeting DMAs where ALC is the least effective  The multi-criteria tool will be implemented in 2024 by SWDE.

	RESULTS -STATISTICAL DISTRIBUTIONS
	Calibration period: 01/01/2018 -31/12/2021 Validation period: 01/01/2022 -31/12/2022 OF OTHER COMPONENTS
	"Backward" selection: 18 models tested
	Full model: 
	 23 covariates, 69 parameters to estimate	 Material (6 levels)
	 20 covariates, 47 parameters to estimate  Area under the curve Best model:  4 covariates non significant  AIC= 62,468 ∆LLI	 Pipe length (meter, log-scale)  Number of connections (log-scale)  Diameter (decimeter) NDFI
	 Geographical sector (14 levels)  Urbanisation  Roads' type in 2022 WSI leaks actually repaired  AIC=62,400 0.84  Inspection state R: 6  The 20% of the  Material (2 levels) leakiest pipe sections  Length of the DMA (meter, log-scale)  Number of past leaks or bursts account for 70% of the 5.6 according to the model 0.36	 Age at inspection (century)  Number of past pipes' failures  Pressure (millibars)  Proportion of pipes under the road  Urbanisation  Loam proportion in the soil  Geographical sector (9 levels) ∆LLI x NDFI x WSI  Clay proportion in the soil
	 Year of the Leak Detection
	Campaign (LDC)
	 Duration of the LDC

 Leakage state Y (PLP): 13 STATISTICAL DISTRIBUTION OF PIWL 5.6 14 CONTRIBUTIONS OF


S André

) S Becker

L Farge

A Delconte

Certifying the intrinsic character of a constitutive law for semicrystalline polymers: a probation test

Keywords: Inverse identification, parameter estimation, constitutive behavior's law, semi-crystalline polymer, tensile test, HDPE

A study of methodological nature demonstrates the efficiency of a probation test allowing for the intrinsic character of a rheological constitutive law to be assessed. Such a law is considered here for Semi-Crystalline Polymers exhibiting necking and for large deformation. In the framework of a (𝜎, 𝜎, 𝜀, 𝜀) behavior's law, tensile experiments conducted at an imposed constant strain rate 𝜀0 bring true stress responses from which constitutive (material) parameters can be identified from Model-Based Metrology concepts. The same experiment repeated at various strain rates gives then access to the dependence of the nonelastic parameters on the strain rate. Then the intrinsic law is tested severely by considering a new set of experiments carried out for constant displacement rates of the grips. In that case, the specimens show local strain rates which evolve strongly during the test (by a factor of 5-10 here). The parameter identification process requires then the introduction of the exact realized input strain and strain-rate command into the model. Accounting for strain rate dependency requires additionally the knowledge of the preliminary identified strain rate dependence of the non-elastic constitutive parameters for good predictions of the experimental response directly. This is what is proven here. The conclusion speaks in favor of a possible upgrade of international standards for the mechanical characterization of polymers based on constant strain-rate tensile tests and properly applied model-based metrology.

Introduction

Model-Based Metrology (MBM) means essentially a whole of procedures, at the heart of which are the principles of parameter estimation theory (identifiability proofs) to treat data with the objective of obtaining measurable quantities of physical interest. Such an approach should be a rule of thumb each time a model is used in conjunction with observed quantities because (i) it allows stating in full transparency the performance of the model and the quality of the estimates, but, most important, (ii) it gives a feedback of possible missing foundations of a model or theory. These are revealed by the shape of after-identification residuals (the gap curve between observed and model-reconstructed signals). Such practices are very common in other fields of science (especially in the heat transfer community which was at the origin of the International Conference on Inverse Problems in Engineering ICIPE, but also of course in the automatic and system identification communities) and many textbooks on this subject are now available [START_REF] Beck | Parameter Estimation in Engineering and Science[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF][START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF] to cite a few references). But it is pretty much absent from mechanical/rheological studies. In this field, an example of the strength offered by these practices can be studied in [START_REF] Maillet | Regularized parameter estimation through iterative rescaling (PETIR): an alternative to Levenberg-Marquardt's algorithm[END_REF] where it is shown how two models equally optimal for experimental data description can be discriminated through mathematical arguments of parameter identifiability. In the paper of [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF], these latter were used to characterize the behavior's law of a semi-crystalline polymer (SCP). It was shown in particular that the instantaneous elastic modulus can be retrieved from the common tensile test with a very good agreement with the values determined more directly by other techniques operating at microscopic level. The present paper aims at giving now the second step of the methodology. It consists in proving the validity of a given constitutive law model with a probation test based on its underlying physical concepts: if a rheological relationship (𝜎̇, 𝜎, 𝜀, 𝜀, 𝛽) is really intrinsic to some material then, any loading path producing strain and strain rate fields heterogeneities in the specimen structure should lead to the same set of model parameters 𝛽. In section 2, information regarding the constitutive model we used and the numerical computation associated within the identification process is given. Section 3 describes the three test cases used for this study. Section 4 discusses aspects of the identifiability of material parameters in the most convenient case where tensile experiments are conducted with the specimen center constrained to follow a constant strain rate path. It additionally produces a so-called 'companion' relation of the constitutive law, in charge of describing strain rate effects on the relaxation spectrum. Section 5 presents the direct computations and identification results that were obtained from a probation test: tensile experiments performed now at imposed displacement rate, where strains and strain rates evolve now in a totally different manner with respect to time.

Constitutive model and numerical treatment

2.1.Operational set of equations

The thermodynamical foundations of the model have been given elsewhere [START_REF] Cunat | A thermodynamic theory of relaxation based on a distribution of non-linear processes[END_REF](Cunat 2001a) ) and we focus directly on the mathematical operational form which consists in a set of modal Ordinary Differential Equations (ODE). Each mode is denoted by the 𝑗 subscript in equation ( 1) and the overall stress response 𝜎(𝑡) is obtained by summation over all modes, meaning that each mode contributes to the global response. Each modal ODE contains the simple description of a relaxation process through a first order kinetic model of characteristic/modal time 𝜏 𝑗 . This can be seen in a first approach as a set of analogical Voigt-Maxwell units which would be connected in parallel. These units are uncoupled one with each other (in the simplest view), but the weight 𝑝 𝑗 0 of a single unit on the global response is connected to the modal time 𝜏 𝑗 through a single "universal" law. As a result and for a given spectrum distribution, a hierarchically recursive scheme applies to the modal weights. More details on this are available in Andre et al. (2003), where the connection of this approach to non-integer differential constitutive laws is explained (see for example [START_REF] Lion | On the Thermodynamics of fractional damping elements[END_REF].

r jj u jj j E  -  =  -  ; N j j1 =  =   (1)
where 𝜀(𝑡) is the strain rate, 𝜎(𝑡) is the overall stress, 𝜎 𝑗 (𝑡) its j th modal contribution.

Regarding the instantaneous or unrelaxed contribution to 𝜎(𝑡), 𝐸 𝑗 𝑢 is defined as

u j u j E p E 0 =
where u E stands for the common elastic modulus (Young modulus denomination will be avoided to favor the thermodynamic foundation of this parameter as the twice differentiated thermodynamic potential with respect to the strain variable) and 0 j p denotes the modal weighting coefficients fulfilling the normalized condition:

N 0 j j1 p1 = =  (2)
Regarding the relaxed (dissipative) contribution to 𝜎(𝑡), we define the modal relaxed stress  refers to the stress in the relaxed state, a thermodynamic state which has been properly defined by Prigogine and Defay (1958). It corresponds to the stationary state for non-equilibrium internal forces and is defined by 𝐴 ̇ 𝑟 = 0 where 𝐴 denotes the affinity in the T.I.P. theory [START_REF] De Donder | Thermodynamic theory of affinity: A book of principle[END_REF]. It implies formally a direct coupling between internal variables reflecting the microstructure evolution and the macroscopic command (strain rate here). Logically, this state can be identified when all viscoelastic modes get through, which we associate to the set-up and further evolution of the fibrillar state in SCP. It can then be simply and very efficiently modeled here with the

hardening law r HT (t) G (t)  =  (3)
suggested by theoretical considerations on the elastomeric state (Treolar 1975;[START_REF] Haward | Strain Hardening of Thermoplastics[END_REF][START_REF] Haward | Strain hardening of High Density Polyethylene[END_REF][START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF][START_REF] Haward | Strain Hardening of Thermoplastics[END_REF][START_REF] Tervoort | Strain-hardening behavior of polycarbonate in the glassy state[END_REF] and experimental evidences (Supplementary Information Fig. S1). The parameter G (called the hardening modulus) introduces the proportionality of the relaxed stress with the 'Haward-Thackray' strain variable exp( ) exp( ) Generally, 𝑁 = 50 dissipative modes are considered to figure out a continuous spectrum.

2 1 HT 2      - = - = -
A number of decades 𝑑 = 6 is generally required for the tested polymers so that the spectrum extends to sufficiently low relaxation times beyond which the model becomes unsensitive.

The spectrum is then determined once the single parameter max  , the longest relaxation time, is fixed.

The key conceptual tool is to connect the dissipative modal weight 0 j p to its corresponding relaxation time j  and this is made with a thermodynamically based argument stating an equipartition of the entropy created by each dissipative (relaxation) mode [START_REF] Cunat | A thermodynamic theory of relaxation based on a distribution of non-linear processes[END_REF][START_REF] Cunat | The DNLR approach and relaxation phenomena : part I : Historical account and DNLR formalism[END_REF][START_REF] Faccio-Toussaint | Modeling of the Mechanical Behavior of a Nickel alloy by Using a Time-Dependent Thermodynamic Approach to Relaxations of Continuous Media[END_REF]. This leads to (1) in terms of a recursive series to an alternative formulation in terms of derivatives of fractional order (exact in the limit of an infinite number of relaxation times). All details about this interesting link between both formulations are extensively described in Andre et al. (2003).

The constitutive relationship expressed through equation ( 1) can be rewritten now as

0 N N j j HT 0u j j j 1 j 1 j p G p E = =    -   =  =  -         (6)
The functional form associated to the behavior's model of Eq.( 6) can be written as =   . The relevance of this 3-parameter model has been proven in Blaise et al. (2016). The sensitivity analysis and parameter identification procedures have shown how well conditioned is this model with respect to the Parameter Estimation Problem (PEP). The analysis of the after-identification residuals proved that the viscoelastic behavior is correctly caught by the modal recursive approach, as well as the hardening stage within the approach of Eq.(3). Evidence of this will appear again when presenting identification results in Section 4.

2.2.Numerical solution

When we consider a tensile test with a constant imposed strain rate 𝜀̇0, namely 𝜀(𝑡) = 𝜀̇0, eq. ( 6) can be solved analytically for each modal contribution as:

( ) ( ) ( ) ( ) ( )( ) j 00 j t t 2 t 0 0 0 j 0 j 0 j j j 0 t 0u j 0 j j j 0 j 0 0 j 3p e p 2 1 e p 1 e t p E 1 e G 1 2 1 / / - -  -    -  +  +  -    =   - +      -  +    for each modal "branch" and N j j1 (t) (t) =  = 
 for the overall response (7)

Tensile tests at constant imposed strain rate are, for this major reason but among other ones, a clear advantage when metrology is the foremost objective.

In the case considered later where the ever-changing true strain rate excitation must be introduced as the input of the model, a numerical scheme is required that computes the solution as time steps proceed. We use then a classical two steps implicit Euler Backward second order scheme (GEAR or bE2) applied on the system of ODE's of (Eq.6). For the 𝒋 𝒕𝒉 modal component, we consider ( )

r j j 0u j j j j j σ (t) σ (t) t,σ σ (t) p E ε(t) τ τ f = = - + + (8) discretized through ( ) n+1 n n-1 n+1 n+1 j j j j 4 1 2h σ σ σ t, σ 3 3 3 f = - + (9)
with the time step

n+1 n h t t =- .
Initialization of the algorithm rests upon the two first time steps: the initial condition 𝜎 𝑗 0 (material in equilibrium state) and a first single step Euler Backward scheme of order 1. An alternative scheme has been used, very precise too, based on a semi-analytical approach of the integration of the system of eqs. ( 6). This algorithm (SA3) is described in [START_REF] Sorvari | Time integration in linear viscoelasticitya comparative study[END_REF]. It is used in Section 4 to ensure proper implementations of both algorithms but it extends pretty much the computation times. Furthermore, in the case of a non-idealized forcing term (experimental data), it amplifies the input noise and appears therefore not appropriate to this study.

Considered test cases

Three different experiments will be considered for this analysis. They all apply to specimens of HDPE (same original material from Röchling Grade Natural) at same room temperature. Reproducibility is not discussed as already proved excellent with our metrological equipments (elements provided for example in [START_REF] Ye | Kinematic study of necking in a semi-crystalline polymer through 3D digital image correlation[END_REF]. Information about the metrological tools involved in these experiments are not discussed here but also available from the above-mentioned reference.

Experiment 1 refers to a tensile test performed at an imposed constant strain rate 1 0.005 s

 - =
. It relies on a VideoExtensometer working with a limited number of markers located in the specimen gauge (4/5/4 dots located on 3 different longitudinal lines), providing a real-time measurement of the local longitudinal strain. This signal drives the actuator of a hydraulic tensile machine through a PID feedback loop. In the case of a ramp input; the realized strain signal perfectly follows the input command as can be seen in Fig. 1: the constant strain-rate is achieved with a small random noise of 0-mean superimposed to it.

Experiment 2 refers to a tensile test performed at an imposed constant displacement rate

0.02 / u mm s =
on the same machine.

Experiment 3 is identical to experiment 2 except for an imposed constant displacement rate now 4 times greater: 0.08 / u mm s = .

For all 3 experiments, a 3D (stereo) DIC system provides full field strain measurements on the 2 front and lateral faces of the specimen gauge, allowing for precise true stress measurements [START_REF] Farge | WAXS study of the lamellar/fibrillar transition for a semi-crystalline polymer deformed in tension in relation with the evolution of volume strain[END_REF]. These 3 experiments will serve the objective of the demonstration. With respect to an intrinsic behavior's law expressed in terms of true stress and strain, they provide very different strain and strain rate excitations on the Representative Element Volume (REV) of the material as shown in Figure 1.

Figure 1: Strain-rate evolution with strain during a tensile test performed at constant imposed strain or displacement rates.

It is of the utmost importance at this point to insist on the fact that only true strain-true stress analysis of the identification results is relevant to certify the intrinsic character of the behavior's law. As a consequence, the case of experiment 1 is more comfortable for the scientist: experimental data are obtained directly with the variables suitable to launch the inverse approach based on the model of Eq.( 6). This case will be specifically studied in section 4. Note that true strain-true stress curves are in that case identical to time-true stress curves (the shift in abscissa is the proportional factor 𝜀). The case of experiments 2 or 3 requires from the scientist a preliminary work of determining the local evolution of 𝜀(𝑡) before going on with the inverse approach (detailed in the Supplementary Information file).

Additionally, as evidenced in Fig. 1, the strong evolution of 𝜀(𝑡) and 𝜀(𝑡) requires to know or to identify the relationship between the relaxation parameters of the behavior's law and the strain rate. This will be the object of section 5. Note that in that case, the only analysis in terms of true strain-true stress makes sense for model identification, hence the adopted representation of the results.

Experiment 1 (constant strain rate): Identification of model parameters.

Applying the model (Eq. 6) to the experimental data obtained for experiment 1 through an inverse estimation of the model parameters will give key elements of discussion to assess our demonstration. Estimation is performed on a least-square criterion optimization and using indifferently a Levenberg-Marquardt or simplex algorithm. Because the problem is relatively well-conditioned (see explanations in [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF] results are always identical.

Accounting for the Levenberg-Marquardt algorithm's fast computation times, this latter will be used preferentially.

In figure 2 below, we show a typical result of the adjustment of the model to the data for experiment 1. We selected a special set of data which, unlike the one shown in Fig. 1, suffers from accidental superimposed low-frequency noise on the input (due to a different set of PID settings). This signature will make clearer our analysis of the residuals obtained after the optimization process. Four different identifications were applied on these data whether the true or idealized strain and strain-rate were considered in the model or the numerical or analytical versions of the model were used (Eq. ( 7) versus eq. ( 9)).

• Identification 1 = analytical model of Eq.7, which implies the assumption of a perfect ramp command effectively realized in the test.

• Identification 2 = Numerical algorithm (bE2) which uses the real (noisy) strain and strain-rate input signals. Considered as the reference case for the study.

• Identification 3 = Numerical algorithm (bE2) which uses "filtered" (un-noisy) strain and strain-rate input signals. Typically, a linear regression is applied to the measured input strain command and is substituted for it in the algorithm.

• Identification 4 = Same as Identification 3 but using the (SA3) algorithm to prove the numerical consistency.

As shown by Figure 2, the least-squares optimization of the model produces absolutely identical results although the identified parameter values can show 2-3% of variations (see Table 1). The model can explain the data very well: the after-identification residuals -i.e. discrepancy between experimental and reconstructed data-produce a control of the quality of the optimization process. In the residuals signature and as a general rule of thumb, two components are present that can be recognized easily in Fig. 2. and difference between input and realized strain ramp command (arbitrarily amplified and shifted to a mean value set at +40).

• The noise at higher frequency is a direct contribution of the initial noise pre-existing on the input excitation. In general, it is essentially the intrinsic noise of the sensors and possible electronic conditioning of their raw signal. Here there is an additional nuisance noise due to the input command as mentioned previously. Estimation parameter theory shows how the noise affects the variance on the estimated parameters. But in case the model is a perfect idealization of the experiment and the sensitivity analysis shows highly sensitive and perfectly un-correlated parameters, it does not bias the estimated parameters. Consequently, it is fully recoverable in the identification residuals i.e. the s.t.d of the residuals should be equal to the s.t.d. of the noise on the data. This is nearly the case here. After-identification residuals (right axis) are nearly centered on or fluctuate around a 0-mean level and corrupted by a noise signal which corresponds exactly to the 0-mean input noise observed on the strain signal. This can be seen from the data set shifted to the arbitrary 40 mean value for clarity of Fig. 2 (second curve, same right axis) which indeed corresponds to the direct difference between the input measured strain and its expected behavior (the command).

• The long-range wavy behavior of the residuals is reflecting the sole bias due to some deficiency in the model/experiment agreement. This bias either results from an uncontrolled bias in sensor measurements (clearly not the case here) but in contemporary science is more likely the result of a more or less pronounced incompatibility between the assumptions of the model (an ideal design) and the objective experimental conditions. From our insight on this specific problem, it is more likely due to the approximate phenomenological description of the relaxed state according to eq.( 3). A bias due to the force signal conversion into a true stress signal, based on DIC surface measurement strains on the three surfaces of the specimen in its central part cannot be excluded totally either. We indicate on Fig. 2 where the local bias is at maximum (of the order of 1.5MPa at a strain of 1.2).

ˆu

E (MPa) max T  (s) Ĝ (MPa)
Relative error

Relative error

Relative error Identification_1 Of course, it is essential to constraint the bias as low as possible. In general, the variance of the noise on the true strain signal is far below the bias, and is of the order of 0.04 𝑀𝑝𝑎 (in our experiments, which corresponds to a noise on the Force sensor of s.t.d less than 2N).

In Table 1, the identification 2 using bE2 algorithm is taken as the reference as a numerical approach is mandatory to treat data obtained in the case of experiments of type 2 and 3.

5. Experiments at constant displacement rate: the probation test.

Because a tensile test performed at constant displacement rate makes the central part of the specimen (where the true stress-true strain curve is recorded) subjected to very different strain rates (Fig. 1) due to necking, the constitutive model which is expressed in terms of a (𝜎̇, 𝜎, 𝜀, 𝜀) law must obviously consider the real excitation 𝜀̇ but not only. The key point is that under evolving strain rate excitations, the specimen's viscoelastic response reflects a change in relaxation times. As a result, the condition for the model to be applied to fit any kind of experimental excitations with pertinent results, is to have available a relationship between the strain rate and -in the framework of our model-the maximum identified relaxation time. This latter was obtained in Blaise et al. [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF], the constitutive law having been applied to a series of stress-strain curves monitored at different constant strain rates. Thanks to the well designed and parsimonious model described in section 2 through Eq.( 7), the identified parameters show a strain rate dependence which is expected to be a real intrinsic property of the material. For the seek of some universal character, the Weissenberg number We was introduced. It corresponds to the ratio between the time constant mat t characterizing the intrinsic 'fluidity' of the material and the time scale of the experiment or of the observer exp t . The fluidity of the material is inversely proportional to the Weissenberg number. In the present case, its maximal value can be calculated by the where max  , the maximum relaxation time of the spectrum, is used for the material characteristic time and 1  for the experimental characteristic excitation time. We and hardening moduli G for the HDPE specimen subjected to tensile tests performed at different strain rates.

5.1.Strain rate dependency of model parameters

• Time spectrum: There are two options regarding the "modeling" of the relaxation times spectrum dependency with respect to strain rate. For such dependence description, we will speak of "companion" relation of the constitutive law, rather than model, law or rule. Indeed, it can not be considered as part of the theoretical modelling. It is different for instance as to consider models for describing plasticity and damage, including the definition of a yield criterion… to complete the formal description of the whole behavior In this case, this latter is built from a set of equations (laws) catching independently each physical mechanisms: Norton's law for the viscous description, a hardening law, role played by eq. 3 in our approach, the GTN model for plasticity and damage phenomena… Such laws are either phenomenologically or theoretically founded and examples of their use to model the behavior or to interpret experimental data of polymeric materials are numerous [START_REF] Laiarinandrasana | Temperature dependent mechanical behaviour of PVDF: Experiments and numerical modelling[END_REF][START_REF] Pardoen | An extended model for void growth and coalescence[END_REF], Negahban et al. 2006).

Option 1 As reported in Table 1, experiments conducted at various imposed strain rates have provided the identification of different maximum relaxation times. The corresponding calculated Weissenberg number varies from more than 50% over 3 decades of strain rates but in an unknown manner (no model available). A plot of We as function of  in logarithmical scale shows a roughly linear dependency. Applying a regression, gives the following relationship with R-square coefficient of 0.976. The spectrum parameter max

We

 can be considered as known (simple division of We by  ) once the local strain rate  imposed locally on the REV is known.

Two controls can be applied to check whether this relationship is valid for our material.

Either the adjustment procedure is reduced to the identification of the only one remaining parameter (Relaxed modulus G ) because the model parameter max  is supposed to be known from (11) (𝜏 𝑚𝑎𝑥 𝑊𝑒 (𝜀));

Or parameter max  will be considered as corresponding to max max We  =  and the adjustment procedure will consider both the relaxed modulus G and the correction coefficient  as parameters to identify. The latter will measure in some ways, whether or not the residuals can be further minimized when authorizing max  to slightly fluctuate around the supposed known value max We  . Of course, if the relationship (11) perfectly holds, this coefficient should be identified equal to 1. As this relation results from measurements also delivered by an identification process, some errors may exist but if rather insensitive to the obtained results, this will prove the intrinsic character of both the constitutive law and its companion material spectrum relation.

Option 2 Without any explanatory concern, a direct fit can be applied to the 𝜏 𝑚𝑎𝑥 (𝜀) data of Table 2. In that case a perfect representation of the data is obtained (R-square coefficient = 1) with the following exponential relation:

𝜏 𝑚𝑎𝑥 𝑓𝑖𝑡 (𝜀) = 27.24 𝑒𝑥𝑝(-2.1 𝑙𝑜𝑔10 3 𝜀) (12)
The 𝜏 𝑚𝑎𝑥 𝑊𝑒 (𝜀), 𝜏 𝑚𝑎𝑥 𝑓𝑖𝑡 (𝜀) and 𝜏 𝑚𝑎𝑥 (𝜀) data are plotted in Figure S2 (Supplementary Information). It is shown that both relations ( 11) and ( 12) capture the apparent strain-rate dependence even if the 𝜏 𝑚𝑎𝑥 𝑊𝑒 (𝜀) approximation underestimates the time spectrum higher bound at very low strain-rates.

• Hardening modulus: In the same way, the Hardening modulus identified from the stressstrain response to an input ramp at constant strain rate evidently depends on the selected strain rate. This is generally observed indirectly by varying the temperature of polymers.

The decrease of 𝐺 when the temperature increases was already observed for both amorphous (Van Melik et al., 2003) and semi-crystalline polymers [START_REF] Na | Inverse temperature dependence of strain hardening in ultrahigh molecular weight polyethylene: Role of lamellar coupling and entanglement density[END_REF]. It is due to thermally activated relaxation mechanisms that can be evidenced by DMA (Dynamic Mechanical Analysis) measurements on pre-deformed specimens. In the case of polyethylene, the identified relaxation process is the 𝛼 relaxation mode associated with the slip of mosaic blocks inside the lamellae [START_REF] Na | Inverse temperature dependence of strain hardening in ultrahigh molecular weight polyethylene: Role of lamellar coupling and entanglement density[END_REF]. Because of the temperaturetime equivalence principle, a strain rate increase has the same qualitative effect as a temperature decreases and results in an increase of 𝐺 values which is exactly what we observe from parameter 𝐺 identifications from experimental tensile curves. But as far as the authors know, there is no reported idea about the mathematical description of this dependency. Unlike relaxation times which can be assumed to shift rather simultaneously with a change in microstructure, a direct instantaneous dependence on strain rate could hardly be defend for the relaxed modulus G . It would be firstly in contradiction with the hypothesis of describing precisely a relaxed state. Secondly, as it describes the elasticity of the SCP's fibrillar network established differently according to the loading path history, it can change physically only in a delayed manner (like for temperature). For those two reasons, we will keep the strict relation (Eq. 3) for the relaxed state, the modulus G representing the apparent hardening elasticity in a somewhat homogenized manner for constant displacement rate tensile tests.

Overcoming the probation test

The objective of the paper is achieved in this section. Figure 3 shows the best model The model accounts for the exact realized input strain command with varying strain rates of Fig. 1 (see explanations of the fitting process of these command signals in Supplementary Information: Curve fitting process for the command/excitation signals). But it does not account yet for the strain rate dependency evidenced in the previous section on the max  model parameter (Eqs.11 or 12). The demonstration is clear. The constitutive law model is unable to reproduce the real behavior which nevertheless has an apparent behavior identical to experiments monitored by strain rate (compare with the stress-strain curve obtained at constant strain rate in Fig. 2). It produces a short time behavior corresponding to a sudden jump in stress with very low yield stress (around 2MPa) and the post-yield behavior is simply unacceptable as well as the identified instantaneous modulus found to be around 60 GPa. An expeditious conclusion would lead to consider the proposed constitutive law as non-intrinsic to the material, which is, of course, what a rheologist tries himself to prevent from. 12) is used to subjugate parameter max  of the constitutive behavior's law (Eq. 6) to the actual strain rate, then the prediction of the true stress-true strain curves for both experiments at constant displacement rate (Fig. 4a when using Eq.11 and Fig4b when using Eq.12) is very reliable with all rheological regimes of the tensile curve clearly recovered (Comparison to be made with curves of Fig. 3). These direct simulations were performed for an unchanged instantaneous modulus of 3220𝑀𝑃𝑎 (Table 1-Line 2). Concerning the hardening modulus G, the best strategy is to determine it in accordance with the theory used for the relaxed state description as was presented in [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF]. The plot of the experimental stress versus the Haward-Thackray strain variable should exhibit a linear behavior. This is nearly the case for both Experiment 2 and Experiment 3 tensile tests. From a linear regression (see Fig. S1 in Supplementary

Information) it can be easily retrieved as G=1.56 MPa and G=1.64 MPa respectively for the experiments at 𝑢̇= 0.02 𝑚𝑚/𝑠 and 𝑢̇= 0.08 𝑚𝑚/𝑠 . It is clear from this direct modelling that the constitutive model (Eq. 6) is able to describe the experiment in an intrinsic manner. Residuals show of course that the bias is now enhanced as a consequence of the approximation of Eqs.11 or 12 to determine the spectrum shift with strain rate variations and of the a priori determination of the overall hyper-elastic modulus G . The extreme sensibility of the spectrum shift is evident from Figs 4a and 4b. Although the discrepancies between the model prediction and the experiment remain at maximum equal to 10 MPa, it happens rather in the yielding regime (Fig. 4b) or in the hardening regime (Fig. 4a) whether Eq.12 or Eq.11 are used respectively.

This suggests that the metrological approach based on a set of tensile tests at constant strain rate should be conducted very carefully.

5.3.Identification of the model with strain rate correction

The analysis can be pushed a little forward by producing model adjustments on the experiments, with the same constitutive model of eq.6, the exact realized input command signal but allowing for the adjustment of parameters max  and G .

Results of identifications are reported in Table 3 for experiment 2 (𝑢̇= 0.02 𝑚𝑚/𝑠).

Again and accounting for the results of Table 1, the instantaneous modulus measured in tensile tests of the type of experiment 1 is considered as known and equal to 𝐸 𝑢 = 3.22 𝐺𝑃𝑎 ('Identifications 2 and 3' in Table 3). For 'Identification 2', the parameter max  of the relaxation spectrum will be determined as either  for providing a better fit. Parameter G will be considered as constant to provide the best overall estimation of the behavior in hardening regime (although it has been observed, as discussed earlier, that the establishment of the fibrillar regime depends on the strain-rate). It can be initialized to the estimated value given in section 4.2. For 'Identification 4', both three parameters 𝐸 𝑢 , 𝜏 𝑚𝑎𝑥 , 𝐺 are identified to check the stability of the estimation process. Regarding identified parameters, it is shown (Identifications 2 and 3) that whether Eqs 11 or 12 are used for initial estimation of ( ) max   is not so important. The 𝛼 value identified to correct (shift uniformly) the relaxation spectrum is of the order of 10% at maximum. When the instantaneous modulus is free to adjust the data, it increases by 6.5% when the max  parameter can be less corrected (𝛼 = +4.5%). This is just the result of a slight correlation effect existing between these two parameters. Note that assuming no bias in the couple (experiment/model), the ideal relative error on these two parameters was calculated of the order of 2% in Blaise et al. [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF]. 'Identification 4' then proves the robustness of the P.E.P. conditioning. Regarding the 'apparent' hardening modulus G, nothing can be concluded except to observe that the initial value estimated from the linear regression at strong deformation is corrected by the estimation process performed on the whole strain range to product better adjustment to the experimental data.

Turning now to the tensile test performed at 𝑢̇= 0.08 𝑚𝑚/𝑠 (Experiment 3), the same analysis can be made and results of identifications are reported in Table 4. Results concerning the model adjustment, the residuals magnitude (lower than 3.4 MPa), and the parameter estimates are similar to Experiment 2 (see Table 4 and Figure 6) except for identification 4.

In this case, the instantaneous modulus is not considered as known but as a parameter to estimate. The less good-conditioning of the P.E.P. becomes obvious. Remember that P.E.P.

are Non-Linear Problems in terms of the parameters (and possibly also the inputs). It is clear here that the wide strain rate range in input introduces a stronger correlation between parameters 𝐸 𝑢 and  (i.e. max  ). These two parameters evolve now strongly (+36% and - 24%) and the model adjusts much further than for Identifications 2 and 3: the cost function is reduced by a factor of 25%. This can be seen on Fig. 6 where in that case, the residuals Concerning Identifications 2 and 3, the accommodation required on max  through  is much weaker in the case where the Weissenberg companion relation is used (+4.5% against -14% for eq. 12). Because in that case the strain rate variations are very important (up to 1.7 10 -2 𝑠 -1 𝑖𝑛 𝐹𝑖𝑔. 1) this may suggest that the Weissenberg approach (eq. 11) is more sounded that the mathematical fit (eq. 12).

Figure 6: Comparison between tensile tests at constant displacement rate and inverse identification of parameters 𝛼, 𝐺 with the model of eq. 6 and companion strain rate relations for correction of the relaxation spectrum (Eq.12).

Still, these last results confirm the intrinsic character of the constitutive proposed in Eq.6 for the HDPE semi-crystalline polymer, which is able to model tensile test in imposed displacement command mode once a complementary information is joint to the model to account for strain rate dependency. Note that this result tends to confirm two postulates about micromechanisms of deformation of HDPE: the relaxation times spectrum is directly coupled to strain rate, the hardening modulus characterizing the fibrillar structure is not so sensitive to strain rate and linked to it in a more complex manner, still to investigate.

Conclusion

The conclusion of this study is two folded. The first conclusion stays in the initial objective of the paper: illustrating how the proof can be made that a proposed constitutive model is really intrinsic to the material. Following previous results reported in [START_REF] Blaise | Advantages of a 3-parameter Reduced Constitutive Model for the Measurement of Polymers Elastic Modulus using Tensile Tests[END_REF], it is shown that a proper of theoretical concepts in P.E.P. allows the 'best as possible' estimation of the model parameters. These parameters, if really intrinsic to the material, should convey strain rate dependency if tensile tests are performed at different strain rates. Once the strain rate effect on the parameters is characterized and a thus named 'companion' relation is derived, the probation test consists in considering tensile tests performed now for varying strain rates all along the tensile test. Because heterogeneity in strain rates can result from the strain localization (or necking) appearing in SCP specimen, simple (usual) displacement rate monitored tensile tests are enough for the probation. The pair of constitutive behavior's model and companion relation describing the strain rate dependency is shown to be able to fairly well predict the material behaviors when the strain and strain-rates vary. Beyond this example, the method could be extended for instance to account for temperature dependence.

The second conclusion should possibly lead to better practices in mechanical characterization. The tensile test at constant strain rate, coupled to both a properly designed theoretical behavior's law and a deep investigation of the parameter identifiability problem is very effective to produce meaningful material parameter estimations. This could produce an evolution of international standards. To speak only about the instantaneous Young modulus, the ASTM and ISO standards for polymers still recommend to calculate the slope of a linear approximation of the stress-strain curve in a normalized strain interval. In polymers, viscoelastic effects start as soon as a mechanical excitation is imposed on the specimen and a modulus estimation should take this into account.



  is chosen as logarithmically distributed over a determined number of decades d below a maximum relaxation time max τ .

  2 where 𝑞 is a constant univocally determined by the parameters of the spectrum of relaxation times. This mathematical recursivity (corresponding physically to a self-similar property of mechanisms involved) makes it possible to relate the mathematical structure given in Eq.

Figure 2 :

 2 Figure 2: True stress-true strain curves for Experiment 1 -Left axis: Experimental data and adjusted model; Right axis: After-identification residuals (magnified by a factor of 4)

  with prior experiments and if really intrinsic to the material, this relation can be seen as the companion relation to be associated to the constitutive law.

  identification on the experimental stress-strain curves obtained for constant displacement rate tensile tests (Experiments of type 2 and 3).

Figure 3 :

 3 Figure 3: Experiments at constant displacement rate and model inconsistency when intrinsic strain-rate dependency is neglected

Figure 4a :

 4a Figure 4a: Comparison between tensile tests at constant displacement rate and direct model computation from eq.6 and companion strain rate dependency for the relaxation spectrum (Eq.11).

Figure 4b :

 4b Figure 4b: Comparison between tensile tests at constant displacement rate and direct model computation from eq.6 and companion strain rate dependency for the relaxation spectrum (Eq.12).

  will present the advantage of showing directly in % the amount of adjustment made from eqs.11 or 12 to max

  Figure5shows how well the model identification compares with the experiment for the 3 identification options. The estimated response of the model is the same, the Least-Square cost function leading to the same value at convergence. Because the noise on the experimental strain-stress signals is very low, the residuals (magnified by a factor of 5) are the pure reflect of the remaining bias. It is at maximum equal to 2.5 MPa which is pretty good (compare with Figs4a and 4b).

Figure 5 :

 5 Figure 5: Comparison between tensile tests at constant displacement rate and inverse identification of parameters 𝛼, 𝐺 with the model of eq. 6 and companion strain rate relations for correction of the relaxation spectrum (Eq.12).

  validate the whole process. This is what is shown in figure S8 below.

Figure S8 :

 S8 Figure S8 : Validation of the strain-rate fit with a numerical time integration to recover the strain signal.

Identification_2 Reference Numerical algorithm bE2 (Eqs 8 and 9)

		3112	-3.5%	5.394	+3.6%	2.552	-0.2%
	Analytical model
	(Eq.7)
		3224	-	5.207	-	2.556	-
	Identification 3	3215	-0.3%	5.22	-0.25%	2.55	0%
	Numerical
	algorithm bE2
	Identification_4	3218	-0.2%	5.22	-0.25%	2.557	0%
	Numerical
	algorithm SA3

Table 1 :

 1 Identified model parameters for the different options of model/algorithm for the identification.

Table 2

 2

	summarizes the findings for HDPE.

Table 4 :

 4 are very minimized at short times, in the early beginning of the tensile test (cross-dot residuals curve of Identification 4). The best strategy here is clearly to consider that parameter 𝐸 𝑢 remains at its known value, previously estimated and to consider only

	Identifications 2 and 3 as well-founded.
			ˆu E (MPa)			Ĝ (MPa)
	( ( (	Reference Identification_2 ( ) max max We  =   ) Cost function =2100 Identification 3 ( ) max max fit  =   ) Cost function=2250 Identification_4 ( ) max max We  =   )	3220 known 3220 known 3220 known 4374	Relative variatio n ---+36%	1 1.045 0.855 0.764	Relative variation +4.5% -14% -24%	1.64 (estimated) 1.85 1.86 1.85	Relative variation +12.6% +12.6% +12.6%
		Cost function =1665

Identifications performed on Experiment 3 (imposed constant displacement rate of 𝑢̇= 0.08 𝑚𝑚/𝑠)

Acknowledgements

The authors wish to acknowledge Jérémy Bianchin for technical assistance at LEMTA lab. The authors gratefully acknowledge financial support from different research programs involving the CPER 'SusChemProc' of Grand-Est Region, the CNRS 'PEPS' programs 2018 and 2019, the CARNOT label for its specific financial support #ICEEL-2019 and the EMPP Research Axis of Lorraine University.

constant strain rate (Experiment1) plotted as function of strain rate, along with the two relations (Eq.11 and Eq.12 in dashed and solid lines) proposed for their description.

Error bars correspond to a 10% uncertainty.

Curve fitting process for the command/excitation signal

When the excitation in the real experiment is applied on the grip displacement ( , ) uu , curve fitting of the excitation variables ( , )  of the behavior's law is required and was performed under Matlab Curve fitting Toolbox 3.5.6. Note that very good fits are sought to in order to make the model computations precise and un-biased but there is no need of any explicative model. Mathematical forms can then be selected without any physical constraints. For this reason the curve fitting process is performed in two steps to enhance the quality.

A first fit is performed of the signal () t  . In many cases, the mathematical function used for the fit is of the form A second fit is performed on these residuals. Generally a series in sine function of a few terms is obviously appropriate.

Figure S4 : Fit of the residuals under a series of sine functions.

Figure S5 plot the residuals after this second step. They are in this example of the order of 0.001 for a signal which varies from 0 to 1.8.

Figure S5 : fit of the 1 st step residuals and residuals after this second step.

Then Figure S6 shows the strain fit that can be obtained for the real experiment at u=0.02mm/s Figure S6 : Final curve fitting of the experimental strain signal.

The process can be repeated now on the experimental strain-rate signal 𝜀̇ which is obtained from a numerical differentiation. Figure S7 gives an example of the experimental strain-rate signal and its 1-step and final (2-step) fitted counterpart. Note that having a good un-noisy signal for the strain-rate realized during the experiment, it is possible to compute the strain signal resulting from the displacement command by


text:
Marius Dufraisse
email: marius.dufraisse@onera.fr

Pauline Trouvé-Peloux

Jean-Baptiste Volatier

Frédéric Champagnat

RGB Photometric Model for Optical/Neural Network Co-Design

We propose to relate a photometric model of imaging sensors to optical parameters, including sensor noise and non-linearities in optics/neural-network co-design. We apply it in extension of the depth of field examples.

Introduction

Imaging systems consist of a lens and processing, usually designed separately. Co-design methods aim at optimizing both components simultaneously. Joint optimization of neural network parameters and a phase plate on the diaphragm of a lens has been shown to improve depth of field extension, super-resolution, and depth estimation [START_REF] Haim | Depth estimation from a single image using deep learned phase coded mask[END_REF][START_REF] Sitzmann | End-to-end optimization of optics and image processing for achromatic extended depth of field and super-resolution imaging[END_REF]. Co-design can also optimize a complete lens and neural network [START_REF] Halé | End-to-end sensor and neural network design using differential ray tracing[END_REF][START_REF] Sun | End-to-end complex lens design with differentiate ray tracing[END_REF][START_REF] Tseng | Differentiable compound optics and processing pipeline optimization for end-to-end camera design[END_REF]. Optimizing simultaneously a neural network and optics relies on differentiable optical simulation based on Fourier optics or ray tracing. Generally, it simulates only the impact of optical parameters on aberrations but not their effects on light intensity.

However, some optical parameters, including the aperture, have a direct influence on the amount of light that passes through the lens: closing the diaphragm reduces aberrations but lowers the amount of light reaching the sensor at the expense of noise. The best possible compromise will depend on the processing's ability to denoise or deblur, so it is interesting to integrate aperture optimization into co-design to find its optimal value. This was first done by considering the aperture as a hyperparameter and testing several values [START_REF] Rostami | On design of hybrid diffractive optics for achromatic extended depth-of-field (EDoF) RGB imaging[END_REF]. To optimize it directly, we have already shown how to link Schottky noise to optical parameters [START_REF] Dufraisse | Deblur or denoise: the role of an aperture in lens and neural network co-design[END_REF] rather than simply adding a fixed variance Gaussian noise as most co-design works do.

We propose here an improved image formation pipeline that incorporates a more complete sensor model with multiple noise sources and non-linearities. It remains fully differentiable and as such can be used for co-design. We validate this pipeline with experiments on the extension of the depth of field (EDOF) and highlight the influence of the scene illumination on the co-design results.

Image formation pipeline

Our image formation pipeline first models the focusing of light by the lens before modifying the image intensity depending on the entrance pupil diameter and applying the sensor model. It is shown in Figure 1.

Optical model

The effect of a lens on an image is described by its Point Spread Function (PSF): the image I o formed by the lens is obtained by convolving the PSF with the sharp ground truth image I gt :

I o = I gt * PSF.
Several models exist for calculating the PSF and its gradient [START_REF] Dufraisse | On the use of differentiable optical models for lens and neural network co-design[END_REF]. Fourier optics is suitable for phase masks, but not for off-axis objects and multi-element lenses. Ray tracing solves these issues by describing a lens as a series of surfaces and using Snell's law to trace light rays through them. We use the Formidable1 library, which performs differentiable ray tracing (DRT) to compute ray trajectories and their jacobian with respect to the optical parameters. Using DRT we trace rays from a given point so that they uniformly sample the diaphragm. Following the path of these rays we obtain their intersection with the sensor plane which gives us a discrete approximation of the energy received by each pixel. To obtain the value of the PSF in each pixel and take into account diffraction, we integrate the convolution between the discrete distribution obtained by DRT and an Airy disk. Indeed, reducing the aperture reduces aberrations but can introduce diffraction. The Airy disk radius is computed in a differentiable manner using DRT tools.

These tools can compute the PSF for objects placed at different distances on and off the optical axis. They simulate how a scene would be projected on the sensor by the lens.

Radiometric sensor model

The PSF is normalized: it only takes into account the influence of optical parameters on the blur of the image but not on its dynamic range. However, some optical parameters, especially the aperture, have a direct influence on the amount of light passing through the lens: closing the aperture reduces aberrations but decreases the amount of light reaching the sensor at the expense of more noise. Our co-design pipeline models the effects of aperture, which is fixed in the literature, in two steps: an intensity change and a sensor simulation.

The intensity change is proportional to the square of the entrance pupil diameter which is the image of the physical aperture stop, as seen through the front of the lens. We use DRT to compute the image of the physical aperture stop and thus get the radius of the entrance pupil as well as its gradient with respect to the optical parameters. As such the intensity change is fully differentiable.

The sensor model adds noise sources and non-linearities. We first used a simple model which only adds Schottky noise to demonstrate the validity of our two step approach [START_REF] Dufraisse | Deblur or denoise: the role of an aperture in lens and neural network co-design[END_REF]. Now we simulate more noise sources (dark current and readout noise) and non-linearities (saturation of the photo-receptors and quantification of the signal) described in [START_REF] Tseng | Differentiable compound optics and processing pipeline optimization for end-to-end camera design[END_REF]. We also expanded the scope of our experiments with RGB images and off-axis objects.

Application to EDOF

We propose several experiments in which we apply our image formation model to an image quality enhancement task. Indeed, neural networks have the ability to deblur or denoise and thus, can have two distinct modes of operation when jointly optimized with aperture. Initially, we study a system with objects on the axis at a fixed distance to validate our model. Then we consider EDOF by varying the distance of objects. Finally, we include off-axis objects. We carried an ablation study to compare the importance of the multiple noise sources. In our experiments, the added sensor noise sources lead to similar optical system but the co-design process benefits from having 3 wavelength with 3 different focus distances to reconstruct the image. We also study how the optimized systems differ based on the level of illumination.

Conclusion

Our differentiable image formation pipeline incorporates the effects of lens parameters on the PSF but also on the image acquisition by the sensor. We applied it to the task of EDOF. We can now more accurately simulate the behavior of imaging systems and we plan to use our co-design methods to a range of imaging problems.

Fig. 1 .

 1 Fig. 1. Our image formation pipeline.

Freeform optics ray tracer with manufacturable imaging design capabilities: https://gitlab.space-codev.org/formidable