optimized regression algorithm. . . . .

List of Figures

Chapter 1

Curriculum Vitae • Reviewer of multiple international journals (International journal of material forming, Advanced modeling and simulation in engineering science...) • Organizing committee member to several international conferences and chairperson in multiple sessions • Program coordinator of the solid mechanics track and instructor of several courses in the same track • Advisory for more than 70 students • Member of the accreditation team with multiple sections of the self study reports written by myself • Recent awards:

-Most distinguished supervisor and best project in LIRA 2016 forum.

-Worldwide top 12 finalist in the ASME Innovative Design Simulation Challenge IDSC competition, 2017.

-Fulbright visiting scholar nomination by USA for the academic year 2020-2021.

• Visiting Scholar positions:

-Visiting Scholar to Ecole Centrale Nantes(France): September 27 -October 4, 2015.

-Visiting Scholar to Ecole Centrale Nantes (France): January 31 -February 9, 2016.

-Visiting Scholar to ETH-Zurich (Switzerland): November 21 -24, 2016.

-Visiting Scholar to Ecole des Mines Paristech (France): September 20 -23, 2017.

-Visiting Scholar to ENSAM -Paris (France): September 24 -28, 2019.

-Visiting Scholar to ENSAM -Paris (France): March 19 -23, 2020.

-Fulbright visiting associate professor to Stanford University (USA): January 15 -May 15, 2021.

-Visiting scholar to Ryerson University, FRAMES lab (Facility for Research on Aerospace Materials and Engineering Structures, Toronto, Canada): June 1 -December 15, 2021.

-Visiting scholar to ENSAM -Paris (France), starting November 1, 2021.

• Ongoing R&D projects:

-AM Valore: ESI funded project to use digital twins and data-driven simulation for building automated tape placement application in the industry 4.0 framework. Participation in the framework of a scientific collaboration between NDU and ENSAM-Paris.

-GDmodeling: Good year funded project aiming to model the die cross section shape required to built a given green tire cross section. Participation in the framework of a scientific collaboration between NDU and ENSAM-Paris.

-Magnetic bearing: SKF funded project aiming to enhance magnetic bearing control in real time to provide zero contact rotation. The project contains a data driven modeling and a model reduction component. Participation in the framework of a scientific collaboration between NDU and ENSAM-Paris.

• Successfully accomplished R&D projects as associate professor at NDU:

-ASALA: European project for simulation and hybrid-twins modeling of metal 3D printing. Delivered a real-time simulation app of the final deformation of the 3D printed part based on the process parameters based on real physics and validated experimentally for steel and Titanium-Aluminum alloy. Supervision of a PhD thesis. Project with GKN Sweden and Tekniker Spain. The project was successfully accomplished in 2021. Personal contribution in building a real-time simulation abacus for real-time simulation-based control of 3D printing. Contribution as a collaborator with the PIMM group -Electric engine simulation: Parametric simulation and optimization of the magnetic field in the electric engine to reduce the magnetic induced vibration. Project performed for ESI Group. The project was successfully accomplished in 2020. Publication has lead to being contacted by SKF in a novel ongoing project. Contribution in the framework of the ESI-ENSAM chari through the construction of the 3D parametric simulation abacus.

-Buckling of aeronautic composite materials reinforced frame structures, fully datadriven modeling performed using a novel technique I created personally, outperforming all other tested methods. Benchmark project by Dassault Aviation and successfully accomplished in 2020. Implementation of the Tabu search optimized regression algorithm.

-Ultrasound welding of PE/PET plastic tapes. Parametric thermal simulation and adhesion built allowing to perform real-time industrial control. Project performed for SEALESTER, France. The project was successfully accomplished in 2021. Participation in the construction of the 3D parametric simulation abacus.

-Real-time ATP: ESI-group contract. The project created a real-time simulation platform of the automated tape placement using optical ray orientation and absorption, as well as heat transfer system. The results are validated experimentally. This was a first ever state of the art project going from the ray tracing to the bounding estimation in automated tape placement, personal contribution in all parts of the project. The project was successfully accomplished in 2021.

-HJIMO: CNRS-L project for the simulation of the mechanical impact on human knee joint, validated experimentally. Application on vertical jumps and Dabke, the Lebanese folkloric dance. Personal contribution in creating a high fidelity simulation which lead to being invited to write a book chapter in the CNRS biomechanics edition. The project was successfully accomplished in 2020.

-SOFIA: simulation of steel 3D printing in collaboration with Michelin and ENSAM-Paris. Achieved path parametrization with real-time simulation of the deformation as a function of the path. Personal contribution in creating the model, the code and the post processing. The project was successfully accomplished in 2021. Cosupervision of the PhD thesis of Khalil el Rai and construction of the thermomechanical coupled model and simulation.

-SATP: simulation for the automated tape placement process, CENTRALE INNO-VATION research contract. Delivered a real-time simulation app for the thermal fields in during tape placement using coupled laser absorption optical and thermal simulations. The project was successfully accomplished in 2019.

-LAY2FORM: Simulation of the hybrid mylti-layers consolidation of hybrid structures composite/aluminum alloys. Project performed in collaboration with ESI Group. The project was successfully accomplished in 2019. Participation in the frameowrk of the collaboration between NDU and ESI-group by creating the Ericksen flow 3D simulator.

-Advanced simulation of composite material forming processes, CENTRALE IN-NOVATION project. In this project I used Ericksen anisotropic flow to simulate composite manufacturing during the squeeze flow as well as void dynamics, which reproduced an experimental effect unseen before and never again until now by any other simulation. Project perfomed in collaboration with Delaware university and successfully accomplished in 2018.

-ESI Chair affiliation between Ecole Centrale Nantes (ECN) and ESI group, focusing in the domain of efficient simulation of composites manufacturing processes. The project was successfully accomplished in 2018. Participation in the framework of a scientific collaboration between NDU and Centrale Nantes.

-FAURECIA Chair between ECN and FAURECIA group focusing in thermoplastic powder consolidation. Delivered a coupled thermo-mechanical simulation with phase change of the process. Data-driven modeling of the process and material parameters. The project was successfully accomplished in 2016. Participation in the framework of a scientific collaboration between NDU and Centrale Nantes.

-SIMUTOOL (European project) focusing in the simulation of electromagnetic problems in laminates; co-supervision of a PhD student. Delivered a high fidelity simulation of microwave absorption in composite microwave curing. Project performed in collaboration with ESI Group and successfully accomplished in 2017. Partici- • Contribution to supervision started after the beginning of the thesis and thus not officially integrated as member of the supervision committee

Khalil el Rai

• University: Arts et Métiers Paris

• Subject: Modeling, simulation and shape optimization in additive manufacturing.

• Year of thesis beginning: 2018

• Supervisor: Francisco Chinesta

• PhD thesis supervision contribution: 50%

Dominic Strobl

• Bundesanstalt für Materialforschung und -prüfung (BAM) -Federal Institute for Materials Research and Testing (BAM), Berlin (Germany). In collaboration with TU Berlin doctorate school

• Subject: Modeling and optimisation of metal 3D printing using model reduction techniques.

• Year of thesis beginning: 2021

• Supervisors: Annika Robens-Radermacher and J org F. Unger

• PhD thesis supervision contribution: 50%

Master thesis supervision

Charbel Gerges

• University: Université Saint Esprit de Kaslik

• Subject: Fluid and thermal simulation to assess the thermal comfort in a controlled environment.

• Year: 2014

• Supervisor: Chady Ghnatios

Georges El Morr

• University: Université Saint Esprit de Kaslik

• Subject: 3D simulation of the linear friction welding using the PGD

• Year: 2015

• Supervisor: Chady Ghnatios

Alexis Khalil

• University: Lebanese International University

• Subject: Proper Generalized Decomposition applications in fluid mechanics

• Year: 2016

• Supervisor: Chady Ghnatios

Alaa El Lakkis

• University: Lebanese International University

• Subject: A P1-P1 stabilized formulation using the PGD for fluid flow

• Year: 2017

• Supervisor: Chady Ghnatios

Instructed graduate level courses

• Finite Elements method: The concepts and fundamentals of the finite element method with applications to problems in solid mechanics, fluid mechanics, and heat transfer.

• Model reduction techniques: This course offers a unified presentation of the concepts and general principles of constructing a model order reduction solution, common to all branches of heat transfer, solid and fluid mechanics. The main applications shown are based on composite material forming processes. The course arises from the need to provide a general preparation in computational mechanics for students who will pursue further works and graduate studies in computational mechanics and composite materials. It will offer the link that helps you pass from the classical finite elements to the generalized laws of reduced order models.

• Mechanics of composite materials: Introduction to composite materials. Lamina and laminate mechanical properties. Micromechanics. Mechanical and hygrothermal behavior of laminae and laminates. Lamina and laminate strength theories

• Numerical optimization with machine learning: Engineering optimization is of upmost importance in design applications. The students are initiated to gradient decent algorithms, and nature inspired techniques. The course ends with an introduction to machine learning based on artificial neural networks and regression techniques.

• Stress analysis: This course covers plane stress and plane strain problems, inelastic axial deformation, inelastic torsion, introduction to membranes, composite beams, curved beams, inelastic bending, shear center, absolute maximum shear stress and strain, theories of failure, fully stressed beams, the secant formula, inelastic bulking, stresses in cylinders, energy methods, Castigliano's theorem, and introduction to plates and shells are covered.

• Advanced heat transfer: This course covers classical and numerical techniques applied to the solution of steady and transient conduction problems. Multidimensional steady state and transient problems are considered. Duhamel's theorem, complex variables solutions, error function and Laplace Transform solution to PDE are further considered to solve conduction problems.

• Continuum mechanics: This course offers a unified presentation of the concepts and general principles common to all branches of solid and fluid mechanics, designed to appeal to the intuition and understanding of advanced undergraduate students in engineering.

The course arise from the need to provide a general preparation in continuum mechanics for students who will pursue further work in specialized fields such as viscous fluids, elasticity, viscoelasticity, and plasticity. It will offer the link that helps you pass from the classical point mechanics to the generalized laws of mechanics.

Instructed undergraduate level courses

• Introduction to thermal fluid sciences: Fundamentals of thermal sciences including thermodynamics, fluid flow and heat transfer for non-ME students: properties of pure substances, energy & power, first law of thermodynamics, systems/ control volumes, fluid statics, Bernoulli equation, momentum equation, introduction to conduction, convection and radiation heat transfer.

• Science of materials: This course presents an introduction to materials science for engineers. Topics include atomic bonding, crystal structure, and crystal defects in their relationship to properties and behavior of materials (polymers, metals, and ceramics), phase equilibrium, and non-equilibrium phase transformations.

• Statics :This course treats only rigid-body mechanics and forms a suitable basis for the design and analysis of many types of structural, mechanical, or electrical devices encountered in engineering. As the course name suggests, this course deals with the equilibrium of bodies that are either at rest or move with constant velocity. Therefore, this Statics course provides the students with the principles that treats the statics of particles and rigid bodies, trusses, frames, machines; centroids, centers of gravity; and friction.

• Mechanics of materials: Equilibrium, force and moment resultants. Stress/strain concepts, generalized Hooke's law, classification of material behavior. Axial loads, torsion of rods and circular bars. Pressure vessels. Bending and shear stresses in beams. Analysis of Statically statically determinate and indeterminate structures. Stress transformation, combined loading.

• Advanced Mechanics of materials: Analysis of more complicated problems in stress and strain. Energy methods, torsion of non-circular members. Shear center concept. Curved beams, thick cylinders and rotating disks. Contact stresses.

• Dynamics: Description of force, position, velocity and acceleration in fixed and moving reference frames. Kinematics and kinetics of particles, of collections of particles and of rigid bodies. Energy and momentum concepts.

• Mechanical vibrations: Free and forced vibrations for undamped and viscously damped single-degree -of-freedom systems. Conservation of energy approach and Rayleigh's method. Vibration of multi-degree-of-freedom systems, eigenvalue problems and mode shapes. Modal analysis.

• Numerical analysis: Numerical integration, non-linear equations, approximations and interpolations of numerical functions, linear systems of equations (direct and iterative methods). Applications with matlab programming and implementations.

• Fluid mechanics: Fundamentals of fluid mechanics: fluid properties, fluid statics and kinematics, inviscid flow, potential flow, simple viscous incompressible flow. Flow analyses: control volume analysis, differential analysis, dimensional analysis and similitude. Flow equations: continuity, momentum and energy equations.

• Supervision of engineering final projects, in average 6 projects per year.

ACADEMIC SERVICES

Committee services

• Academic advisory for more than 70 students

• Member of the Department committee (2015-2021)

• Member of the department curriculum committee (2015-2021)

• Member of the ABET outcome assessment committee (2015-2020)

• Member of the faculty of engineering graduate committee (2016-2019)

• Member of the University Council (2017-2019)

• Member of the Faculty of engineering strategic plan ad-hoc committee (2017-2021)

• Member of multiple other ad-hoc committees.

Review services

• Reviewer for the International Journal for Numerical Methods in Engineering

• Reviewer for Computer methods in advanced mechanics and engineering

• Reviewer for composites part B

• Reviewer for the international journal of material forming

• Reviewer for COMPEL journal (The international journal for computation and mathematics in electrical and electronic engineering)

• Reviewer for Advanced modeling and simulation in engineering science AMSES journal

• Reviewer in multiple international conference

• Session chair in multiple international conference

INVITED TALKS

• "Modélisation et simulation de la dépose de bandes thermoplastiques". Groupement de Recherche : mise en oeuvre des composites et propriétés induites, Nantes (France), 2011.

• "Utilisation de la PGD pour la simulation des problèmes thermiques non-linéaires". Invited lecture at LTN, polytechnic university of Nantes, Nantes (France), 2013 (Plenary lecture).

• "Modélisation de la consolidation des matériaux composites thermoplastiques". Forum Plasturgie Composites (FPC), Paris (France), 2013.

• "Advanced Modeling and simulation of Materials and Processes". TCOP thematic conference, Piriac sur mer (France), 2013 (Plenary lecture).

• " Non-intrusive model order reduction technique for industrial applications". 5 th model reduction, POD, PGD conference (Mortech 2019), Paris (France), November 2019 (Plenary lecture).

• " Non-intrusive ePGD file construction through Morphing of space". ESI flagship conference PUCA 2019, Tokyo (Japan), November 2019 (Keynote Speaker).

• " Advanced interpolations: morphing vs. optimal transport". First international workshop on technologies enabling digital twins, Santo Dominigo (Dominican Republic), December 2019 (Keynote Speaker).

• " Real-time simulation of WAAM metal 3D printing". ESI scientific days 2019, Paris (France), December 2019 (Keynote Speaker).

• " Additive Manufacturing Hybrid-Twin". Second international conference on simulation for additive manufacturing SIM-AM 2019, An ECCOMAS thematic conference and IACM special interest conference, Pavia (Italy), September 2019 (Keynote Speaker).

• " A Proper Generalized Decomposition method for real-time parametric simulation and digital twin construction". Invited to the mechanical engineering additive manufacturing team meeting at Federal Institute for Materials Research and Testing (BAM), Berlin (Germany), June 2021.

EDITORIAL WORK

• Editor in chief for Frontiers in Materials -Switzerland special issue entitled: "Advanced Materials Modeling Combining Model Order Reduction and Data Science".

• Organizer of the minisymposium MS72: Model Order Reduction -Challenges in Engineering and Industrial Applications, at the 8th ECCOMAS 2022, Oslo -Norway.

• Reviewer for the International Journal for Numerical Methods in Engineering

• Reviewer for Computer methods in advanced mechanics and engineering

• Reviewer for composites part B

• Reviewer for the international journal of material forming

• Reviewer for COMPEL journal (The international journal for computation and mathematics in electrical and electronic engineering)

• Reviewer for Advanced modeling and simulation in engineering science AMSES journal

• Reviewer in multiple international conferences

• Session chair in multiple international conferences

CERTIFICATION

• Certified in Machine Learning by Stanford University, 2020.

• Specialization in Mathematics for machine Learning by Imperial College London, 2020 (3 courses certification including multivariate calculus and principal components analysis...).

• Certified in Data Sciences by Johns Hopkins University, 2020 (8 courses certifications including Practical machine learning and data scientists toolbox...).

MEMBERSHIPS

• Member of "The Order of Engineers and Architects" in Lebanon (OEA)

• Member of "European scientific association for material forming" (ESAFORM)

• Member of "association des centraliens de Nantes" (ACN)

• Former president of the lebanese students association of Nantes (ANADYL)

General Introduction

Model reduction techniques have been shaping the advancement of sciences since the early beginning of this century. Many industrial fields benefit from the abilities offered by such approach to achieve near real-time, efficient and high fidelity simulations. Moreover, with the recent democratization of artificial intelligence and data-driven techniques, many industries and academicians are moving towards the adoption of artificial intelligence and data-driven modeling to improve their knowledge and existing models of different industrial processes. My area of research is therefore mainly focused on model reduction and artificial intelligence/data-driven techniques, while creating technologies that facilitate their combination, which is currently known as the "digital twins" paradigm. My research activities are markedly divided into three main subjects:

• Developing model reduction methods in computational mechanics

• Modeling and simulation of different materials behavior as well as material forming and characterization processes

• Coupling Data-driven modeling and artificial intelligence to simulations with application to digital twins in mechanical engineering

In fact, currently data-driven modeling is increasingly popular on a worldwide level. Datadriven and digitization of systems and processes is a compulsory path to improve the knowledge and the current models. However, in engineering, experimental data is expensive to produce, many models are available and highly perform, while real-time simulation is in high demand in the industrial sector. Thus, I use data to correct and improve already existing models, which is currently known as a hybrid-twins or digital twins approach. This approach has increased the use of the term industry 4.0, where the aim is the development of technologies facilitating the digital twins training, modeling and applications in real-time, for industrial applications.

The scientific part of this work is divided therefore into three chapters, covering the main aspects of each of the three tracks of research mentioned above, while taking into consideration that these three tracks are not independent from each other but rather complementary as shown in this document.

In chapter 3, I review the main contributions to the development of model reduction techniques, mainly through the use of the PGD. Many developments are preformed including the use in stabilized mixed formulations for fluid problems [START_REF] Ghnatios | A stabilized mixed formulation using the proper generalized decomposition for fluid problems[END_REF][START_REF] Ghnatios | A stabilized pgd mixed formulation of the navier-stokes equation, application to a lid-driven cavity[END_REF], the use of the PGD and domain decomposition in general for the simulation of complex shapes [START_REF] Ghnatios | Computer methods in applied mechanics and engineering[END_REF][START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF][START_REF] Ghnatios | Shape optimization of complex domains using the proper generalized decomposition[END_REF][START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF] and the extension of the PGD to non-intrusive complex-shapes applications [START_REF] Ghnatios | Spuriousfree interpolations for non-intrusive pgd-based parametric solutions: Application to composites forming processes[END_REF][START_REF] Ghnatios | Towards parametric rtm processes: The interpolative mapping[END_REF]. Other ongoing works aim to create a PGD improvement through the use of machine learning [START_REF] Ghnatios | Artificial intelligence based space reduction of structural models[END_REF]. For instance, the aim of the method is to find the 3D solution of a problem using only one solution of a 2D partial differential equation. Preliminary results are impressive as shown in chapter 3. Several publications are selected and shown in appendix A, related to the main developments in the PGD framework performed within the last few years.

In chapter 4, many applications are shown using the PGD method. First of all, composite materials manufacturing processes consists one of the most extensive applications in my research. For instance, automated tape placement is one of the most promising out-ofautoclave manufacturing techniques and has taken a large part of my research leading to many publications including sensibility thermal analysis [START_REF] Perez | Sensitivity thermal analysis in the laser assisted tape-placement process[END_REF] as well as deformation and adhesion [18,37,[START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Poulahon | A numerical approach for the evaluation of residual stresses in the automated tape placement process[END_REF]. Other simulations related to the consolidation of tapes (in either automated tape placement settings or other processes) are studied. For instance the squeeze flow in composite laminates in deeply researched with many publications resulting from this work [START_REF] Ghnatios | 3d modeling of squeeze flow of unidirectionally thermoplastic composite inserts[END_REF][START_REF] Ghnatios | 3d modeling of squeeze flow of multiaxial laminates[END_REF][START_REF] Ghnatios | 3d modeling of the squeeze flow of unidirectional and fabric composite inserts[END_REF][START_REF] Ghnatios | 3d modeling of squeeze flows occuring in composite laminates[END_REF]. Also void closing dynamics during squeeze flow is studied and a continuous model is derived leveraging the PGD domain separation abilities [START_REF] Ghnatios | Modeling consolidation and void dynamics during automated tape placement process[END_REF][START_REF] Ghnatios | A non-local void dynamics modeling and simulation using the proper generalized decomposition[END_REF].

Other composite manufacturing processes are also studied including for example industrial blade manufacturing processes (project funded by General Electric) [6], the LCM process [START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF], ultrasound compaction, friction welding and consolidation [START_REF] Justo | Study of the ultrasonic compaction process of composite laminates -part ii: Advanced numerical simulation[END_REF][START_REF] Tannous | A new friction welding process for soft packaging applications[END_REF] and electromagnetic curing of composites [17,[START_REF] Tertrais | Simulation of microwave heating of a composite part in an oven cavity[END_REF][START_REF] Tertrais | Simulation of microwave heating of a composite part in an oven cavity[END_REF][START_REF] Tertrais | On the proper generalized decomposition applied to microwave processes incvolving multilayered components[END_REF][START_REF] Tertrais | On the proper generalized decomposition applied to microwave processes involving multilayered components[END_REF]., with high fidelity solutions to the Maxwell equations [START_REF] Ghnatios | On the high-resolution discretization of the maxwell equations in a composite tape and the heating effects induced by the dielectric losses[END_REF]. Recent works aim on using the ultra sound compaction for plastic tapes welding and bag forming (ongoing PhD thesis collaboration) and improving microwave curing of composite materials.

Other applications leverage the ability of the PGD to simulate multidimensional problems in a convenient framework and realistic calculation time [27]. For example, the use of the PGD to simulate SMC processes [START_REF] Perez | Advanced modeling and simulation of smc processes[END_REF], magnetic field in rotating synchronous machines [START_REF] Sancarlos | Fast computation of multi-parametric electromagnetic fields in synchronous machines by using pgd-based fully separated representations[END_REF], stir friction welding [40], simulating defects and holes in metallic structures [30,[START_REF] Ghnatios | A reduced-order model manifold technique for automated structural defects judging using the pgd with analytical validation[END_REF][START_REF] Ghnatios | A reduced order model based manifold technique for automated defect judging, application to structural holes[END_REF] and real-time control of processes [START_REF] Ghnatios | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: Application to pultrusion[END_REF][START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven of thermal processes[END_REF][START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven control of material forming processes[END_REF][START_REF] Masson | Numerical tools for the control of the unsteady heating of an airfoil[END_REF].

Moreover, advanced material characterization and homogenization is addressed using the PGD algorithm, such as XLPE modified concrete [START_REF] Himo | 2d modeling of the thermal conductivity of xlpe-modified concrete[END_REF][START_REF] Zehil | Electromagnetic field propagation in a composite laminate and induced thermal field: Application to microwave composites processing[END_REF], micro-structure homogenization including effective conductivity and effective hydraulic diameter [10], and soft biphasic materials such as hydrogels and cartilage replacements [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]. Indeed, cartilage simulation and human joints were also addressed in [START_REF] Ghnatios | Knee joint injury risk assessment bymeans of experimental measurements and proper generalized decomposition[END_REF].

Recently, I started working on the simulation of 3D printing of metals and the potential shape optimization, possible through the use of 3D printing. Promising and impressive results are published in [START_REF] Ghnatios | Reduced order modeling of selective laser melting: from calibration to parametric part distortion[END_REF][START_REF] Rai | Stress-based on intrusive shape optimisation for additive manufacturing parts[END_REF]. Further development in this field, to address the 3D printing of continuous fiber reinforced plastics, are currently underway in collaboration with colleagues from Ryerson university [START_REF] Ghnatios | Warping estimation of continuous fiberreinforced composites made by robotic 3d printing[END_REF]. Appendix B shows two selected publications from the different applications related to the 3D printing and to biphasic materials studies, as ongoing efforts are being developed to extend these applications and improve their accuracy.

The last part of my research, which focuses on the use of data-driven techniques and machine learning to model new processes, enhance simulations and incorporate variability into the simulation, is addresses in chapter 5. Depending on the application, different techniques were created and used to model data where existing models fail, using either Tabu search optimized regressions [START_REF] Ghnatios | An efficient tabu-search optimized regression for data-driven modeling[END_REF][START_REF] Hage | An optimized statistical model for predicting composite modulus of elasticity[END_REF][START_REF] Hage | Statistically validated and optimized tabu search estimation of cutting tool life in turning[END_REF][START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF], classical fitting methods [START_REF] Ibanez | On the data-driven modeling of reactive extrusion[END_REF], and a newly created separated representation, incremental fitting technique, able to operate on low data limit [START_REF] Reille | Incremental dynamic mode decomposition: A reduced-model learner operating at the low-data limit[END_REF]. Recent efforts aim to incorporate physical laws into fitting techniques. One of the used techniques is the GENERIC method, where a fitting of a system of differential equation is performed satisfying the thermodynamics constraints. An application on biphasic materials is shown in [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF].

Digital twins applications are also created using either a simulation enhancement when experimental ground truth can be defined or do exist [START_REF] Ghnatios | A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials[END_REF], or incorporating variability into the simulation when experimental variability is large, along with correcting the simulation to model with higher fidelity the experiment [START_REF] Ghnatios | A nonparametric probabilistic method to enhance pgd solutions with data-driven approach, application to the automated tape placement process[END_REF]. Appendix C comprehend some selected publications related to the modeling of data using data-driven techniques solely, to digital twins application for simulation correction and to digital twins applications for incorporating variability into the simulation.

In the rest of this document, chapter 3 will review the main efforts related to the development of the PGD method, chapter 4 shortly reviews the developed applications of the PGD and chapter 5 reviews the works related to the advanced data-driven techniques, digital twin technologies and their applications.

Chapter 3

Development of model reduction techniques 1 Chapter Introduction

Despite the recent progress in computational power and hardware tools, many physical problems remain intractable. When it comes to high dimensionality problems, the number of degrees of freedom increases drastically, even while using a coarse mesh in each dimension/space. This phenomenon is known as the combinatory explosion or the curse of dimensionality. Thus, the industrial applications' need of faster and reliable simulations, at affordable computational cost, has lead to the creation of new set of solvers known as the model order reduction techniques, currently known as MOR. In fact, MOR techniques tend to circumvent the curse of dimensionality by using a separated representation of the physical domain. For example, a 3D solution is computed as a sum of products of 1D vectors, or 2D×1D vectors [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF].

The development of model order reduction techniques is an ongoing process. Many families of these techniques are available, and can be classified into two big categories: (i) An "a posteriori "MOR, which has to solve the partial differential equation several times before building the reduced basis of the problem. (ii) An "a priori "model reduction technique, which can built the reduced basis or the reduced representation of the problem without having to solve any direct problem [START_REF] Ryckelynck | On thea priori model reduction: Overview and recent developments[END_REF]. The Proper Generalized Decomposition also known as PGD is an "a priori "model reduction technique, which builds a separated representation of the solution of the problem, without having to compute it first [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]. Interest in the PGD and the number of applications where it is used is in perpetual increase, and is the main subject of this chapter. Several developments are presented, which lead to broadening the applications of the PGD in industrial and digitization frameworks.

Developments of a stabilized PGD method for fluid problems

It is commonly known that mixed formulation in fluid simulation problems requires an extra stability condition to satisfy within the mesh in a finite element framework. This condition is known as the inf-sup condition of the Ladyzhenskaya-Babuska-Brezzi, shortly LBB condition. The PGD method is used for mixed formulations using quadratic polynomial interpolations (P2) in mesh elements for fluid flow, while using a linear interpolation (P1) in the mesh elements to represent and approximate the pressure field [START_REF] Ibanez | Simulating squeeze flows in multiaxial laminates: Towards fully 3d mixed formulations[END_REF]. In [START_REF] Ghnatios | A stabilized mixed formulation using the proper generalized decomposition for fluid problems[END_REF], a stabilized mixed formula-tion for fluid problems, based on the use of Galerkin Least Squares, is derived for solutions in separated form and applied to multiple problems using the PGD. A P1/P1 formulation is derived using linear shape functions for interpolating both the pressure and fluid velocity fields.

The results show a stable formulation even for parametric problems involving geometrical parameters and requiring a mapping into an apparent domain to efficiently solve the problem. Interested readers can review the details of the algorithm as well as the examples in A section 1.

Later on, the work was extended into a Navier-Stokes formulation, with an application to a lid-driven cavity in [START_REF] Ghnatios | A stabilized pgd mixed formulation of the navier-stokes equation, application to a lid-driven cavity[END_REF] using the variational multiscale method along with the PGD formulation. The work addresses the steady-state solution in a lid-driven cavity using a P1/P1 formulation along with a stabilized weak form of the problem. The weak form of the problem is modified by adding an extra stabilization term and thus becomes:

   a(w h , v h ) + b(P h , w h ) = l(w h , f ) b(Q h , v h ) - n el e=1 γ e ∇Q h , ∇P h Ωe = - n el e=1 γ e ∇Q h , f Ωe (3.1)
The work shows a promising solution addressing the LBB conditions for P1/P1 meshes using a separated representation framework like the PGD. The work will be extended into industrial applications, and also will be further developed to address transient phenomenons in fluid flow problems. These developments are to be performed while using the simple P1/P1 formulation along with a stabilization method based on the use of the variational multiscale method.

Separation of hardly separable domains

The separation of dimensions in the PGD is one of its most appealing applications. The possibility of obtaining a full 3D solution as a sequence of 1D problems using a separation of 3D into 1D/1D/1D domains, or the separation of 3D into an in-plane-out-of-plane decomposition 2D/1D alleviates the computation cost. Therefore, multiple applications are derived in this field [16,[START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Ghnatios | Towards parametric rtm process: The interpolative mapping[END_REF][START_REF] Ghnatios | A reduced order model based manifold technique for automated defect judging, application to structural holes[END_REF][START_REF] Perez | Sensitivity thermal analysis in the laser assisted tape-placement process[END_REF][START_REF] Tertrais | Simulation of microwave heating of a composite part in an oven cavity[END_REF]. In this framework, the geometrical separation of the domains is intuitively possible when the material and geometrical properties are easily separable. However, the problem can have multiple separation complications, including geometrical ones, material properties, multiphase flow... In the rest of this section we review the main developments performed in each of the considered topics and the achieved benefits.

Geometrical complexity

In this section, we address the geometrical separation complexity, which is studies in several publications. In [START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF], the geometrical complexity was addressed using two methods, (i) an extension of the domain to an easily separable one, (ii) A change of coordinates into a more suitable system where the geometry is easily separable.

The problem addressed in [START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF] is transformed using the extension of domains into an ill-posed Cauchy problem, solved parametrically using the PGD in a large hexahedral domain.

Once the solution available in the large domain, we then identify the solution in the original small domain embedded in the large one, using an optimization algorithm along with the diffuse approximation stabilization technique.

The second method uses a change of coordinates to a more suitable framework. The approach used in this work fueled further development of the geometrical transformation The microstructure interface defined by a sinusoidal form, separated by mapping the domain into a prismatic one method, deemed a morphing of space approach. Several publications were derived from the use of the morphing of space to solve the microstructure, geometrical and material properties separation as discussed in section 3.2.

Geometrical and material complexity

Problem of material properties separability are in general addressed by using a singular value decomposition of the tensor representing the material properties. This approach appears to be cumbersome and reduces the efficiency of the PGD algorithm by increasing much the number of matrix multiplications involved and the number of linear operators. In [START_REF] Ghnatios | Advanced separated spatial representations for hardly separable domains[END_REF], a novel method for addressing the separability of the domain both geometrically and from material properties point of view was addressed. The results are based on the use of physically based morphing of space to address the domain separability. For example, if we consider figure 3.1, the real domain (x, y) is mapped into an apparent domain (r, s) such as the interface is prismatic. More details are found in 2.

In [10], the authors used a morphing of space to derive the effect conductivity and apparent or effective diameter of a rough tube, mimicking the blood veins with debris and partial clotting. The results showed a effective diameter lower than the lowest diameter found in the The morphing of space was also used in an electric engine simulation to achieve a parametric simulation of the magnetic field, a currently very interesting simulation to allow the understanding and control of electrical machine's vibrations. Figure 3.2 reviews the work published in [START_REF] Sancarlos | Fast computation of multi-parametric electromagnetic fields in synchronous machines by using pgd-based fully separated representations[END_REF], where the simulation is detailed. The simulation includes 4 extra process parameters with a total number of degrees of freedom exceeding the Avogadro number 6.02214076 × 10 23 . The simulation was performed on a normal portable PC within few hours of calculation.

Addressing the non-intrusive approach for the PGD framework

PGD is known to be intrusive in the original approach published in [45]. However, to address industrial challenges, a non-intrusive approach of the PGD is mandatory. For that aim, a non-intrusive version of the PGD is derived and used in several applications. For example, in [START_REF] Ghnatios | Spuriousfree interpolations for non-intrusive pgd-based parametric solutions: Application to composites forming processes[END_REF], the proposed method uses a Sparse Subspace Learning (SSL) interpolation technique based on multilevel methods to address the solution between computation grid points. The SSL-PGD is derived in [START_REF] Ghnatios | Towards parametric rtm processes: The interpolative mapping[END_REF]. The algorithm build a response surface with respect to a set of parameters µ using:

U(x, t, µ) = i=N i=1 X i (x, t) • M i (µ) (3.2)
Now knowing the solution field as a function of (x, t) at two different points of the parametric domain U 0 (x, t, µ 0 1 ) and U 0 (x, t, µ 0 2 ) one may write a level 0 approximation as:

U 0 (x, t, µ) = U 0 (x, t, µ 0 1 ) • ψ 0 1 (µ) + U 0 (x, t, µ 0 2 ) • ψ 0 2 (µ) (3.3)
Once a new value of U is known, for instance U 1 (x, t, µ 1 1 ), an error can be computed using: Now a correction can be addressed to reach the level 1 approximation:

Ū1 1 (x, t, µ 1 1 ) = U 0 (x, t, µ 1 1 ) -U 1 (x, t, µ 1 1 ) (3.4)
U 1 (x, t, µ) = U 0 (x, t, µ) + Ū1 1 (x, t, µ 1 1 ) • ψ 1 1 (µ) (3.5)
The algorithm was coupled to a smart morphing of space to address the issue of localization and the impossibility of interpolation. Such problems appear when we have a strong dependence of the solution on the topology of a given parameter. For example, injection point location in RTM can't be interpolated for such interpolation results in artifact like the appearance of different injection locations weighted with a coefficient ∈ [0; 1]. For example, the reader may refer to figure 3.3 where an interpolation for different values of the permeability is performed. We reader can clearly see the artifact in the classical interpolation to the top right of the figure.

The algorithm was successfully used in several application including RTM permeability identification, SMC simulation, crack propagation and many others [START_REF] Ghnatios | Spuriousfree interpolations for non-intrusive pgd-based parametric solutions: Application to composites forming processes[END_REF].

The use of machine learning to enhance the PGD and reduce intrusivity

With the recent development of artificial intelligence and machine learning algorithms, some combinations of model reduction techniques with machine learning are addressed to account for error reduction and reliability improvements, either by error correction [4] or by improvement of modal selection using the quantities of interest [START_REF] Hartman | A deep learning framework for model reduction of dynamical systems[END_REF]. However, to the best knowledge of the authors, the model reduction techniques were only used as a black box, externally coupled to the machine learning algorithms used. The inclusion of artificial intelligence routines inside a model reduction techniques to enhance its performance has never been addressed before.

In this work, we use the machine learning algorithms to enhance the performance of the Proper Generalized Decomposition (PGD) model reduction technique. In fact, PGD uses a fixed point iterative algorithm to built a reduced basis "a priori", without any previous knowledge of the solution. The employed rank-one update greedy algorithm requires the solution of several lower-dimensions problems than the original one. For example, a 3D problem is reduced to several 2D×1D problems. The use of artificial intelligence will enhance the performance of the classical PGD algorithm into the so called "Artificial Intelligence enhanced PGD" or "AI-PGD", where for example only one 2D solution is required to fully solve a full scale 3D problem. The required 1D solutions will be obtained from a previously trained neural network.

The work starts with a description of the test problem, and several 3D solutions used to train the neural network. An artificial intelligence neural network is trained from the computed 3D solutions to allow the real-time identification of the out-of-plane functions on the fly. Later on, an enhanced PGD algorithm uses the computed out-of-plane functions at every in-plane location, to compute the in-plane 2D solution. The result is a full scale 3D solution with the calculation time of only one 2D problem.

Description of an example thermal problem

The steady state heat diffusion problem will be solved over a thin square plate with a width of a = 1m and a height of h = 0.1m. For instance the in-plane dimension x = (x, y) ∈ [0, a] × [0, a] and the out-of-plane one z ∈ [0, h]. The thin plate lateral surfaces will be insulated, while the bottom surface will be subjected to a constant temperature and the top surface will be subjected to a linearly variable temperature. The used boundary conditions and the variation of the temperature on the top surface occurs along the y-direction and is represented using Equation 3.6 and shown in Figure 3.4.

     ∇U = 0 at x = 0, x = a, y = 0, y = a U(x, z = h) = 350y + 50 U(x, z = 0) = U 0 (3.6)
The differential equation governing the heat diffusion in a thin plate is Poisson's equation shown in Equation 3.7,

∆U = - q K (3.7)
where U is the temperature along the plate, q is the internal heat generation and K is the thermal conductivity. The temperature on the bottom surface, U 0 ∈ [0; 400] A parametric solution data set is divided into multiple independent sets to either train the neural network or to test it and ensure no over fitting is occurring. Moreover, the set of values of each parameter is normally distributed along its range to avoid the presence of outliers in the resulting sets.//

The solution of the problem can be approximated in classical PGD as a sum of product of an in-plane functions, G i (x), times out-of-plane functions, H i (z),

U(x, z) = m i=1 G i (x)H i (z) (3.8)
In the new technique, the approximation will consider only one function R(x), multiplying a function of the thickness z changing as a function of the in-plane position x:

U(x, z) = R(x) • S(x, z) (3.9)
The shape of the out of plane function, or also will be referred to as the generic function S i (z), known at every location x i , differ from a node to another, depending on the location of the node and the problem statement. If the shape of the generic function S i (z) is known at each node of the in-plane domain x i , it will be enough to solve the 2D problem to find the in-plane function, R(x), to obtain an accurate solution of the 3D problem.

A similar method to the one utilized in the 2D-1D PGD method to obtain the solution of a 3D problem will be implemented in the following problem while keeping in mind that the 1D function is known from previous solutions [25,[START_REF] Tertrais | On the proper generalized decomposition applied to microwave processes incvolving multilayered components[END_REF].

Non-homogeneous Dirichlet boundary conditions on any surface of the plate will be added to the solution using previously defined terms to avoid deterioration of the solution on boundary surfaces while solving for the in-plane function, as classically performed within PGD framework [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]. Thus, the solution, U(x, z) is described using Equation 3.10:

U(x, z) = R(x, y)S i (z) + N j=1 X j (x, y)Z j (z) (3.10)
Where N is the number of predefined terms required to describe the non-homogeneous Dirichlet boundary conditions. We assume the weight function, U * , is written as [45]:

U * (x, z) = R * (x)S i (z) + R(x)S * (z) (3.11)
Since the out of plane functions along z are all known, thus the variation of S i (z) becomes equal to zero and the weight function reduces to,

U * (x, z) = R * (x)S i (z) (3.12)
Without loss of generality, the previously described heat problem will be utilized, to describe the procedure of obtaining the solution.

Multiplying Equation 3.7 by the weight function, U * (x, z), and integrating by parts to obtain the weak form depicted in Equation 3.13,

Ω×z ∇U * • ∇UdΩdz = Ω×z U * q K dΩdz (3.13)
where

Ω = x × y = [0, a] × [0, a] and z = [0, h].
Replacing Equation 3.10 and Equation 3.12 in Equation 3.13 and evaluating the 1D integrals at every single node in the in-plane domain results in the following 2D equation,

Ω α(x) ∂R * ∂x ∂R ∂x + ∂R * ∂y ∂R ∂y ) + β(x)R * R dΩ = q K Ω γ(x)R * dΩ - N j=1 Ω a j (x) ∂R * ∂x ∂X j ∂x + ∂R * ∂y ∂X j ∂y + b j (x)R * X j dΩ (3.14)
where α(x), β(x), γ(x), a j (x) and b j (x) are all vectors resulting from the evaluation of integrals involving the 1D functions, S i (z) and Z j (z), at every single in-plane node. The 1D terms, that are introduced to Equation 3.14 as linear fields distributed along the in-plane nodes, are obtained as follows,

•

α i = z (S i (z)) 2 dz • β i = z dS i (z) dz 2 dz • γ i = z S i (z)dz • a j,i = z S i (z)Z j (z)dz • b j,i = z dS i (z) dz dZ j (z) dz dz Equation 3.
14 is a 2 nd order degree equation with only the in-plane function, R(x), as unknown. The 3D solution of the problem, U(x, z), is thus assembled using the computed in-plane function, obtained from solving Equation 3.14, and the already known out of plane functions at every node.

Training using the out-of-plane functions

To solve the 2D equation for the in-plane function, Equation 3.14, evaluation of the 1D integral terms is required. To evaluate these terms, the generic functions, S i (z), are obtained from previously computed solutions. A shallow neural network is trained on these solutions to predict the required functions for the newly defined problem.

Parameters

Values q(W/m 3 ) 0, 100, 300, 500, 700, 1300, 1500, 1700, 2000 K(W/m.

Data sets

The training of the neural network is performed using 52500 generic function extracted from 500 parameter combinations of the 6D defined problem U(x, y, z, q, K, U 0 ). The solution of such heat problem can be obtained using many different methods. In this work, the 2D-1D proper generalized decomposition model order reduction method, PGD, is utilized to obtain fast and accurate solutions [49]. The different values of parameters used to obtain multiple solutions belong to the ranges specified in previous list and are shown in Table 3.1.

Although the set of solutions, generated using the set of parameters shown in Table 3.1, is divided into three sets (training, validation and testing set), a new problem leading to an independent set constituting of 13125 generic function, obtained using values of parameters shown in Table 3.2, is solved and used to ensure that no overfitting is occurring in the network.

Neural network architecture

The generic function extracted at every node from the computed solutions possess a parabolic shape and thus can be approximated accurately with a 2 nd degree order polynomial, as shown in Equation 3.15. The maximum error occurring between the actual function and the 2 nd degree approximation is limited to 2.7 % where the error is computed as shown in Equation 3.16.

S h i (z) = c 0 + c 1 z + c 2 z 2 (3.15) Error = S i (z) -S h i (z) S i (z) (3.16)
The trained regression vanilla type feedforward multilayer perceptron consists of an input layer, one hidden layer and an output layer.

The input layer nodes are the variables required to describe the problem and specify the location of the node where the generic function is to be predicted. Since the top surface temperature varies only along y as shown in Equation 3.6 and Figure 3.4, the shape of the generic functions does not change with the x-ordinate and thus only the y-ordinate of the in-plane node is required. Thus, the input neurons are selected to be:

• Location of in-plane node specified by the y-ordinate.

• Value of parameters: the internal heat generation (q), material conductivity (K) and the bottom surface temperature (U (x,z=0) ). The simplification of the in-plane dependency to y coordinate only simplifies the illustrative example, wihtout any loss of generality. The neural network is later on required to predict the generic function at each node. Knowing that S i (z) is approximated using 2 nd order degree polynomial, the output layer of the network will consist of the approximation coefficients, c 0 , c 1 and c 2 .

A shallow multilayer neural network with only a single dense hidden layer and a finite number of neurons is able to predict accurately the solution of any continuous mapping, non-linear problem given that an appropriate activation function is utilized in [START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF].

A 10-fold cross validation is performed to compare between different neural network architectures, activation functions, performance evaluation, training algorithms and data processing. The results obtained from below described network were satisfactory and thus ceased our search for a better one.

The hyperbolic tangent activation function is one of the most commonly used activation functions due to its ability of approximating a non-linear function, its symmetry around zero and fast convergence [START_REF] Kalman | Why tanh: choosing a sigmoidal function[END_REF]. Thus, the hyperbolic tangent or tanh is utilized as our activation function. Moreover, the number of neurons in the hidden layer is limited to 30 in order to:

• Avoid overfitting by keeping the number of trainable parameters much less than the number of training equations.

• Obtain good accuracy while keeping the training time as short as possible.

Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], a variation of the standard backpropagation algorithm [START_REF] Demuth | Neural Network Design[END_REF], is used to train the network since it is 10 to 100 times faster than the usual gradient descent backpropagation method [START_REF] Hagan | Training feedforward networks with the marquardt algorithm[END_REF]. Normalization of the data features is also performed. The performance is evaluated using the mean of the squared error, MSE. A summary of the neural network utilized to predict the generic functions is shown in Figure 3

Evaluating the neural network performance

The trained neural network is tested on the independent data set obtained using the values of parameters shown in Table 3.2. The performance of the network is evaluated using the mean of the mean of the residual, over all nodes and samples, between the target generic function and the predicted generic function. The residual for generic function i at out of plane node j is computed as shown in Equation 3.17.

R i,j = S h i,j -S p i,j (3.17)
S h i is the generic function approximated using the target coefficients and S p i is the generic function approximated using the predicted coefficients. Table 3.3 summarizes the obtained mean along with the maximum obtained residual between all predicted functions. The normalized values are obtained by dividing the values by a reference temperature chosen to be the highest temperature at the top surface of the plate, i.e. 400 • C. Figure 3.6 compares the predicted generic function with maximum residual with the target generic function.

AI-PGD results

After obtaining the 1D generic functions using the neural network described in subsection 5.2, the 2D finite element equation, Equation 3.14, can be easily solved for the in-plane function, R(x). Solutions for all problems governed by values of parameters stated in Table 3 5.84% Table 3.4: Summary of obtained maximum error on all the nodes of the in-plane domain assembled and each obtained solution is compared with the reference PGD solution.

To compare both solutions, the error at each node is computed using Equation 3.18. The mean square error, for every problem, is computed and the maximum between these values is shown in Table 3.4. Also, the infinity norm error, for every node in the problem, is computed and the maximum between these values is shown in the same table.

Error = U pgd -U generic U pgd (3.18)
The resulting error, with a maximum mean error of 0.36%, is very satisfactory compared to the decrease in computation time and complexity. Solving the problem using PGD would require the computation of multiple enrichment terms while iterating between the in-plane function and the out of plane function for each term using the fixed point iterative algorithm [44]. In the following method, only one term is computed and no iteration is required. Figure 3.7 shows the solution with the maximum mean error. Showing the reference PGD solution is meaningless due to the similarity between both solutions. A plot of the maximum error at each in-plane node for the above solution is shown in Figure 3.8 where a maximum of 2.99% occurs at the lower edge of the plate.

Discussion

The AI-PGD solution allows the reduction of a full 3D simulation into a single 2D simulation. However, this comes with the price of previously training a neural network to predict the outof-plane shape function at every single node of the 2D domain. One may note that if the functions S i (z) are exact, the correct 3D solution can be found by only multiplying the S i (z) by a value equal to one, and thus R(x) = 1 would yield the exact 3D solution of the problem. However, since S i (z) is not exact, a solved R(x) using a 2D in-plane problem enhances the final 3D solution by reducing the maximum error. This hypothesis is tested and we found that, in case of a very well-trained neural network, assembling the generic functions, S i (z), at each node, might be enough to obtain a solution for the problem and hardly noticeable enhancement occur by solving for the in-plane function, R(x). However, in case the prediction of the generic functions is not accurate, solving for the in-plane function, using Equation 3.14 will reduce the error and enhance the solution.

In fact, the generic functions for the same set of problems are obtained using a less trained network consisting of 10 hidden neurons. Solving for the in-plane function, R(x), succeeded in reducing the mean error for 96.8% of the tested problems as compared to using R(x) = 1. Comparison of the accuracy of both networks is shown in Table 3.5 along with the rate of success in reducing the mean error along a whole set of 125 problems.

Chapter conclusion

This chapter illustrates the developments achieved in the PGD framework, leading to the possibility of simulating a mixed formulation of fluid flows using a stabilized P1/P1 formulation. This formulation shall be extended to higher Reynolds number applications and therefore to industrial and aeronautics problems.

On the other hand, this chapter illustrates several improvements on the PGD leading to the possibility of constructing parametric solutions when interpolation is apparently impossible and even geometry is inseparable at first site. This leads to the construction of offline computation abacus, even for possibly non separable and non interpolable solutions. The development is also extended into the SSL-PGD non-intrusive framework, where the applications are numerous. An enhancement of the non-intrusive version of the PGD using autoencoders and machine learning algorithms in currently an ongoing work.

Finally a current and recent development aims to find the PGD solution of a full 3D problem using only one 2D solution. This is achieved through learning the out-of-plane 1D functions by using a neural network, previously trained to learn these 1D functions. The method is very promising and will be developed further in different applications.

Chapter 4

Advanced simulation of materials and processes 1 Chapter introduction

In this chapter, I review few applications revolving around the second axis of my research. This axis focuses on leveraging model order reduction and advanced simulation techniques, to simulate complex materials as well as advanced material processing methods. The selected applications are of current major interest to the scientific community. For example, biphasic materials are considered solely as the most suitable components for use as a replacement of several human body tissues, like contact lenses and cartilage replacements for example. On the other hand, 3D printing of metallic components is being actively researched to cope with the growing need of fast manufacturing of complex parts.

The two sections of this chapter are based on the modeling of complex materials involving more than one phase, coupled mechanism and the need to have an extremely fine mesh in the thickness direction to capture accurately the simulated material behavior. Other applications revolve around additive manufacturing, automated tape placement simulations, electromagnetic propagation in composite materials, optimization of stresses and assessment of stresses in plates around holes, human joints modeling... In the rest of the chapter, we review briefly some of the accomplishments in this axis of research and discuss future development plans.

Modeling and simulation of biphasic materials

In this section we review the work performed on the modeling and simulation of biphasic materials, mainly hydrogel brushes, believed to be one of the most promising human cartilage replacements. Hydrogels are also currently used in contact lenses among other applications in the human body. Biphasic materials, as their name indicates, are made of two phases, solid and fluid. It is believed that the fluid pressurization is the main load carrying phenomenon in hydrogel brushes [START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF][START_REF] Li | Poly(acrylamide) films at the solvent-induced glass transition: adhesion, tribology, and the influence of crosslinking[END_REF]. However, the exact estimation of the fluid pressurization inside a porous medium requires a good estimation of the permeability in that medium. In [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF], a micro indentation test (indentation on a micro-scale) is used to estimate the hydrogel reaction force as a function of the indentation depth, for different indentation velocities ranging from 0.5µm/s to 80µm/s.

• indentation depth H n • indentation speed U n • permeability K n • Domain radius R n
Therefore, the result is a generic abacus containing the pressure fields P as a function of the process parameters written in a separated form: Once the solutions available, one should mention that the separated form facilitated the computation of the gradients, which are now readily available. In fact, the derivatives of P with respect to K n can be written as:

P = i=N i=1 X i (x, y) • Y i (y) • H i (H n ) • R i (R n ) • U i (U n ) • K i (K n ) (4.1)
∂P ∂K n = i=N i=1 X i (x, y) • Y i (y) • H i (H n ) • R i (R n ) • U i (U n ) • ∂K i (K n ) ∂K n (4.2)
The gradients are therefore coupled with a gradient decent algorithm to identify the permeability in the porous medium. The identification in instantaneous using a Newton algorithm. More details can be found in B of this manuscript.

3D printing of metals

In this section we review the accomplishments performed in 3D printing of metals, where a thermal simulation, a mechanical deformation and a parametric path simulations are performed.

Thermal simulation

A thermal simulation is performed a priory to achieve good understanding of the process and the coupling involved in the mechanism. The thermal simulation uses a non linear variation of the material properties along with a phase change heat generation using a classical model by [21,[START_REF] Pinheiro | Mould heat transfer with mould flux lubrication[END_REF][START_REF] Radovic | Numerical simulation of steel ingot solidification process[END_REF][START_REF] Riede | Material characterization of aisi 316l flexure pivot bearings fabricated by additive manufacturing[END_REF]. The governing equation is therefore given by:

ρC P ∂T ∂t -∇ (k • ∇T ) = Q -(k • ∇T ) • n = Q laser T base = 293K -(k • ∇T ) • n = h(T -T ∞ ) T (t = 0) = 293K; T ∞ = 293K (4.3)
The result of the thermal simulation shows that a "steady state" is achieved along a path few millimeters after beginning a path deposition, figure 4.3 shows a single path deposition with element activation (the elements beyond the deposition point have zero conductivity). An abacus corresponding to the maximum reached temperature as a function of the process parameters is also illustrated in figure 4.4.

Mechanical deformation, a meso-scale

Once the thermal simulations available, the study of the mechanical deformation on a mesoscale is possible. The used model accumulates the deformations through the use of phase changing materials during solidification [START_REF] Ghnatios | Reduced order modeling of selective laser melting: from calibration to parametric part distortion[END_REF][START_REF] Rai | Stress-based on intrusive shape optimisation for additive manufacturing parts[END_REF]. The model uses an incremental building approach to simulate the deformation of the building block until reaching a shape allowing to extract a relevant representative volume element. Such a construction is illustrated in figure 4.5, with a zoom showing the quality of the mesh in 4.6. This refined mech wouldn't be possible to use when the solution is computed on a normal portable PC, within realistic computation time, without the use of the PGD and the in-plane-out-of-plane decomposition.

More details can be found in section 2.3 of this manuscript.

Mechanical deformation, a macro-scale

To illustrate the deformation on a macro-scale, the simulation uses the stress tensor obtained from the meso-scale simulation and applies it to the large part, where a desired final shape is to be obtained [START_REF] Ghnatios | Reduced order modeling of selective laser melting: from calibration to parametric part distortion[END_REF]. The results are comparable to the experimental ones as shown in section 2.3. The model is used to study the deformation on a large cylinder, on a wall and on a cylinder with ring, as illustrated in figures 4.7, 4.8, 4.9, 4.10.

Composite materials manufacturing and simulation

The composite materials manufacturing methods fill a large part in my research portfolio. Many processes are simulated, using different methods and techniques. To cite a few, I can mention for example an advanced simulation of the squeeze flow in Newtonian and non Newtonian regime [START_REF] Ghnatios | 3d modeling of squeeze flows occuring in composite laminates[END_REF], anisotropic fluids subjected to squeezing [START_REF] Ghnatios | 3d modeling of squeeze flow of multiaxial laminates[END_REF][START_REF] Ghnatios | 3d modeling of the squeeze flow of unidirectional and fabric composite inserts[END_REF], void closing dynamics in automated tape placement [START_REF] Ghnatios | Modeling consolidation and void dynamics during automated tape placement process[END_REF][START_REF] Ghnatios | A non-local void dynamics modeling and simulation using the proper generalized decomposition[END_REF], automated tape placement thermal and mechanical deformations [37,[START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Perez | Sensitivity thermal analysis in the laser assisted tape-placement process[END_REF][START_REF] Poulahon | A numerical approach for the evaluation of residual stresses in the automated tape placement process[END_REF], electromagnetic simulation of composite manufacturing [17,[START_REF] Ghnatios | Towards parametric rtm processes: The interpolative mapping[END_REF][START_REF] Tertrais | On the proper generalized decomposition applied to microwave processes incvolving multilayered components[END_REF][START_REF] Tertrais | On the proper generalized decomposition applied to microwave processes involving multilayered components[END_REF], RTM [START_REF] Ghnatios | Towards parametric rtm processes: The interpolative mapping[END_REF], modified concrete [START_REF] Himo | 2d modeling of the thermal conductivity of xlpe-modified concrete[END_REF][START_REF] Zehil | Electromagnetic field propagation in a composite laminate and induced thermal field: Application to microwave composites processing[END_REF], among others. In what follows is a brief review of some selected examples of my contributions in the composite manufacturing field.

Automated tape placement

The automated tape placement has been a center of interest for many researchers as it is a promising composite manufacturing technique, which has the potential to create a oneshot composite material, without the need of an autoclave. Since my PhD thesis [START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF], many developments are preformed based on the simulation of the composite manufacturing using the automated tape placement. Recently, a an extension of the method is performed for fiber filament winding. In winding, the equations are written in cylindrical coordinates, to ease the modeling of the process, and the use of the in-plane-out-of-plane decomposition. Therefore, the weak form of the problem is written as:

U(r, θ, z, T in , T b , w, Q 1 , Q 2 ) = i=N i=1 R i (r)Z i (θ, z)T i (T in )U i (T b )Q i (Q 1 )P i (Q 2 ) (4.5)
The boundary conditions of the problem are shown in figure 4.11. The solution is obtained within few minutes, less than an hour, using a normal portable PC. This solution contains the equivalent of 32 × 10 12 3D finite element solutions, each having 36100 degrees of freedom.

The solution for different positions are sections is illustrated in figures 4.12.

squeeze flow of multiaxial stands and fabrics

We assume that the the fiber directions is noted by p and q, and their relative rotation creates a distributed torque all along the continuous fibers, that depends on the the relative rotary velocity through a friction coefficient µ. The fibers in direction p are assumed to have a length 2L p whereas those parallel to q of length 2L q , as depicted in Fig. 4.13.

The resulting torque can be represented by a couple of forces applying at the fiber ends F f (pL p ) and F f (-pL P ) = -F f (pL p ), where the super-index • f indicates its frictional nature. These forces act in the perpendicular direction to p, designed by p ⊥ , with p × p ⊥ = (0, 0, 1) T . Thus, being M the z-component of the torque vector M, frictional forces write for both family of fibers:

F f (pL p ) = M 2Lp p ⊥ F f (qL q ) = -M 2Lp q ⊥ , (4.6)
that added to the hydrodynamic forces gives:

F h (pL p ) = ξ(v 0 + ∇v • pL p -v G -ṗL p ) F h (qL q ) = ξ(v 0 + ∇v • qL q -v G -qL q ) , ( 4.7)
results

F(pL

p ) = ξ(v 0 + ∇v • pL p -v G -ṗL p ) + M 2Lp p ⊥ F(qL q ) = ξ(v 0 + ∇v • qL q -v G -qL q ) -M 2Lp q ⊥ , ( 4.8)
being F(-pL p ) = -F(pL p ) and F(-qL q ) = -F(qL q ). The linear momentum balance equation applied to each rod implies again v 0 = v G . The angular momentum balance implies again:

F(pL p ) = ξ(v 0 + ∇v • pL p -v G -ṗL p ) + M 2Lp p ⊥ = λ p p F(qL q ) = ξ(v 0 + ∇v • qL q -v G -qL q ) -M 2Lp q ⊥ = λ q q , ( 4.9)
that multiplying the first by p and the second one by q result in:

λ p = ξL p ∇v : (p ⊗ p) λ q = ξL p ∇v : (q ⊗ q) , ( 4.10)
that implies the rotary velocities are written by:

   ṗ = ∇v • p -(∇v : (p ⊗ p))p + M 2ξL 2 p p ⊥ q = ∇v • q -(∇v : (q ⊗ q))q -M 2ξL 2 q q ⊥ . (4.11)
The frictional couple M is assumed proportional to the relative rotary velocity:

M = -µ( ṗ -q ). (4.12)
Equations (4.11) and (4.12) define the nonlinear kinematic model to be solved at each location x ∈ Ω, from which directions p and q defining the fabric microstructure can be updated at time t + ∆t from the known values at time t.

The problem was solved for large composite parts, patches in the domain either unidirectional strands or woven patches. Figure 4.14 illustrates the solution for a given woven patch compressed inside a non Newtonian fluid domain.

Electromagnetic heating high fidelity simulations

Electromagnetic heating is currently being investigated to study the possibility of ensuring through-volume uniform heating of the composite preform through microwaves. The simulations focused on the high fidelity representation of the heterogeneous composite materials. Therefore, a mesh of the domain showing different material phases is performed using the domain decomposition. Such as a 3D domain is decomposed into 1D × 1D × 1D, allowing therefore the use of very fine mesh while using limited computational resources, and obtaining the results within a reasonable time frame [START_REF] Ghnatios | Advances in Computers and Information in Engineering Research, ASME ACIER book series[END_REF].

Electromagnetic heating, effect of fiber orientation with respect to the wave length

In this part we review the heating modeling of composite materials when using an electromagnetic heat source. The current work focuses on the use of microwave heating, while considering the potential influence of changing the electromagnetic frequency ω.

Electric simulation and process parameters

In this part we simulate the electric field propagation at the fiber scale level using two approaches: (i) an extruded 2D section and (ii) a full 3D scanned part obtained by using micro tomography.

Simulation of a composite part with extruded in-plane section

In this section, we consider the material cross section shown in figure 4.20. The depicted cross section has a height of L y = 0.2mm and a width of L x = 1.25mm. The simulated part is an extrusion of this 2D cross section with a depth of L z = 50mm.

The material properties in-plane are separated using a singular value decomposition. Therefore, for any material property p we can write:

p = X p i (x)α i Y p i (y) (4.13)
with X p i the normalized x components of the singular value decomposition, Y p i the normalized y components and α i the singular values. The truncation is performed for α i /α 1 < 10 -6 .

To identify the fiber orientation effect on the electric field propagation in the part, an extremely fine mesh consisting of 7000 nodes along x direction, 4000 nodes along y direction and 4000 nodes along the z direction is used, representing the equivalent of 342 × 10 9 degrees

                        E y (x = 0, x = L x ) = 0 E z (x = 0, x = L x ) = 0 E x (y = L y ) = cos 2πL x x Lx E z (y = L y ) = 0 E x (z = 0, z = L z ) = 0 E y (z = 0, z = L z ) = 0 (4.14)
Γ * is the part of Γ where no Dirichlet boundary conditions are prescribed. On the other hand, set two writes:

                           E y (x = 0, x = L x ) = 0 E z (x = 0, x = L x ) = 0 E x (y = L y ) = 0 E z (y = L y ) = cos 2πL z z Lz E x (z = 0, z = L z ) = 0 E y (z = 0, z = L z ) = 0 (∇ × E) • n × Ē * | Γ * = 0 (4.15)
In equations (4.14) and (4.15), L x and L z are respectively the length scale along x and z computed using the wave frequency ω and its wave length λ = c ω , using:

L x = Lx λ L z = Lz λ (4.16)
being c the light speed. We therefore simulate the application of the real wave propagating along a unidirectional axis at the top surface of the part y = L y . To identify the dependence on the frequency, two frequencies are selected for illustrative purposes: ω 1 = 2.45GHz and ω 2 = 24.5GHz. The material properties are given in table 4.1.

Parameter Value // to the fibers Value ⊥ to the fibers Value in the matrix

ǫ r 20 + i 20 + i 1 µ r 1 1 1 σ 10 3 10 2 10 -4
Table 4.1: Properties of the fiber and matrix phases in the parallel and perpendicular direction to the fibers axis With the values indicated in table 4.1, the set 1 of boundary conditions, considering a propagation along x direction and thus perpendicular to the fibers are illustrated in figures 4.21 and 4.23 for ω 1 and ω 2 respectively, while figures 4.22 and 4.24 illustrate the electric field propagation along the thickness of the domain for ω 1 and ω 2 respectively.

The same results for the second set of boundary conditions are given in figures 4.25 and 4.27 for ω 1 and ω 2 respectively, while figures 4.26 and 4.28 illustrate the electric field propagation along the thickness of the domain for ω 1 and ω 2 respectively.

From the results given in figures 4.21 to 4.28, one can clearly state two main conclusions:

• The electric field propagation inside the part is deeper for lower frequencies, and the electric field exhibits therefore a higher intensity for lower frequencies

• The electric field propagation inside the part is higher when the wave propagates parallel to the fiber orientation

The first conclusion is quite intuitive as the higher the frequency the lower becomes the penetration depth, while the second conclusion seems more related to the possibility of the wave to propagate between parallel fibers when oriented in the same direction as these fibers.

Simulation of a real 3D scanned composite part

In this section we consider a 3D micro tomography scan of a composite part and we identify the different phases using k-means algorithm. After selecting 3 phases, we identify the phases as shown in figure 4.29. The k-means algorithm identified three phases illustrated in red, dark blue and light blue. The red is the carbon fibers phase, the dark blue is the matrix phase and the light blue represents the air bubble entrapment. The air is considered to have ǫ air r = 1, µ air r = 1 and σ air r = 0. With these assigned properties we can use an in-plane-out-of-plane decomposition using the PGD method.

First of all, the material properties in the 3D domain are separated using singular value decomposition. Again, for any property p, one may write:

p(x, y, z) = np i=1 U T i α i V i (4.17)
with U i (x, y) being the normalized singular vectors decomposition in the in-plane domain (x, y), and V i (z) the normalized singular vectors of the out-of-plane domain z. The values α i are the singular values or the weights of the i-vectors. We truncate the material properties representation at i for α i /α 1 < 10 -6 . When having the material properties in the in-planeout-of-plane separated form, one can easily apply the proper generalized decomposition for the considered domain [24,[START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF].

The considered specimen has the dimensions L x = 1.6mm, L y = 0.512mm and L z = 2.1mm. The results are given in figures 4.30 and 4.31 for the same boundary conditions given in equations (4.14) and (4.15) in section 4.4, respectively.

From figures 4.30 and 4.31, one can reinforce the conclusion from section 4.4: the propagation inside the plate is higher when imposing a wave propagating in the direction parallel to the fibers, even if the fibers are not fully aligned. Another conclusion can be inferred: when the fibers are not aligned, the in-depth propagation is smaller than the one seen in the case of fully aligned fibers.

Thermal simulation

The electric field is leveraged to compute the volumetric heat source Q as described in equation (??). Without loss of generality, the simulation in this section is performed on L z = 50mm, the extruded fibers domain.

The propagation of the electric field operates differently in the domain, and consequently the heat generation will also depends strongly on it and on the material parameters. For example, if we consider the second term of equation (??): From the heat generation, the temperature field is computed using the heat equation and the boundary conditions given by equations 4.34 and 4.35 at t = 50s. The time discretization consists of 500 nodes within the considered 50 seconds time-frame. The heat generation terms Q are multiplied by 380 × 10 3 V in both cases, to consider a more realistic heat source.

1 2 E T • σ • Ē = 1 2    E x E y E z    T •    σ ⊥ 0 0 0 σ ⊥ 0 0 0 σ / /    •    Ēx Ēy Ēz    = 1 2 E x .σ ⊥ . Ēx + E y .σ ⊥ . Ēy + E z .σ / / . Ēz

Conclusion

In this work, we used the PGD to simulate microwave heating of composite parts with a fiberlength resolution. The simulation clearly proves the ability of the PGD to solve high fidelity problems with the equivalent of more than 10 10 degrees of freedom, within about an hour, on a standard laptop. The results also proves the effect of the frequency on the wave penetration, as well as the effect of the relative orientation between the external imposed wave and the fibres. Waves parallel to the fiber orientation can penetrate deeper into the material, while waves normal to the fibres are absorbed at the surface neighborhood.

Chapter Conclusion

In this chapter, we reviewed few selected applications of material and process modeling. Several other applications are also available including the simulation of the squeeze flow of composite materials and composite inserts using anisotropic fluids, many other simulation of microwave curing of composite materials [START_REF] Ghnatios | Electrical fields simulation in heterogenous domains using the proper generalized decomposition[END_REF], and others.

The future plan in this axis would be first increasing the number of possible applications in composite materials, biphasic materials, non destructive testing and characterization... Current works are focusing of the simulation of automated placement of fibers in the framework of 3D printing of continuous fiber reinforced plastics. Other applications are focusing of RTM optimization and permeability detection. Moreover, future plans will also include a review of the simulated applications to enhance the simulation using digital twins paradigm, enriching the simulation with data-driven modeling, using techniques detailed in chapter 5.

Chapter 5 Digital twins and data-driven techniques 1 Chapter introduction

Data sciences and artificial intelligence are nowadays the newly booming combination in research. With the impressive progress in computing power and numerical algorithms, handling big data through regression and machine learning is getting popular. With the ease of use and high accuracy, machine learning-based models are competing with well established constitutive models. Some recent works tend to replace long known constitutive models with data-based identification [START_REF] Ibanez | A manifold learning approach to data driven computational elasticity and inelasticity[END_REF][START_REF] Ibanez | Data-driven non-linear elasticity: constitutive manifold construction and problem discretization[END_REF], while others correct the uncertainty or ignorance in the models using regression on error values [START_REF] Moore | A hybrid point-particle force model that combines physical and data-driven approaches[END_REF]. However, some other materials constitutive models are highly non linear and not well established. In such case, the data-driven approach constitutes an appealing route. In this chapter, we explore several possibilities of creating digital twins and data-driven applications:

• A data-driven modeling, based on the use of fitting and regression techniques, when modeling is poor, hard to accomplish or non existing. A model method for creating optimized regressions is illustrated in section 1. An application of this technique on the modeling of the drilling of basalt rocks is illustrated in [START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF]. Another method based on incremental separated representation is discussed in [START_REF] Reille | Incremental dynamic mode decomposition: A reduced-model learner operating at the low-data limit[END_REF], which can operates on low data limits. An application is shown in [START_REF] Ibanez | On the data-driven modeling of reactive extrusion[END_REF].

• A hybrid modeling or digital twin modeling based on the correction of existing models using error approximations from data [START_REF] Ghnatios | A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials[END_REF]. Section 2 reviews the method for a biphasic hydrogel behavior modeled using nanoindentation. The used method is suitable when models are good enough to have a small margin of error with respect to the original model size.

• A digital twin modeling using non parametric probabilistic approach for separated solutions applications [START_REF] Ghnatios | A nonparametric probabilistic method to enhance pgd solutions with data-driven approach, application to the automated tape placement process[END_REF]. This method is suitable when the experimental data includes high variability and therefore a correction is not suitable. Instead, incorporating the variability in the simulation results, while correcting the simulation if required using an "averaging" method is suitable. Section 3 illustrates the method in the framework of an continuous unidirectional fiber automated tape placement process.

In what follows, we review shortly some of the selected methods for the creation of datadriven fitting and classification techniques, as well as some hybrid twins techniques used for incorporating variability in the solution, or for enhancing the representation of a solution.

Creating data-driven techniques 2.1 A Tabu-search optimized regression

Regressions are hard to optimize, as selecting novel terms and enhancing the fitting may lead to over fitting, while reducing the number of terms in a representation may decrease the ability to correctly represent a data set. The classical trade-off problem between eh high variance representation and the high bias one is therefore to be tackled through an enhanced and automated optimized representation through the use of a Tabu-search meta-heuristic algorithm [START_REF] Ghnatios | An efficient tabu-search optimized regression for data-driven modeling[END_REF]. This algorithm is coupled with the use of statistical relevance criterion to omit or accept a newly added term to the regression.

A Tabu-search operates on the possible data representation terms. Hence the algorithm starts by creating a pool σ, which consists the potential candidates to the fitting representation σ i . A potential representation of σ can be a set of functions (polynomials, cosine, sine...) of the fitting variables.

Later on, one shall define a neighborhood of variation of the existing regression F. The neighborhood can be for example adding one or more extra terms σ i to F, removing one or more terms from F and substituting terms. Every iteration consists of a neighborhood testing and assigning a cost J(F). The best cost is recorded as S best and would correspond to the best selected solution. At every iteration, the terms are tested for their pvalue and assess their statistical relevance to the solution. The algorithm is illustrates in figure 5.1.

A possible cost function that is used in our work is selected as [START_REF] Ghnatios | An efficient tabu-search optimized regression for data-driven modeling[END_REF]:

J = w * R 2 (F) + (1 -w) * MAP E (5.1)
where R 2 is the regression's coefficient of estimation and MAP E is the mean absolute percentage error. w ∈ [0; 1] is a weight value, selected equal to 0.3 in previous works [START_REF] Ghnatios | An efficient tabu-search optimized regression for data-driven modeling[END_REF].

The selected regression algorithm has shown its superior performance against published regressions in the literature [START_REF] Hage | Statistically validated and optimized tabu search estimation of cutting tool life in turning[END_REF], but also its high performance in real-life problems [START_REF] Hage | An optimized statistical model for predicting composite modulus of elasticity[END_REF][START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF].

An incremental dynamic mode decomposition

In general, many processes operates on the low data-limit. In fact, experiences and measurements in mechanical engineering are expensive and must be extremely controlled to avoid environmental factors. Therefore, these considerations renders many measurements as either very expensive or not feasible on the required scale. Therefore, the need of developing an approximation tool as good as possible, able to generate accurate approximations for the low available amount of data has surfaced.

A method operating incrementally by correcting the previous error at every appearance of a new set of data is therefore created. The method operates incrementally such as it starts from the basic physical principles:

KU = F => K(UU T ) = FU T (5.2)
Later on, one can create a lower rank approximation, K LR , of the linear relations between the inputs F and the outputs U for the n known set of measurements, such as:

K LR = n i=1 C i × R i (5.3)
For a one set of known couples (U 1 , F 1 ) one may write:

C 1 R T 1 U 1 = F 1 (5.4)
The problem shown in equation ( 5.4) can be solved in a least-square setup, but also an acceptable and fast solution can be written as:

R 1 = F 1 ; C 1 = F 1 F T 1 U 1 (5.5)
This technique apparently works for linear models only. However, it can be easily extended into non-linear models by defining K LR per given neighborhood of the domain for example. The method is used in several linear and non linear cases. Interested readers can refer to [START_REF] Ibanez | On the data-driven modeling of reactive extrusion[END_REF][START_REF] Reille | Incremental dynamic mode decomposition: A reduced-model learner operating at the low-data limit[END_REF].

GENERIC formalism

One of the main drawbacks of classical fitting techniques (regression, neural networks...) is the potential violations of physical laws, as no physical constraint is applied in such a process. Entropy production, mass and energy conservation... are basic laws of thermodynamics that can be violated when using classical data-driven modeling techniques. One of the solution is the GENERIC formalism. GENERIC stands for "general equation for non-equilibrium reversible-irreversible coupling ". The method was originally established to model the rheological behavior of complex fluids [START_REF] Grmela | Bracket formulation of diffusion-convection equations[END_REF][START_REF] Grmela | Thermodynamics of driven systems[END_REF], and later on extended into a larger general framework for modeling material behavior [START_REF] Onate | Computation of the stabilization parameter for the finite element solution of the advective-diffusive problems[END_REF][START_REF] Ottinger | Beyond equilibrium thermodynamics[END_REF].

The method ensures a solution satisfying the basic thermodynamics laws through the fitting of several differential equations and establishing constraints on the fitted parameters [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF][START_REF] Gonzalez | Learning corrections for hyperelastic models from data[END_REF][START_REF] Gonzalez | Thermodynamically consistent datadriven computational mechanics[END_REF][START_REF] Gonzalez | Thermodynamically consistent data-driven computational mechanics[END_REF]. For instance, the general form of any problem is written in GENERIC formalism as:

ẋt = L(x t )∇E(x t ) + M(x t )∇S(x t ), (5.6)
supplemented with the following initial conditions:

x(0) = x 0 (5.7)
Now we add the degeneracy conditions relative to the GENERIC formalism [START_REF] Ottinger | Beyond equilibrium thermodynamics[END_REF]:

L • (∇S) = 0 M • (∇E) = 0 (5.8)
The basic thermodynamics constraints are satisfied by the following constraints:

Ė(x t ) = 0 Ṡ(x t ) ≥ 0 (5.9)
Now considering L a skew-symmetric operator and M a symmetric, positive semi-defined operator we can ensure the conservation of energy [START_REF] Ottinger | Beyond equilibrium thermodynamics[END_REF]:

Ė = ∇E • ẋ = ∇E • L∇E + ∇E • M∇S = 0 (5.10)
Moreover, the satisfaction of the second law of thermodynamics can be achieved through considering [START_REF] Ottinger | Beyond equilibrium thermodynamics[END_REF]:

Ṡ = ∇S • ẋ = ∇S • L∇E + ∇S • M∇S ≥ 0 (5.11)
Once the constraints set, we can start a discreet identification process. First of all, we consider the discreet form of the variable x such as:

x n+1 -x n ∆t = L(x n+1 )DE(x n+1 ) + M(x n+1 )DS(x n+1 ) (5.
12)

The matrices L, M, DE and DS can be identified using any minimization algorithm, for instance a gradient decent algorithm would yield an acceptable solution.

Sometimes L can be easily identified for a given problem, the unknowns are therefore U = (M, DE, DS). The cost function to minimize can be written in general as:

U = argmin (x (U) -x meas ) (5.13)
The selection of the problem variables x is not always trivial. Interested readers can refer to [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF] for an application of data-driven modeling using the GENERIC formalism, with an application on a poro-viscoelastic material. In that case, the selected set of problem variables x was given by:

x n =    w n U n σ n    (5.14)
x n being the material coordinate, w n the indentation depth and U n the indentation velocity. Considering the discreet form of the problem, we can write:

DE n = A • x n DS n = B • x n (5.15)
The minimization problem becomes therefore:

     U = (M, A, B) = argmin (x (U) -x meas ) L • Bx = 0 M • Ax = 0 (5.16)
The minimization problem is solved using an interior-point algorithm. The error between the found solution and the experimental measurement did not exceed 2.24% in the worst indentation case [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF].

A digital twin for correcting and enhancing existing models

Data-driven techniques are increasingly popular. Data is becoming more abundant and both the improvements in computational techniques and power have allowed scientist to overcome classical computation problems. However, in engineering, data is still expensive to collect in many fields and experiments, and is not as abundant as one would wish for. Moreover, in engineering many models are already very well performing in reproducing reality, with high fidelity in several situation. Such models are the fruit of centuries of research and should not be discarded easily [47].

Based on the aforementioned arguments, the possibility of combining the previously well established models, to a data-driven error modeling appears to be an appealing approach. The well established models would therefore predict a physical behavior, while the model ignorance can be fitted or evaluated using a data-driven technique. Such approach would combine the possibility of certifying classical models for aeronautic and transportation use for example, while enhancing them with data-driven correction terms, when the fitting is not very accurate for example.

The modeling effort starts by creating a classical simulation, and then enhancing it using a data-driven technique. For example, in [START_REF] Ghnatios | A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials[END_REF], a biphasic material is indented and modeled using a combination of a solid reaction force and a fluid reaction one [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]:

1 S • A P dA + F solid S = k (z(t) -w(t)) S (5.17)
S is the indented surface, A is the surface of contact with the indenter, P the pressure fieds in the fluid, F solid is the solid reaction force, k the indenter stiffness, z(t) the displacement of the base of the indenter and w(t) the indentation depth. To identify material properties, the permeability K for example, the pressure field P is required. The model identifying the permeability is enhanced using a neural network, fitting the error of the model as illustrated in figure 5.2.

To train the network, the permeability is required as the neural network target is the error E of the model defined by:

E = F exp -F solid -F f luid (5.18)
with F exp is the experimentally measured force and F f luid = A P dA. Noting that Y is an estimation of the total experimental force F exp , while Y 1 is an estimation of the error E, we can write:

Y = Y 1 + F solid + F f luid (5.19)
Considering figure 5.2, one can clearly see that every iteration of the training algorithm requires an update of the permeability value. Therefore, the use of the PGD to simulate the pressure fields becomes the best approach to take in this situation. The pressure is therefore simulated "a priori "before taking any optimisation step. We obtain the pressure therefore in the following form:

P = i=N i=1 F i (l, m) • N i (n) • H i (H n ) • R i (R n ) • U i (U n ) • K i (K n ) (5.20)
Interested readers can refer to [START_REF] Ghnatios | A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials[END_REF] of to section 2 for more details.

A digital twin for incorporating the experimental variability into existing models

Digital twins are very useful for correcting solutions indeed. However, in some situations, a ground truth can't be established or doesn't exist. This is the case for example when experimental variability is high, or when a large number of experiments lead to fluctuating results [START_REF] Ghnatios | A nonparametric probabilistic method to enhance pgd solutions with data-driven approach, application to the automated tape placement process[END_REF]. In this case, enhancing the solutions with data-driven modeling of the error would not be the best approach to take, as the target value is not clearly set. Instead, incorporating the experimental variability into the simulation becomes a suitable approach, defining therefore a margin of error for the simulated results.

Non-parametric probabilistic method review

The non-parametric probabilistic method (NPM) is one of the suitable results for the incorporation of the data variability into the simulated solutions. The method starts from a reduced basis defined through the use of POD-like methods and defines a probabilistic reduced basis, enhanced with variability based on the maximum entropy theory [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. The method also creates a reduced basis varying on the tangent manifold of the aforementioned initial reduced basis, conserving therefore mathematical properties like orthogonality for example. The NPM algorithm builds a stochastic reduced basis [W], by computing a variation of the original reduced basis [V ] through:

                     [W] = ([V ] + s[Z])[H s (Z)] [H s (Z)] = ([I n ] + s 2 [Z] T [Z]) -1/2 [Z] = [A] -[V ][D] [D] = ([V ] T [A]+[A] T [V ]) 2 [A] = [U] -[B]([B] T [U] [U] = [G(β)][σ]
(5.21)

[G(β)] is the only stochastic entity in the algorithm, built using the maximum entropy theory [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. The problem's parameters are identified using (s, β, [σ]) = argmin(J). The parameters are fitted by minimizing the cost function J defined by:

J = w E(Y ) -E( Ŷ ) + (1 -w) σ(Y ) -σ( Ŷ ) (5.22)
Eventually the minimization of the NPM method requires a large amount of direct solutions of the differential equation and therefore, even with the use of POD-like model reduction techniques, the enrichment of the simulation with experimental variability requires a large computing facility for parallel computations.

NPM applied for the PGD

The PGD solutions are computed in a separated form such as [45]:

y(q 1 ; q 2 ; • • • ; q D ) = i=N i=1 j=D j=1 Q ij (q j ), (5.23)
with q j the separated coordinates of the problem, and Q ij (q j ) the recurrent solutions in the domain q j . In a discrete form,

Q ij becomes a matrix [Q] j ∈ R[N j , N ],
with N j the number of degrees of freedom in q j domain, and N the number of products of functions required to converge the PGD solution of the problem. Using NPM, we choose to enrich some (or all) of the PGD reduced coordinates solution using the same algorithm depicted in equation (C.58). The results for a reduced basis of a domain [q j ] results in the following:

                       [U] = [G(β)][σ] [A] = [U] -[B]([B] T [U]) [D] = [Q] T j [M][A] + [A] T [M][Q] j /2 [Z] = [A] -[Q] j [D] [H] = [I n ] + s 2 [Z] T [M][Z] -1/2 [Q j ] = ([Q] j + s[Z]) [H] (5.24)
The algorithm illustrated in equation (5.24) represents an enrichment of one dimension in the PGD solution domain. The same algorithm can be designed to enrich as much dimensions as required, by generating the random basis [G(β)] as a higher dimensionality array. The number of hyper-parameters of the problem with increase accordingly.

In the optimization process, evaluating the quantities of interests O([Q], t) do not require any solution of the PDE problem, as the PGD defines a surrogate model replacing the PDE with a product of functions defined in every dimension of the domain. For instance, the evaluation of the quantities of interest for any solution requires only the knowledge of the solution y defined by: y(q 1 ; q

2 ; • • • ; q D ) = i=N i=1 j∈D d Q ij (q j ) j∈Dp Q ij (q j ), (5.25)
where D d are the deterministic, non-enriched basis, while D p are the probabilistic basis, enriched using NPM. Therefore, one solution and thus one evaluation of the quantities of interests, is computed only by using the product and sum of matrices. the method is applied in [START_REF] Ghnatios | A nonparametric probabilistic method to enhance pgd solutions with data-driven approach, application to the automated tape placement process[END_REF] to enhance the simulation of the automated tape placement in winding by using experimental measurements and their variability. The experimental data was available at several points of the domain, and has lead to the construction of the variability and confidence intervals in all the domain. Interested readers can refer to [START_REF] Ghnatios | A nonparametric probabilistic method to enhance pgd solutions with data-driven approach, application to the automated tape placement process[END_REF] or to section 3.

Chapter conclusion

In this chapter we reviewed several technologies for facilitating and creating a digital twin and data-driven modeling. Each technology is suitable for a given set of data application. Other techniques are previously used, including for example the GENERIC formulation creating a PDE modeling satisfying thermodynamics constraints such as the conservation of energy, the generation of entropy... The modeling is fully data-driven and is able to satisfy automatically basic engineering constraints [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF][START_REF] Gonzalez | Thermodynamically consistent datadriven computational mechanics[END_REF][START_REF] González | Learning corrections for hyperelastic models from data[END_REF][START_REF] Grmela | Thermodynamics of driven systems[END_REF].

Other methods are created, based on optimized regressions and incremental mode decomposition learners. These methods can easily operates on low data limits as they built the model incrementally, step by step, through the addition of terms to correct previously established models.

Finally, digital twins are introduced based on the possibilities offered by model reduction techniques, especially the PGD. Having parametric solutions allows the use of gradient decent algorithms on the fly, without the need to compute any additional solution of the partial differential equation as illustrated in sections 3 and 4.

In the current projects, an autoencoder technique is being developed to enhance the NPM method and reduced the number of variables to identify in the optimization algorithm. The method boosts therefore the convergence speed by about one order of magnitude.

Another ongoing plan is to extend the use of the GENERIC formalism to create digital twins corrections, correcting therefore the simulated results, with data modeling preserving physical laws. Many more applications are also awaiting these novel techniques to be applied, in composite manufacturing, materials processing, characterization, residual stress estimations...

Chapter 6

General conclusion

In this HDR proposal I reviewed the main research axis and achievements of my post-doctorate work. The main achievements for each axis of research are resumed in what follows:

Development of reduced order models tools and methods

In this part of my research, several achievements are performed with much impact in the use of model reduction techniques in the scientific community:

• A stabilized mixed formulation is derived and applied for several example benchmark problems. The method is being extended for flow with high Reynolds numbers, as previous applications revolved around flows with low Reynolds numbers.

• A physics based morphing of space is used to allow the separation of apparently nonseparable domains in the PGD framework. This morphing of space results in faster and more accurate solutions when separating material properties using classical singular value decomposition.

• A non intrusive PGD approach coupled to space morphing allowing an interpolation in apparently non-interpolable geometries, like the propagation of injection front in resin injection molding for example.

• Improvement of the PGD physical space decomposition (3D to a sequence of 2D/1D for example) using machine learning algorithm. The algorithm learns the 1D dimension functions from previously computed solutions, then applies it to new problems. Therefore, the complexity of a full 3D problem is reduced to only one 2D problem [START_REF] Ghnatios | Artificial intelligence based space reduction of structural models[END_REF].

In this axis of research, the future developments with a high potential of applications in the industrial world, especially with the digitization era, industry 4.0 and the industry of the future framework appears as follows:

• The development of the stabilized mixed formulation method for PGD in flows with high Reynolds numbers and transient solutions.

• The improvement of physics-based morphing of space into the PGD framework to become an automatic/learnable feature using machine learning algorithm based on Natrual Language Processing (NLP) techniques for example.

• Develop 1D dimension learners for different physics as previous applications revolved around thermal fields only.

Applications in the materials and processes field

This part is extremely rich with different state of the art applications. Several applications are illustrated in this work, however much more are achieved for different materials and processes:

• A creation of a permeability identification method based on PGD parametric simulation, coupled with model learners and optimization algorithms

• An advanced modeling and simulation of the 3D printing of metallic components using PGD and parametric path. The results can be used to optimize the deposition path based as they offer a real-time thermal field and deformation field in the final printed part. Based on the quantities of interest, a gradient decent algorithm can be used to optimize the deposition path as well as the process parameters.

• A modeling and simulation of the microwave curing of continuous fiber reinforced plastics. The modeling used heterogeneous material properties, an extremely refined mesh and the in-plane-out-of-plane decomposition to allow realistic computation times. The simulated electric fields are leveraged into a thermal heat generation term, and therefore the thermal fields are obtained. The results are a unique state of the art solution that nothing similar is produced until now [16].

• Advanced modeling of the compression of composite materials including the squeeze flow of plates, strands and inserts, Ericksen flow formulation of the squeeze flow, closing of preexisting voids using a one-shot continuous simulation... All these topics are a state of the art results never previously shown, used or produced. The results reproduce experimental observations never previously illustrated including for example the rotation of unidirectional composite laminae during compression. Also recent works involved the compression molding of hybrid materials including aluminum components along with unidirectional continuous fiber reinforced laminae.

• The modeling and simulation of real-human knee joint cartilage and meniscus during different types of jumping with experimental measurements. The results consists a first attempt ever to model and simulate the knee joint during impact [START_REF] Ghnatios | Knee joint injury risk assessment bymeans of experimental measurements and proper generalized decomposition[END_REF].

The developments in this axis of research are endless. All processes and material characterization techniques will make use of the simulation capabilities. Mainly the following points can be addressed:

• The improvement of material characterization, especially biphasic materials, by making use of advances achieved in axis 1 of this manuscript, as well as machine learning capabilities.

• Creation of advanced homogenization and quantities of interest "transfer function" approach for the improvement of the electromagnetic curing of composites.

• Improve the understanding of the squeezing of composite materials and void closing, leading to the optimization of the end results by controlling the process parameters.

• Gaining extensive knowledge of human tissue properties, which will lead to better modeling and manufacturing of effective replacements.

Technologies for data-driven and digital twins

This part novel and is an active research topic worldwide. Several achievements are highlighted below, all achieved in the recent few years:

• A creation of a statistically validated and optimized regression technique based on the use of meta-heuristic optimization algorithms to select only the most relevant features to a regression.

• A digital twin correction technique is used to create a better understanding of the error involved in the modeling, using machine learning algorithms

• Using a non parametric probabilistic approach to incorporating the experimental variability into the simulation.

• Improving the Non parametric Probabilistic Method reducing the computation time by one order of magnitude while improving the result. This achievement is performed by using convolution autoencoders based on convolution neural network layers.

• Creating a data-driven modeling technique based on the GENERIC formalism to achieve full satisfaction of thermodynamic constraints, while using only data-driven formalism.

This third axis of research is novel and therefore the development horizon is far from sight. Some ongoing developments are :

• The incorporation of physical constraints in the selection of optimized regression models • improving the Non parametric Probabilistic Method to achieve a final result within one minute, and therefore improve its feasibility and its ability to attract funds in the industrial framework and the industry 4.0 applications.

• Couple the GENERIC formalism to previous existing models to generate digital twins satisfying realistic thermodynamic constraints.

Next steps

As a general point of view, the research plan and future perspectives presented in this work are specific and the recruitment of the right persons to move forward consists a real challenge not to be underestimated. Moreover, funding is mandatory to achieve the aforementioned points and research goals.

Definitely, existing national and international collaborations are to be consolidated and novel ones are to be established. Industrial collaborations are of utmost importance to successfully achieve the proposed goals and work plan. However, as the industrial world is moving inevitably into digitization and industry 4.0 framework, I am confident that funding opportunities will be widely available in this proposed research area, for all the 3 axis of research, axis 1, 2 and 3. Moreover, composite and advanced materials manufacturing are in high demand, as aeronautical and transportation industry seek high performance at low weight. On the other hand, drones are coming at high speed as potential replacement of cars, and will require novel material, novel artificial intelligence based control and novel modeling. The current funding origins will be consolidated and more are yet to come, especially in these tackled emerging fields.

All proposed research projects must include a fast exploration tool, a creation of an industrial digital tool for real-time use, based on digital twins and advanced simulations, as well as an enrichment of the real-world manufacturing scenario, leading to the ability to optimize materials and processes in real-time, by minimally educated workers. The proposed work will therefore take all the knowledge and power of advanced modeling, artificial intelligence and simulation, to put it at hand of everyone in need, in real-time.

Introduction

The proper Generalized Decomposition or PGD is one of the MOR techniques with the particularity to be an "a priori "model order reduction technique. In fact, it requires no prior knowledge of the problem's eigenvectors to generate the results. Instead, an approximation of the eigenvectors is computed on the fly, through solving problems in lower dimensionality spaces [45]. The PGD is currently used in different computational mechanics fields. In fact, it made its proofs in heat transfer [START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF], elasticity [23], even in coupled multidimensional and multi-physical problems as addressed in [START_REF] Ghnatios | Computational vamedecum of the coupled thermal/mechanical behavior of composite materials during ultrasonic welding[END_REF]. When it comes to fluid mechanics, the PGD is used with the penalty formulation to solve the Stokes equation in different composite materials manufacturing related applications, the squeeze flow and Ericksen flow for instance [START_REF] Ghnatios | 3d modeling of squeeze flow of multiaxial laminates[END_REF]. Some other works tackled the Navier-Stokes equation and Rayleigeh Bernard non linear fluid flow, using different approximations [5,[START_REF] Dumon | Proper generalized decomposition (pgd) for the resolution of navier-stokes equations[END_REF] to avoid the tricky mixed formulation. In fact, the simulation of incompressible flow is known to be complicated since it requires to deal with subgrid phenomena and the imposing of the conservation of mass equation. The problem becomes even more complicated when dealing with separated representation, specially encountered in the model order reduction techniques. Recent works solved the Ericksen flow with mixed formulation using the PGD, with an attempt to define an LBB condition in the separated representation framework by increasing the degree of the interpolation functions of the velocity with respect to the pressure ones [START_REF] Ibanez | Simulating squeeze flows in multiaxial laminates: Towards fully 3d mixed formulations[END_REF]. To the best knowledge of the authors, no stabilized mixed formulation, using the same interpolation functions for both the velocity and the pressure fields, is implemented within the PGD framework.

We propose a new proper generalized decomposition for mixed formulation. Indeed, we refer to the use of Galerkin least Squares method (GLS) to deal with equal order velocitypressure interpolation [START_REF] Donea | Finite element method for flow problems[END_REF]. The Galerkin Least Squares method is successfully used to derive stabilized finite elements formulations in different areas [19,[START_REF] Coupez | Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF][START_REF] Hughes | The variational multiscale method -a paradigm for computational mechanics[END_REF][START_REF] Hughes | A multiscale discontinuous galerkin method with the computational structure of a continuous galerkin method[END_REF]. In this work, we derive first of all the governing equations related to the GLS method. Later on, we apply the PGD iterative algorithm to solve the resulting problem as a sequence of solutions in lower dimensionality domains. Finally, we show some numerical applications on some benchmark problems. Finally, we show the solutions of multidimensional problems including material, process and geometrical parameters as extra coordinates of the problem.

Governing equations

In this section, we illustrated the governing equation derivation used in this paper. For sake of simplicity, we propose to deal with the Stokes flow problem, as a general representation for mixed formulation. The resultant equation is written in a domain Ω by:

-µ∇ 2 v + ∇P = f ∇ • v = 0 (A.1)
where v is the velocity field in the domain Ω, P the kinematic pressure, f the applied body force and µ the dynamic viscosity. For the rest of the derivation, to simplify the notations and without any loss of generality, we use Dirichelet boundary conditions depicted on the domain boundary Γ = ∂Ω, where v = v Γ on the domain's boundary Γ.

The weak form of the problem depicted in equation (A.1) is written as:

Ω 2µ • ∇ s w : ∇ s vdΩ -Ω P ∇ • wdΩ = Ω w • f dΩ Ω Q∇ • vdΩ = 0 (A.2)
With ∇ s being the symmetrical component of the gradient tensor:

∇ s = 1 2 ∇ + ∇ T (A.3)
Now using the following notations for the bilinear forms:

a(w, v) = 2µ • Ω ∇ s w : ∇ s vdΩ b(Q, v) = -Ω Q∇ • vdΩ (A.4)
and using the linear form notation:

l(w, f ) = Ω w • f dΩ (A.5)
we can easily rewrite the weak form depicted in equation (A.2) by:

a(w, v) + b(P, w) = l(w, f ) b(Q, v) = 0 (A.6)
Equation (A.6) convergence under two conditions. First of all the bilinear form a should be coercive in the whole velocity space [34]. Second, the well-known LBB of infsup condition on the operator b should be addressed. The first condition is automatically satisfied, while the second condition is not satisfied using any random finite element interpolation. It is well known that low order elements do not satisfy the LBB condition [START_REF] Donea | Finite element method for flow problems[END_REF]. Many methods exist to satisfy these conditions [START_REF] Norburn | Stabiliezd vs stable mixed methods for incompressible flow[END_REF]. In this work, we are particularly interested in the Galerkin Least Squares method to satisfy the LBB condition. In the Galerkin Least Squares, we consider the following functional [START_REF] Franca | Convergence analysis of galerkin-least squares methods for symmetric advective-diffusive form of the stockes and incompressible navier-stockes analysis[END_REF]:

L(Q, w) = -µ∇ 2 w + ∇Q -f , -µ∇ 2 w + ∇Q -f (A.7)
Obviously, the argument minimizing the functional L solves the strong form of the Stokes problem illustrated in equation (A.1). Minimizing the first variational L in the element interior (to avoid adding extra continuity equations [START_REF] Liu | A variational multiscale stabilized finite element method for the stokes flow problem[END_REF]) leads to the following discreet weak form of the problem results after projecting on the domain mesh:

       a(w h , v h ) + b(P h , w h ) + n el e=1 γ e -µ∇ 2 w h , -µ∇ 2 v h + ∇P h -f Ωe = l(w h , f ) b(Q h , v h ) - n el e=1 γ e ∇Q h , -µ∇ 2 v h + ∇P h -f Ωe = 0 (A.8)
For low order interpolation functions, P1P1 for instance, the second derivation of the velocity components vanishes and the weak from (A.8) is reduced to:

   a(w h , v h ) + b(P h , w h ) = l(w h , f ) b(Q h , v h ) - n el e=1 γ e ∇Q h , ∇P h Ωe = - n el e=1 γ e ∇Q h , f Ωe (A.9)
The exact value of the parameter γ e can only be determined empirically [START_REF] Donea | Finite element method for flow problems[END_REF][START_REF] Onate | Computation of the stabilization parameter for the finite element solution of the advective-diffusive problems[END_REF]. For a thorough discussion on the Galerkin Least Squares approximation, one can refer to reference [START_REF] Franca | Convergence analysis of galerkin-least squares methods for symmetric advective-diffusive form of the stockes and incompressible navier-stockes analysis[END_REF] and its references therein.

The PGD algorithm in the mixed formulation framework

Considering the GLS stabilized mixed formulation developed in section 1.2, we develop the PGD algorithm corresponding to the solution of the mixed formulation. In the PGD framework, we use the domain decomposition, such as a 3D problem can be solved as a sequence of 2D × 1D problems or even a 1D × 1D × 1D. For the sake of simplicity, we are illustrating the domain decomposition in a 2D domain, where the solution is performed as a sequence of 1D × 1D problems. Let's consider the weak form obtained after GLS stabilization illustrated in equation (A.9) in a 2D prismatic domain (x, y). To solve the problem using the PGD, we suppose that the solution (v, P ) is written in a separated form such as:

v P ≈     N i=1 X v i (x) • Y v i (y) N i=1 X P i (x) • Y P i (y)     (A.10)
where • is the Hadamard product, or point wise term by term multiplication of the velocity components vectors X v i and Y v i . X P i and Y P i are the pressure P components. In general, such decomposition is possible even if the (x, y) domain is not prismatic [START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF][START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF]. To compute the solution as defined by a sequence in equation (A.10), we first of all suppose that the solution is know until an order n -1 such as n < N, and we search to enrich by a new product of vectors or a new "mode", such as:

v n P n =     n-1 i=1 X v i (x) • Y v i (y) + R(x) • S(y) n-1 i=1 X P i (x) • Y P i (y) + T (x) • U(y)     = v n-1 + R(x) • S(y) P n-1 + T (x) • U(y) (A.11)
Using the notation of equation (A.11) and writing the velocity v = (u, v), we can find the velocity gradient for example written as:

∇v = ∂u ∂x ∂u ∂y ∂v ∂x ∂v ∂y ≈ n-1 i=1 ∂X u i ∂x X u i ∂X v i ∂x X v i •   Y u i ∂Y u i ∂y Y v i ∂Y v i ∂y   + ∂R u ∂x R u ∂R v ∂x R v • S u ∂S u ∂y S v ∂S v ∂y = n-1 i=1 X v i • Y v i + R • S (A.12)
And the pressure gradient is written by:

∇P = ∂P ∂x ∂P ∂y ≈ n-1 i=1 ∂X P i ∂x X i •   Y i ∂Y P i ∂y   + ∂T ∂x T • U ∂U ∂y = n-1 i=1 X P i • Y P i + T • U (A.13)
The simplest choice of the test functions w and Q would be:

w = R * • S + R • S * Q = T * • U + T • U * (A.14)
where R * and T * are test functions that depends only on x only, while S * and U * are functions of y only. Replacing equations (A.12), (A.13) and (A.14) into the weak form of the problem (A.9) leads to a non-linear equation with four unknowns: R and S in the velocity field on one side, T and U in the pressure field on the other side. The resulting problem is therefore a non-linear one, requiring a linearization scheme. The simplest choice is a fixed point iterative algorithm where we first consider S and U as known, and we compute R and T . Later on, we use the just computed R and T to compute S and U. Then, the iterations proceeds until the convergence of the four unknown functions. One may note that the convergence of the PGD algorithm reaches a solution u tensor different than the finite element one unless satisfying some conditions as highlited in [START_REF] Hamidi | On the convergence of altrnating minimization methods in variational pgd[END_REF], non the less the sabilization term γ e in this work is the same used for the finite element stabilization method. The algorithm is illustrated in what follows.

Computing R and T

In this section, we consider S and U as known and we compute the functions R and T . Since S and U are known we choose the test functions S * and U * as equal to zero, resulting in the following test functions:

w = R * • S Q = T * • U (A.15)
Replacing equations (A.12), (A.13) and (A.15) into the weak form (A.9) of the problem, leads to the following bilinear form a:

a(w, v) = Ω ∇ s w : ∇ s vdΩ = 1 4 Ω ∇w : ∇v + ∇w : ∇v T + ∇w T : ∇v T + ∇w : ∇v T dΩ = 1 4 Ω (R * • S) : n-1 i=1 X v i • Y v i + R • S + (R * • S) : n-1 i=1 X v i • Y v i + R • S T + (R * • S) T : n-1 i=1 X v i • Y v i + R • S + (R * • S) T : n-1 i=1 X v i • Y v i + R • S T dΩ (A.16)
While the terms b(•, •) and l(•, •) of equation (A.9) are written by:

b(Q h , v h ) = Ω Q • ∇ • vdΩ = Ω Q • ∂u ∂x + ∂v ∂y dΩ = = Ω (T * • U) • i=n-1 i=1 ∂X u i ∂x • Y u i + ∂R u ∂x • S u + i=n-1 i=1 X v i • ∂Y v i ∂y + R v • ∂S v ∂y dΩ (A.17) l(w h , f ) = Ω w • f dΩ = Ω (R * • f x • S • f y ) dΩ (A.18)
f x and f y being the x and y components of the second hand side vector f . Finally the stabilization terms of equation (A.9) can be deveopped in a separated form such as:

n el e=1 γ e ∇Q h , ∇P h Ωe = = n el e=1 γ e ∂T * ∂x • U • i=n-1 i=1 ∂X P i ∂x • Y P i + ∂T * ∂x • U • ∂T ∂x • U+ +T * • ∂U ∂y • i=n-1 i=1 X P i • ∂Y P i ∂y + T * • ∂U ∂y • T • ∂U ∂y Ωe (A.19)
and:

n el e=1 γ e ∇Q h , f Ωe = n el e=1 γ e ∂T * ∂x • U • f x + T * • ∂U ∂y • f y Ωe (A.20)
The final form obtained from replacing the separated form operators into the weak form (A.9) consists of a 1D equation where the only unknowns are R and T which depends only on x. Interested reader may refer to [START_REF] Ghnatios | 3d modeling of squeeze flow of multiaxial laminates[END_REF] and its references therein for more details.

Computing S and U

After computing R and T , we now consider these functions as known and compute the functions of y (S and U) from the already computed R and T . Therefore we choose the test functions as:

w = R • S * Q = T • U * (A.21)
Replacing into the weak form of the problem (A.9) leads again to a 1D problem which can be project on the mesh in the y direction and solved using any classical numerical technique.

Continue the iterative process

After obtaining the (S, U) values, we repeat the steps detailed in section 1.3, then 1.3 until the convergence of the unknowns (R, S, T, U). Once this convergence achieved, we can add another product of functions until the convergence of the residual of the differential equation of the problem [45]. The convergence is controled by the precentage of reduction in the residual R of the differential equation at an iteration N:

R n=N R n=1 < e (A. 22
)
Where e is the error tolerance set to 10 -4 in general [START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF], and R n=1 is the initial residual before starting to add the PGD product of functions to the soluton.

Numerical examples

In this section, we show some numerical examples solved by applying the PGD algorithm along with the GLS stabilization method. First we show some benchmark problems along with their solutions. Later on, we tackle more complex problems with additional parameters rendering the problems multidimensional and intractable using classical simulation techniques.

A lid-driven cavity

In this section, we use the PGD formulation illustrated in section 1.3 along with the GLS stabilization method to solve the Stokes flow problem in a lid-driven cavity. The geometrical domain is illustrated in figure A.1. If the velocity field in the domain is written by v = (u; v) and the pressure by P , the boundary conditions of the simulated domain are resumed by:

         u = v = 0 at x = 0 and x = L x u = v = 0 at y = 0 u = u 0 at y = L y v = 0 at y = L y (A.23) u = v = 0 u = v = 0 u = u 0 L x L y u = v = 0 x y Figure A.1:
The simulated lid driven cavity Using the PGD algorithm, the physical domains x and y are separated such as the obtained solution is written as

v P =         u = n i=1 X u i (x) • Y u i (y) v = n i=1 X v i (x) • Y v i (y) P = n i=1 X P i (x) • Y P i (y)         (A.24)
In this example, we choose L x = L y = 1m and the viscosity µ = 1P a.s The mesh of the x and y domains consists each of 1000 nodes uniformly distributed with linear shape functions. The same mesh and linear shape functions are used for both the pressure and the velocity fields. The simulation is performed on a core i7 normal portable PC. The convergence is achieved within 31 seconds for an equivalent of 1000 × 1000 × 3 degrees of freedom in the domain, with a total number of computed products of vectors n = 25. Figures A.2, A.3, A.4 illustrate respectively the velocity components u and v as well as the pressure P in the cavity. Figure A.5 illustrates the velocity vector fields in the simulated domain.

The residual reduction factor is shown in figure A. 6. In this figure we can clearly note that the first 10 to 11 product of vectors are the most relevant to reduce the residual error. It would have been wise to stop the iterations at 11 computed product of vectors which could save about 70% of the calculation time without much influence on the final solution. One may also define the error using the final converged solution such as:

e P GD = v P N - v P n v P N (A.25)
The error e P GD is shown in figure A.7 as a function of the added PGD products of functions n. We can clearly see that the first 10 to 11 modes are the most relevant in the solution.

v P =              u = N i=1 U 1 i (x, y) • U 2 i (z) v = N i=1 V 1 i (x, y) • V 2 i (z) w = N i=1 W 1 i (x, y) • W 2 i (z) P = N i=1 P 1 i (x, y) • P 2 i (z)              (A.28)
The test function is also written in a separated form such as:

V * P * = V * xy (x, y) • V z (z) + V xy (x, y) • V * z (z) P * xy (x, y) • P z (z) + P xy (x, y) • P * z (z) (A.29)
The problem is solved using an iterative fixed point algorithm, considering on a first stage (V * xy ; P * xy ) as known equal to zero and computing the functions of z, (V z ; P z ) for instance. On a second stage, the functions of z are known and we compute the in-plane functions of xy. The adopted mesh consists of 3721 nodes in the (x, y) 2D domain and 1000 nodes in the z domain. The resultant is a 3D solution with 14 884 000 degrees of freedom. The solution is computed within 144.8 seconds on a normal portable PC. The solution of the 3D velocity vectors is illustrated in figure A.12 for a viscosity µ = 1P a.s. The velocity components u, v and w at z = h/2 are illustrated in figures A.13, A.14 and A.15 respectively. Moreover, the pressure P at a section normal to z taken at z = h/2 is illustrated in figure A. 16.

Introducing process and material parameters as extra coordinates

Often the process parameters are unknown or should be optimized. Therefore, disposing of the solution for all possible process parameters may be an appealing approach to reduce the optimization time [6]. In this section, we illustrate the solution of the squeeze flow using a In the formulation illustrated in section 1.3, the stabilization coefficient depends on the process parameters. Therefore, this term depends on the process coordinates in this problem. In this section, we seek the solution of the problem as a function of the problem's viscosity µ and the imposed squeeze flow compression velocity U as extra coordinates of the problem. Thus, the solution has the following form:

v P =     N i=1 X v i (x, y) • Z v i (z) • M i (µ) • U i (U) N i=1 X P i (x, y) • Z P i (z) • M i (µ) • U i (U)     (A.30)
The problem is solved using the following boundary conditions:

     v = (0; 0; 0) at z = 0 v = (0; 0; -U)m/s at z = h P = 0 elsewhere (A.31)
The problem is now a 5D problem, the mesh consists of 3721 nodes in the in-plane (x, y) domain, 10 000 nodes in the z domain, 100 nodes in the viscosity µ mesh and 100 nodes in the imposed velocity U mesh. The solution contains the equivalent of 10 4 × 3D problems, each 3D problem having 148 840 000 degrees of freedom. The solution of the 5D problem is obtained within 320 seconds on a normal portable PC, core i7 with 8GB Ram. Such solution is even unimaginable using classical calculation techniques. One may also wish to illustrate the velocity fields at a given point as a function of the process parameters for example. Therefore we can fix the in-plane coordinate at a given node for example and illustrate, as shown in figure A.21, the results as a function of the viscosity µ and the imposed velocity U for example.

Introducing process, material and geometrical parameters as extra coordinates

In this section, we illustrate the solution of the squeeze flow using the process and material parameters as extra coordinates of the problem. Therefore, we consider the thickness of the domain h, its width t and the process parameters µ and U as extra parameters of the problem.

The solution leads to a 7D abacus where the velocity and pressure are a function of the 3 physical coordinates (x, y, z) and the 4 extra coordinates added to the problem. Thus, we are seeking a solution of the form:

v P =     N i=1 X v i (x, y) • Z v i (z) • M i (µ) • U i (U) • H i (h) • T i (t) N i=1 X P i (x, y) • Z P i (z) • M i (µ) • U i (U) • H i (h) • T i (t)     (A.32)
With the following boundary conditions:

     v = (0; 0; 0) at z = 0 v = (0; 0; -U)m/s at z = h P = 0 elsewhere (A.33)
The absence of the extra geometrical parameters from the weak form of the problem, or the boundary conditions adds an extra layer of complication to the problem. In fact, to introduce the extra geometrical parameters into the weak form of the problem, a geometrical mapping is mandatory [START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF][START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF].

The resultant 7D problem is solved in a plane 2D domain (x, y) containing 3721 nodes, while the 1D domains include 10 000 nodes on z, 100 nodes in each of U and µ and 300 nodes in each of the geometrical parameters h and t. The solution contains the equivalent of 9 × 10 8 × 3D solutions, each 3D solution having 148 840 000 degrees of freedom. The solution is obtained within 406.52 seconds on a normal portable PC, core i7 and 8GB Ram.

Conclusion

In this work, we illustrates the PGD solution of the Stokes problem using a full mixed formulation velocity-pressure. The PGD algorithm is used along with the stabilized Galerkin least squares weak formulation. The PGD development in the GLS framework is illustrated, later on some benchmark problems are solved including the lid cavity problem, Poiseuille flow and squeeze flow. Moreover, the PGD is used to solve parametric Stokes problems in the same GLS stabilized weak form framework. The results as well as the computation cost of each test case are illustrated.

As a conclusion, the PGD shows its appropriate use for mixed formulations; this opens the door for different new formulation as well as applications as illustrated the many proposed test cases.

Introduction

When looking for an approximation of the solution u(x, t) of a given PDE, here assumed, without loss of generality, to be scalar and linear, the standard finite element method considers the approximation

u(x, t) = N i=1 U i (t)N i (x), (A.34)
where U i denotes the value of the unknown field at node X i and N i (x) represents the so-called shape function associated to the i-th node. Here N refers to the number of nodes employed to approximate the field u in the domain Ω in which the physical problem is defined. This approximation results in an algebraic problem of size N in the linear case, or the iterative solution of many of them in the general transient and nonlinear case. In order to alleviate the computational cost, model order reduction techniques were proposed and are nowadays intensively used.

In the framework of POD-based model order reduction [48], a learning stage allows extracting the significant modes φ i (x) that best approximate the solution. In many cases, a reduced number of modes R (R ≪ N) suffice to approximate the solution of problems similar to the one that served to extract the modes in the learning stage.

Thus, the solution u(x, t) is projected onto the reduced basis composed of functions

{φ 1 (x), • • • , φ R (x)}, according to u(x, t) ≈ R i=1 ξ i (t)φ i (x), (A.35)
that now requires the solution of linear systems of size R instead the ones of size N characteristic of finite element solutions. This often allows us to obtain impressive computing time savings. Approximations (A.34) or (A.35) imply a finite sum of time-dependent coefficients and space functions. The last are assumed known: they consist of the usual finite element shape functions or the modes extracted by applying the POD. A generalization of this procedure could consist on assuming that space functions are also unknown and therefore to compute both, time and space functions, on the fly [START_REF] Ladevèze | The large time increment method for the analyze of structures with nonlinear constitutive relation described by internal variables[END_REF]. Thus, this approximation reads

u(x, t) ≈ M i=1 T i (t)X i (x).
(A.36)

Since both functions involved in approximation (A.36) are unknown, the problem solution becomes nonlinear, and its solution requires an appropriate linearization strategy. The interested reader can refer to [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF][START_REF] Cueto | Proper Generalized Decompositions: An Introduction to Computer Implementation with Matlab[END_REF] and the references therein for practical details on the computer implementation of separated representations.

Expression (A.36) evidences that the solution procedure requires the solution of about M problems, with M ≪ N and M ∼ R (slightly more, in fact, because of the nonlinearity induced by separated representations) involving the space coordinates (in general three-dimensional and whose associated discrete systems are of size N) to compute the space functions X i (x) and about M one-dimensional problems to calculate the time functions T i (t). Due to the fact that the computing cost related to the solution of 1D problems is negligible with respect to the solution of 3D problems, the resulting computational complexity reduces drastically, scaling with M instead of P (here, P represents the number of time steps considered in the time domain discretization and that corresponds with the complexity of standard incremental time-integration techniques).

Another step forward consisted in assuming model parameters as extra-coordinates. Thus, space-time-parameters separated representations allowed constructing the so-called computational vademecums (also known as abacus, virtual charts, nomograms, ...) efficiently employed for multiple purposes: simulation, optimization, inverse analysis, uncertainty propagation and simulation-based control, all them under real-time constraints [48,[START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF]. When the unknown field involves space, time and a series of parameters µ 1 , . . . , µ Q , its associated separated representation reads

u(x, t, µ 1 , . . . , µ Q ) ≈ M i=1 X i (x)T i (t) Q j=1 M j i (µ j ). (A.37)

Space separation

The separation of space coordinates was also intensively considered in our former works referred later. Space separation allowed addressing multi-physics problems defined in degenerated geometries in which at least one of its dimensions results to be much smaller that the other ones (e.g. beams, plates, shells, laminates) or processes involving additive layers (e.g., automated tape placement, 3D printing, or additive manufacturing). Thus, if domain Ω can be decomposed as Ω = Ω x × Ω y × Ω z , the solution u(x, y, z) could be approximated by using the separated representation

u(x, y, z) ≈ M i=1 X i (x)Y i (y)Z i (z), (A.38)
that allows calculating the 3D solution from a sequence of 1D problems. For some geometries, as the ones concerning plates or shells, in-plane-out-of-plane separated representations become specially appealing,

u(x, y, z) ≈ M i=1 X i (x, y)Z i (z), (A.39)
where the 3D complexity is reduced to a 2D complexity, related to the calculation of in-plane functions X i (x, y). As discussed in [48], the cost savings provided by the use of these separated representations are potentially impressive when the spatial domain is fully separable. Indeed, the complexity of the simulation now scales with the one-dimensional meshes used to solve the BVP's in Ω x , Ω y and Ω z , associated to the computation of functions X i (x), Y i (y) and Z i (z) or with the two-dimensional ones associated with the calculation of functions X i (x, y) in the case of in-plane-out-of-plane separated representations [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF].

When the domain is not intrinsically separable or, in other words, expressible from a direct cartesian product, fully separated representations require the use appropriate geometrical mappings or the immersion of the non-separable domain onto a fully separable one. The interested reader can refer to [START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF][START_REF] Gonzalez | Recent advances in the use of separated representations[END_REF] for additional details on both approaches.

In-plane/out-of-plane separated representations are particularly useful for addressing the solution of problems defined in plate and shell geometries, [24] and [28] respectively, or extruded domains [START_REF] Leygue | Towards a framework for non-linear thermal models in shell domains[END_REF]. A parametric 3D elastic solution of beams involved in frame structures was proposed in [30]. The same approach was extensively considered in structural plate and shell models in [START_REF] Gallimard | Coupling finite element and reliability analysis through proper generalized decomposition model reduction[END_REF][START_REF] Vidal | Composite beam finite element based on the proper generalized decomposition[END_REF][START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF][START_REF] Vidal | Explicit solutions for the modeling of laminated composite plates with arbitrary stacking sequences[END_REF][START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF][START_REF] Vidal | Assessment of variable separation for finite element modeling of free edge effect for composite plates[END_REF] and [START_REF] Pruliere | 3d simulation of laminated shell structures using the proper generalized decomposition[END_REF].

Space separated representations where enriched with discontinuous functions for representing cracks in [START_REF] Giner | The proper generalized decomposition (pgd) as a numerical procedure to solve 3d cracked plates in linear elastic fracture mechanics[END_REF], delamination [START_REF] Metoui | The proper generalized decomposition for the simulation of delamination using cohesive zone model[END_REF] and thermal contact resistances in [START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF]. The in-plane/out-of-plane decomposition was then extended to many other physics like squeeze flows of Newtonian and non-Newtonian fluids in layered domains in [START_REF] Ghnatios | 3d modeling of squeeze flows occuring in composite laminates[END_REF].

Outline review

The present paper aims to propose a new simple procedure for treating, in a fully separated representation manner, problems in which the expression of the parametric dependency in a separated form was considered difficult, or where the separated representation performed by invoking SVD-like techniques had a significant impact in the performance and efficiency of the resulting numerical procedures.

After the just addressed introduction on separated representations, with special emphasis in those involving space coordinates, the next section will introduce the main idea of the present work. Namely, a very efficient geometrical mapping leading to a separable description of the problem. Then, in Section 3, the proposed approach is applied to the solution of thermal problems defined in thin domains, but in which the material conductivity has not a compact separated representation. Section 4 addresses the issue of 1D thermal problems containing an inclusion of a different thermal conductivity and whose position in the domain is assumed to be the model parameter. Consequently, the temperature field for any position of the inclusion is found. The paper finishes by summarizing the most valuable conclusions.

Domain mapping

We introduce the proposed methodology, for the sake of simplicity and without loss of generality, by making use of a heat transfer problem. Thus, we consider the temperature field u(x), x = (x, y), defined in a thin domain Ω = (0, L) × (0, H), with H ≪ L, where a non-planar internal boundary (as the one sketched in Fig. A.28) separates the upper domain Ω u , with thermal conductivity K u , and the bottom domain, Ω b , with conductivity K b . This boundary is defined by a function h(x).

The temperature field is assumed to be the solution of the steady-state heat transfer problem defined in Ω, whose weak form reads

Ω K(x)∇u * • ∇u dx = 0, (A.40)
where u is the problem's variable such as u ∈ H 1 space and u * being the test function existing 94 in H 1 space. The problem is subjected to the boundary conditions

                           ∂u(x,y) ∂x x=0 = 0 ∂u(x,y) ∂x x=L = 0 u(x, 0) = 0, u(x, H) = u g , (A.41)
where u g denotes the prescribed, temperature on the top surface. are faced to the difficulties relative to the presence of the non-planar boundary, that implies a hardly separable conductivity field. The solution of u is considered converged once the residual of the differential equation is small enough. Thus, the representation illustrated in Eq. (A.42) involves a large M value to converge. Of course, conductivity can always be separated by invoking, for instance, the singular value decomposition (SVD), i.e., by expressing the conductivity as

Ωu Ωb

K(x, y) ≈ G k=1 F k (x)G k (y). (A.43)
In the case of a planar boundary Γ = Ω u ∩Ω b , characterized by h(x) = h, with 0 < h < H, a single term suffices for separating the conductivity, i.e.,

K p (x, y) = F p (x)G p (y), (A.44)
where the superscript (•) p refers to the planar interface configuration,

F p (x) = 1 and G p (x) is defined by G p (y) = K u -(K u -K b )χ(y), (A.45)
with

χ(y) = 1 y < h, 0 y ≥ h. (A.46)
However, when internal boundaries deviate from the planar configuration, the number of modes G increases prohibitively, and with it the operators involved in the weak form (A.40). Thus, both the performance and the efficiency of the solver are seriously compromised.

In order to circumvent this issue, we define two mappings transforming Ω u and Ω b into two rectangular domains respectively, R u and R b . (A.48)

The components of the Jacobian matrix J b result

           ∂x ∂r = 1,
B b =   1 -s h ′ (r) h(r) 0 1 h(r)   , (A.52)
and the gradient operator ∇ r = ( ∂ ∂r , ∂ ∂s ) T , such that

∇• = B b ∇ r •, (A.53)
that allows writing the weak form as

R b K b (∇ r u * ) T B T b B b ∇ r u det(J b )dr = 0. (A.54) 2.
Mapping Ω u into R u . Equivalently, for Ω u we define

x = r, y = (s -1)(H -h(r)) + h(r), (A.55) with (r, s) ∈ R u = (0, L) × (1, 2).
The gradient operator and the Jacobian are then calculated for this mapping following the same rationale considered in the previous approach. The final weak form in R u reads

Ru K u (∇ r u * ) T B T u B u ∇ r u det(J u )dr = 0, (A.56)
where B u is now given by

B u = 1 h ′ (s-2) H-h 0 1 H-h . (A.57)
Parabolic interface In this first numerical example we consider an extremely thin rectangular domain Ω = (0, 1) × (0, 10 -3 ), i.e., H = 10 -3 and L = 1, with the internal boundary h(x) defined by the parabolic curve When expressing the conductivity in a separated form according to Eq. (A.43) by invoking the singular value decomposition, G = 96 modes are needed for a reduction of six orders of magnitude of the associated eigenvalues. The solution using that G-mode separated representation of the thermal conductivity requires a non negligible effort, first for performing the SVD and then for constructing the separated representation of the temperature. In the last, more It is worth mentioning that the usual separated representation constructor (see [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]) that proceeds by computing rank-one updates, hardly converged, and each new mode produced only a very slight reduction of the residual.

h(x) = H 2 0.5 + 2x -x 2 . (A.

Wavy boundary In this second example the boundary h(x) is defined from

h(x) = H 4 sin (4xπ) + H 2 . (A.65)
Meshes and model parameters are the same that the ones considered in the previous case study. Despite of poor separability of the thermal conductivity field, again by using the mappings discussed in Section 2, the solution contained only M = 6 modes. It was obtained in 0.36 seconds, with a similar residual reduction factor ( 10 -9 ).

The solutions in the reference and real domains, R and Ω, are depicted respectively in Figs. A.32 and A.33. Again, the interface location is highlighted. When using a SVD-based separated representation of the conductivity field, the solver hardly converges, and every new mode reduced very slightly the residual. Consequently, the calculation required more than half an hour.

Three-dimensional layered domain with an inclined flat interface

In this section the procedure described in Section 2 is extended to the 3D case considering the domain sketched in

for s ∈ (0, 1)      x = r, y = s h(x, z), z = t, for s ∈ (1, 2)      x = r, y = (s -1) (H -h(x, z)) + h(x, z), z = t, (A.66)
where h(x, z) is given by

h(x, z) = (H/4) + (H/2) x L x . (A.67)
The thermal problem is solved by prescribing temperatures at the top and bottom surfaces, u(x, 0, z) = 0 and u(x, H, z) = 25, while enforcing a null heat flow on the lateral surfaces x = 0, x = L x , z = 0 and z = L z . The thermal conductivities of the subdomains located above and below the interface, K u and K b , are the same that the ones previously considered. The meshes considered in this simulation consisted of 100 nodes along the r-and tcoordinates, and 1000 along the thickness s, to reach a resolution equivalent to 10 7 finite element degrees of freedom. The most remarkable fact was that the calculation was accomplished in 3 seconds in a standard laptop, and that M = 15 modes were enough for reducing the residual in 8 orders of magnitude, that is, a reduction by factor of 10 -8 .

When expressing the conductivity using both, an in-plane/out-of-plane separated form or a fully separated form, by invoking in both cases the SVD, the usual rank-one constructor [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF] did not reach convergence after computing 100 modes (note that 15 modes were enough when using the procedure based on the mapping described in Section 2).

Parametric inclusion

In this last example we analyze an issue of major relevance, related to a thermal problem in a one-dimensional domain x ∈ [0, L], with L = 1, assumed to be composed by a material The parametric temperature field is sought. In other words, we look for the expression of the temperature field for any possible position of the inclusion X, i.e., u(x; X). Within the PGD rationale, the inclusion location X is considered as an extra-coordinate, so that the temperature field becomes two dimensional, i.e. u(x, X), while the complexity due to the increase of dimensionality is compensated by the use of the separated representation.

The parametric weak form of the heat equation now reads

Ω ∂u * ∂x K(x; X) ∂u ∂x dx = Ω u * Qdx, (A.68)
with Q the heat source term. We consider the mapping:

s ∈ [0, 1] X = r x = sh 1 (X) = sh 1 (r) s ∈ [1, 2] X = r x = (s -1)h g + h 1 (X) = (s -1)l + h 1 (r) s ∈ [2, 3] X = r x = (s -2)h 2 (X) + h g + h 1 (X) = (s -2)h 2 (r) + l + h 1 (r) (A.69)
where h 1 (X) is the position of the inclusion left boundary, i.e. h 1 (X) = Xl/2, l is the inclusion length and h 2 (X) the distance between the right boundary of the inclusion and the domain right boundary, i.e. By using the expression of the Jacobian, the heat problem weak form (A.68) can be expressed in the reference (r, s) domain and discretized without major difficulties, thus leading the parametric solution in the separated form

u(r, s) ≈ M i=1 R i (r)S i (s).
(A.70)

First numerical example In this first numerical example we consider Q = 50 with the thermal conductivity K(x; X) expressed form:

K(x; X) = 5 if x ∈ [X -l/2, X + l/2] 1 if x ∈ [0, X -l/2] ∪ [X + l/2, L] (A.71)
and homogeneous boundary conditions, u(x = 0) = 0 and u(x = L) = 0.

Meshes contained 1500 and 100 nodes along the s-and r-domains, respectively. The parametric solution u(r, s), and consequently its counterpart u(x, X), were obtained in less than 1s again by making use of a standard laptop, with a reduction factor of 30 with respect to the 100 1D calculations needed for computing the equivalent parametric temperature field, one for each of the 100 positions of the inclusion. Of course, when increasing the problem dimensionality (2D or 3D) the computing time savings are much higher. The constructed separated representation involves M = 10 modes. It is worth noting, with respect to the discussion addressed in previous sections, that if we try, as discussed in [START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF], to proceed by separating the thermal conductivity K(x, X) according to

K(x, X) ≈ G i=1 F i (x)G i (X), (A.72)
by invoking the SVD, the parametric solution u(x, X) directly constructed in the domain Ω involves 97 modes, M = 97 and requires 10 times more computing time.

Second numerical example

In order to better emphasizing the effect of the inclusion, we consider Q = 0, the thermal conductivity expressed now from

K(x; X) = 1 if x ∈ [X -l/2, X + l/2] 5 if x ∈ [0, X -l/2] ∪ [X + l/2, L] (A.73)

Thermal source moving along a one-dimensional domain

In this example, we consider L = 1m and the thermal source location X defined in X ∈ Ω = [0.035L, 0.965L]. X thus depends on the time from X = x 0 + νt. In this example, the thermal source velocity ν corresponds to a movement from X = 0.035L to X = 0.965L in 1s., leading to x 0 = 0.035 and ν = 0.

     h 1 (X) = XL -hg 2 , h g (X) = L 20 , h 2 (X) = L -hg 2 -XL. (A.76)
We define a new mapping:

For s ∈ [0, 1] X = x 0 + νt x = sh 1 (X) For s ∈ [1, 2] X = x 0 + νt x = (s -1)h g + h 1 (X) For s ∈ [2, 3] X = x 0 + νt x = (s -2)h 2 (X) + h g + h 1 (X) (A.77)
that allows transforming the problem weak form into the reference space (x, X). The same problem was solved using the standard PGD procedure. First the source term was expressed in a space-time separated form by invoking the SVD. Because of its poor separability, many modes were required, with the eigenvalues involved in the SVD decomposition decreasing very slowly. Then, the standard rank-one greedy algorithm computed the temperature solution. However, its poor separability required again more than 100 modes and more than 50 seconds calculation to find a solution that exhibits small oscillations.

Conclusions

In this paper we addressed the space separation in layered domains Ω where interfaces are not planar (or, even being planar, are inclined with respect to the in-plane-coordinate). In these circumstances, former works stressed that space-separated representations lose their expected effectiveness. In fact, that conclusion was obtained from the fact that these geometries involved too many terms in the material property-separated representations when invoking the SVD or its high order counterpart, the so-called HOSVD [START_REF] Kolda | Tensor decompositions and applications[END_REF]. In these circumstances, the separated representation constructor converges too slowly, with the subsequent impact on accuracy and computational efficiency.

In this work we proved that when the domain Ω is mapped into a fully separable one R standard separated constructors recover their efficiency, as the numerical results here described and discussed proved.

The technique here proposed opens numerous possibilities, such as the flow in thin ducts with rough surfaces, where fully 3D discretization fails, whereas the use of separated representations can lead to levels of resolution never attained before.

Spurious-free interpolations with non-intrusive PGDbased parametric solutions

Introduction

Computational abaci or vademecums can be built either through the use of numerical simulation or data-driven modeling [START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF][START_REF] Gonzalez | Computational Biomechanics for Medicine[END_REF][START_REF] Lu | Datadriven hopgd based computational vademecum for welding parameter identification[END_REF]. They consist of a compilation or a handbook of solutions precomputed "offline", to be used on the fly, "online". The possibility of performing real time simulation in an industrial framework at the expense of offline heavy computations is motivating the industrial companies to move in this direction. The development of such computational vademecums became possible thanks to the development of model order reduction techniques [30], especially the Proper Generalized Decomposition, PGD, which is probably a very suitable "a priori" model order reduction technique. However, PGD is by nature the most intrusive model reduction technique due to its "a priori" formulation. Recent efforts have been undertaken to formulate the PGD in a non-intrusive manner as much as possible [31,46,[START_REF] Zou | Pgd for constrained parametric space with a non-intrusive implementation[END_REF]. Some non-intrusive versions of the PGD rely on external in-house flow processing codes that controls validated open-source finite element (FE) codes or commercial software [START_REF] Leon | Non-intrusive proper generalized decomposition involving space and parameters: application to the mechanical modeling of 3d woven fabrics[END_REF][START_REF] Quaranta | A minimallyintrusive fully 3d separated plate formulation in computational structural mechanics[END_REF][START_REF] Zou | A nonintrusive proper generalized decomposition scheme with application in biomechanics[END_REF]. The external, home-made code, manages the assigned FE problems to the commercial software. The assigned problems are managed by classical PGD decomposition procedures [6,[START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]. Another way of constructing a computational vademecum through the use of non-intrusive PGD relies on the calculation of snapshots to later on construct a manifold of solutions using the locally linear embedding (LLE) technique for example or other manifold learning techniques [7,31,46,[START_REF] Giraldi | To be or not to be intrusive?the solution of parametric and stochastic equations -the "plain vanilla" galerkin case[END_REF][START_REF] Ibanez | A multi-dimensional data-driven sparse identification technique: the sparse proper generalized decomposition[END_REF][START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Tsiolakis | Nonintrusive proper generalised decomposition for parametrised incompressible flow problems in openfoam[END_REF].

Despite numerous developments during the recent years towards the development of computational vademecums in both intrusive and non-intrusive manners [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF], the calculation of solution handbooks using any of the aforementioned methods faces a major challenge when changing a parameter induces drastic changes in the solution map. The changes at a given point can be very aggressive such that interpolation, a finite element mesh in a parametric domain, for instance, leads to unrealistic numerical artifacts [START_REF] Ghnatios | Towards parametric rtm process: The interpolative mapping[END_REF]. Thus, in resin transfer molding, the interpolation of the resin flow front between solutions computed at different permeabilities may lead to multiple fronts appearing in the solution. For example, interpolating the solutions shown in We can clearly identify that a change in the permeability of the domain will compromise the solution of the resin front propagation and therefore computing a solution handbook becomes useless. Computational vademecums in resin transfer molding appear to be a cornerstone of interest when identifying anisotropic permeabilities [START_REF] Morren | 2d permeability tensor identification of fibrous reinforcements for rtm using an inverse method[END_REF] and even more importantly, for parametric simulation of resin injection in complex molds.

Other challenging problems with localized non-interpolable effects are path-dependent problems. Path dependent problems appear in different domains ranging from economy [START_REF] Dohi | A note on portfolio optimization with path-dependent utility[END_REF] to 3D printing, structural analysis [START_REF] Tsay | Nonlinear structural design sensitivity analysis for path dependent problems. part 1: General theory[END_REF][START_REF] Zhu | Experimental research on transition from scale 3d printing to full-size printing in construction[END_REF], time-dependent vehicle routing problems [START_REF] Huang | Time-dependent vehicle routing problem with path flexibility[END_REF] and welding [START_REF] Wang | Spot welding robot path planning using intelligent algorithm[END_REF] among many others. Recently, the use of experimental and machine learning techniques to optimize path dependent problems in welding were considered, for example [START_REF] Ojo | Loop travel-path of fibre laser welded alclad aa2219-o alloy[END_REF][START_REF] Wang | Spot welding robot path planning using intelligent algorithm[END_REF]. However, the problem of simulating and optimizing a path-dependent problem involves an infinite number of parameters and is impossible to simulate and optimize using classical techniques.

In this part, we aim to tackle problems with non-interpolable localized effects, through proposing an appropriate interpolation technique. The suggested method starts by using a "physics-based" geometrical mapping, to transform different offline computed snapshots into another geometrically similar domain. In the new domain, which will be named "the reference domain", all localized effects overlap, thus rendering the interpolation possible. Once the interpolation is performed, an inverse transformation is used to transfer the constructed solution back to the real domain. Eventually the transformation and interpolation are both performed online, as fast as possible, to allow real time exploration and optimization of the solution. Several examples are illustrated in this paper using those adopted mappings, in Resin Transfer Molding (RTM), path dependent heat transfer problems and sheet molding compound (SMC) involving flow fronts, starting with an illustrative 1D example.

An illustrative example in 1D

As a first example, we consider an elastic rope fixed at its two ends A and B, with an applied displacement at a point P i of coordinate x = X as illustrated in Fig. A.52. For the sake of simplicity, the vertical displacement is assumed to be equal to 1 unit. The generated deformation of the rope is assumed linear between the two fixed support and the position of application of the load.

v = -x X if x < X, v = -L-x
L-X elsewhere.

(A.78)

The solution is illustrated in the 2D domain (x, X) in Fig. A.55(a) for L = 1. With the help of a geometrical transformation, we can transform the point of application of the displacements P i to the same position, for example the center point of the domain. Thus, we can transform the domain x ∈ [0, L] into another domain s ∈ [0, 2] such that:

x = sX, if s ∈ [0, 1], x = (s -1)(L -X) + X, if s ∈ [1, 2]. (A.79)
Using the transformation given by Eq. (A.79), the displacement is always applied at s = 1, and therefore the solution in the (s, X) domain results the one illustrated in As one may notice, in the s domain, displacements are always applied at the same position, and therefore the solution can be interpolated without major issues. To check the described approach, the solution is computed for four different positions of the applied displacement in the x domain, named 4 snapshots, at X = 0 + , X = L/3, X = 2L/3 and X = L -.

Geometrical mapping in 2D

In this section we illustrate the geometrical mapping used to transform the domains. For the sake of simplicity, we are first illustrating the process on a 2D example of a localized heat source and homogeneous Dirichlet boundary conditions. Thermal fields are computed in the

-k ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 = Q, (A.80)
where k is the thermal conductivity and Q the heat generation term. Q is defined as a localized heat generation inside a domain G such as: ) can be computed using any direct numerical solver. For instance, in this work we used the finite element method.

Q = 1 if (x, y) ∈ G,
Once the snapshots are available in the physical (x, y) domain, we can transform the domain into a "reference" domain (r, s) where all heat sources apply at the same location, as illustrated in Equation (A.82) is defined per triangle in the domains illustrated in figure A.60, φ i being the Lagrange interpolation functions defined in the reference (r, s) domain, while ξ i the Lagrange shape functions defined in the (x, y) real domain [START_REF] Reddy | An introduction to the finite element method[END_REF]. The element connectivity is eventually conserved while mapping between the real and reference elements. The values X i and R i are the coordinates of the nodes of the corresponding triangles in the (x, y) and (r, s) domains respectively.

X = X 1 φ 1 (r, s) + X 2 φ 2 (r, s) + X 3 φ 3 (r, s) R = R 1 ξ 1 (x, y) + R 2 ξ 2 (x, y) + R 3 ξ 3 (x, y) (A.82)
In the (r, s) domain, the interpolation can be performed using any classical interpolation technique, for instance the SSL-PGD technique is used [START_REF] Ghnatios | Towards parametric rtm process: The interpolative mapping[END_REF].

Flow fronts interpolation

Parametric domain properties

The propagation of a front remains a tricky issue when using again interpolation. In this section we tackle the propagation of a front in resin transfer molding (RTM) processes. The front propagation can be computed with a non intrusive approach using any classical software. We compute first different snapshots using different permeability values. The technique introduced in Section 3.3 does not apply in the present case, being the injection point located at the same position and remaining unchanged. In that situation a possible solution consists of transforming the 2D domain (x, y) into the 2D polar coordinates domain (r, θ) attached to the injection point. In (r, θ), for the considered problem, the front of snapshot one is defined from R 1 (θ), while the front in snapshot two reads R 2 (θ). We can therefore interpolate the level sets R 1 and R 2 at each value of θ to define R i related to k i = 0.5k. The resulting solution is depicted in figure A.65. Figure A.66 shows the error maps in the domain with respect to the exact solution, for both cases using classical interpolation and the one here proposed, spurious free, interpolation technique. An excellent accuracy can be noticed. Only very few nodes on the outer interface of the propagation region are not exactly captured, while classical interpolation exhibits completely wrong results. This idea can be further generalized for different values of k x and k y .

Flow front prediction with the injection point location as parameters

In this section, the techniques discussed in Sections 3.3 and 3.4, can be combined for simulating the flow front evolution for any permeability and any position of the injection point using only few snapshots. To check the method, we generated a matlab GUI application to modify the permeability values and the injection point coordinates. The solution is therefore obtained in real time. Moreover, we can compute the associated pressure fields at different points and at different time steps. However a slight adaptation for interpolating the pressure is needed. In fact, we first extract the pressure pattern P (r, θ) at every θ j , with θ j ∈ [0, 2π]:

P 1 j = P (r 1 , θ = θ j ) for R = R 1 , P 2 j = P (r 2 , θ = θ j ) for R = R 2 . (A.83)
Later on, we scale the pressures P 1 j and P 2 j to act in the domain with the flow front located at R j :

   P 1 j scaled (r 1 • R j R 1 , θ j ) = P 1 j (r 1 , θ j ), P 2 j scaled (r 2 • R j R 2 , θ j ) = P 2 j (r 2 , θ j ).

Conclusion

In this work we illustrated the possibility of constructing computational vademecums by adapting the interpolations of snapshots with highly localized effects such as front propagation and localized deformations. This approach relies on the use of physically-based mappings. The approach is later leveraged for different applications allowing real-time simulations in previously unreachable parametric simulations. The approach does not require any re-meshing of the domain since all the mathematical procedures apply in the reference domain. we decompose the 7D problem into a sequence of 2D/1D problems which allows the curse of dimensionality to be circumvented [23,[START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF]. In the following sections, the model has been adjusted according to the experimental AFM nanoindentation results and used to identify the depth-and rate-dependent fluid permeability through the solid mesh, consisting of polymer brushes.

The obtained "bank of 3D solutions" or "computational vademecum" has subsequently been coupled to an identification algorithm, in order to identify the parameters of the poroelastic material, for example the permeability. The reader should note that having the 7D parametric solution or computational vademecum, we dispose explicitly of all the sensitivities of the solution with respect to the model parameters. The parameters' identification can therefore be performed within a few seconds on a standard portable PC using classical algorithms such as the Newton method or Levenberg-Marquart algorithm [26,[START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF][START_REF] Ghnatios | Strategies avancées d'optimisation et analyse inverse basées sur l'utilisation de la pgd[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF].

Modeling of AFM nanoindentation of soft matter

AFM nanoindentation of the poroelastic brushes is performed using a rigid spherical probe attached to a cantilever. The system is illustrated in Figure B.2. The indentation process is controlled by the motion of the base of the cantilever, z(t) as shown. An optical sensor measures the displacement of the cantilever right above the colloidal probe (the indenter).

For the rest of this work, the cantilever is modeled as a spring. The cantilever adopted in this work has a resulting stiffness k s = 1.69N/m. An increasing force is applied on the spring to indent the specimen.

z = 0 z(t) Specimen Indenter w(t) Figure B.2:
The reference for z = 0 is taken at the point of contact of the indenter with the specimen, w(t) is the penetration depth of the probe into the specimen Assuming a quasi-static approach, the equilibrium equation for the force acting on the indenter can be written as a sum of liquid and solid contributions to the reaction force from the indented material. Therefore the equilibrium equation can be written as:

F f luid + F solid = F spring (B.1)
F f luid being the fluid contribution to the reaction force overcoming the spring compression force, F solid being the solid contribution to the reaction and F spring is the force imposed by the spring on the specimen. Knowing k s , the spring force can be written as:

F spring = k s (z(t) -w(t)) (B.2)
where w(t) represents the indentation depth into the specimen. The reference of z(t) is considered as the position of the spring base at the initial contact of the indentation probe with the specimen.

Note that the force of the spring is time dependent during the penetration phase. In the following, for the sake of notational simplicity, the force in the spring is simply denoted F (t). Next, we model the individual contributions of the fluid and the solid in the biphasic material.

The fluid contribution to the reaction force

The force in the fluid F f luid can be found by integrating the pressure of the fluid over the area of the indenter. By simple projection onto the z axis of the pressure distribution, we can demonstrate that the fluid force on a spherical indenter is equal to that on a flat indenter with the same contact radius. With ρ as the radius of the indenter, and using the geometrical considerations shown in Figure B.3, the contact radius between the fluid and the indenter can be written as:

R(t) = ρ 2 -(ρ -w(t)) 2 (B.3) ρ ρ -w(t) w(t) R(t)

Figure B.3: The contact radius R(t) of the indenter with the specimen

As z(t) increases, w(t) also increases as well as the contact radius R(t) between the probe and the sample. In the model, one should take in consideration the changes in the geometry of the domain.

The fluid in the specimen is confined in a porous elastic solid. Therefore, the pressure can be computed using the Darcy equation:

v = - K µ ∇P (B.4)
v being the fluid velocity, µ the fluid viscosity and K the effective permeability, which takes into account the material porosity. Moreover, neglecting the density variation, the conservation of mass can be written by [15,[START_REF] Wu | Dynamic compression of highly compressible porous media with application to snow compaction[END_REF]:

∇ • v + ∂φ ∂t = 0 (B.5)
Where φ is the porosity of the domain. Equation (B.5) takes in consideration the trasient variation of the porosity through the term ∂φ ∂t . Moreover, φ is defined by [15,[START_REF] Wu | Dynamic compression of highly compressible porous media with application to snow compaction[END_REF]:

φ = 1 - H 0 H (1 -φ 0 ) (B.6)
H 0 being the initial height of the domain, φ 0 the initial porosity and H the actual height. In the tested specimen, the initial porosity φ 0 is measured via swelling to be 0.9614 [START_REF] Li | Surface-grafted, covalently cross-linked hydrogel brushes with tunable interfacial and bulk properties[END_REF]. In this study, the maximum indentation velocity does not exceed 80 µm/s and, for the given initial porosity, neglecting the porosity change generates a maximum error of 0.52% in the computed pressure field for the most critical situation. Therefore, the conservation of mass can be simplified to:

∇ • v = 0 (B.7)
Combining equations (B.4) and (B.7), assuming K and mu homogenous in the domain, we may find the pressure equation to solve:

K µ ∆P = 0 (B.8)
complemented with the boundary conditions specified later in Section 1.2. The permeability of the domain K, the radius R, the height of the domain H and the compression velocity U are time dependent, as shown in Figure B.4. Since we assume a quasi-static approach of the probe, we may consider at a generic time step n:

U n =
w nw n-1 ∆t (B.9)

H n = H n-1 -U n • ∆t (B. 10
)
where ∆t is the time step. In what follows, the initial height of the specimen H(t = 0) will be denoted by H 0 . The pressure field, given by Equation (B.8) will be solved again at each time step n, updating the values of the three parameters R, H and U, as well as identifying the permeability K from the experimental data. It is here that the power of the PGD comes into play, since it considers all four parameters to be extra coordinates of the problem, and therefore it can solve the problem for any value of these parameters as detailed in Section 1.2.

F (t) z(t) R(t) H(t) w(t) U(t)

Solid contribution to the reaction force

Two approaches were considered for modeling the reaction force in the solid phase of the biphasic poroelastic material. Depending on the comparison with the experimental results, one of the two approaches shall be considered as reliable.

First approach: A rods-based model

The first tentative model considers the solid domain to be a system of rods, as shown in The "rods" approach to modeling the solid part Using the mechanics of materials elementary theory, we found the solid contribution to the reaction force as presented in Appendix A. The identification of the modulus of elasticity of the solid rods E rods can be achieved in the steady-state phase, when the probe indentation depth is fixed at w steady . In fact, after a sufficiently long period of time, the velocity in the fluid phase, and consequently the pressure gradients, tend to zero. Considering the surrounding pressure to be zero, one may find the pressure in all the fluid phase to be zero and only the solid contribution to the reaction remains in the steady state phase. Thus we may estimate the modulus of elasticity of the rods, E rods .

Second approach: Rate dependent model

Another way to define the reaction of the solid component to the indentation force would be the reaction while indenting the poroelastic material with a very slow indentation velocity U. As U tends to zero, the velocity of the fluid components, and consequently the pressure gradients, can be considered as being negligible. In this case, and again considering the surrounding pressure to be zero, the only contribution to the reaction force would be the solid component. We may therefore define the solid contribution to the force as:

F solid = lim U →0 F spring (B.11)
To experimentally measure the force in the solid at any indentation depth w, we set the indentation velocity ż to the machine minimum velocity 0.5µm/s. The spring force obtained at this velocity is then considered as being the contribution of the solid to the reaction force, which is supported by the fact that no significant hysteresis between the indentation and retraction force of the colloidal probe was observed. Both approaches are compared to experimental results in Section 1.4 before identifying the permeability values.

PGD solution of the multidimensional Darcy equation

The problem to be solved in this section is given by the pressure equation (B.8) in a cylindrical domain of height H and radius R with the following boundary conditions:

     v = 0 at z = 0 v = (0, 0, -U) at z = H P = 0 at r = R (B.12)
which consists of a zero velocity on the bottom support, the indentation velocity U at the top of the domain and zero pressure at r = R.

The resulting problem can be solved using a steady-state simulation, however the same problem shall be solved again (more than once) for each time step, at each change in the parameters U or the permeability K (to be identified), and at each change in the dimensions of the domain R and H. At this stage, using the PGD to solve the problem using all geometrical and configurational parameters as extra coordinates of the problem appears to be a suitable approach [START_REF] Ghnatios | Poroelastic properties identification through mico indentation modeled by using the proper generalized decomposition[END_REF].

Introducing the parameters as extra-coordinates of the problem

The geometrical parameters do not appear explicitly in the differential equation and can be introduced using geometrical transformations [START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF][START_REF] Ghnatios | On the space separated representation when addressing the solution of pde in complex domains[END_REF]. Therefore, the following mapping is defined:

     x = R • ξ y = R • γ z = H • φ (B.13)
The domain in which we solve the problem is (ξ, γ, φ) such as (ξ, γ) is a 2D circle of radius one and center O(0, 0) and φ = [0, 1], a 1D domain. The geometrical transformation defined in (B.13) transforms the reference (ξ, γ, φ) domain into a cylindrical (x, y, z) domain of radius R and thickness H.

The Jacobian matrix J of the transformation defined in (B.13) can be written as:

J =    R 0 0 0 R 0 0 0 H    (B.14)
Therefore the determinant of J is given by the product of the diagonal, which can be written as R 2 • H. Substituting the transformation (B.13) into the governing equation (B.8), we can find the final equation to solve in the (ξ, γ, φ) domain:

ξ γ φ K µ H ∂ 2 P ∂ξ 2 + H ∂ 2 P ∂γ 2 + R 2 H ∂ 2 P ∂φ 2 dξdγdφ = 0 (B.15)
Equation (B.15) constitutes the final equation to solve using the PGD. The boundary conditions of the problem defined in (B.12) shall now be projected in the (ξ, γ, φ) domain. Thus the new boundary conditions in the solution domain can be written:

       v = 0 => ∇P = 0 at φ = 0 v = (0, 0, K µ ∇P • N) = (0, 0, -R 2 • U) at φ = 1 P = 0 elsewhere (B.16)
N being the outward normal to the cylindrical domain (ξ, γ, φ) on the upper surface defined by φ = 1. It should be noticed that the permeability K appears naturally in the final

Simulation results

A triangular mesh is used in the in-plane (ξ, γ) circular domain composed of 363 nodes, while the thickness mesh along φ contains 1000 nodes. The mesh along the parameters consists of 10,000 nodes for each coordinate, logarithmically distributed. The obtained solution is defined in a high-dimensionality space, therefore advanced visualization tools need to be used [6]. In Moreover, one may wish to illustrate the variation of the pressure at a given point of the domain as a function of the parameters. This is possible using the PGD solution. In fact we can specify the coordinates of the point in the (x, y, z) domain and illustrate the variation of the pressure in the parametric representation, as shown in Figure B.12, where the variation of the pressure field is illustrated as a function of R, H and K at a given point (x,y,z) and for a given compression velocity U.

Verification of the solution

To verify the obtained results, we compared the PGD solution for a combination of the 4 parameters to a finite-element (FE) solution computed for the same values of the parameters H, R, U and K. The cross-sections at z = H/2 for the PGD and FE are shown in Figures B. 13 The maximum relative error does not exceed 0.05%. This small relative error is generated only by the convergence of the PGD, since when increasing the number of terms involved in the PGD separated representation, the PGD solution tends to the finite elements one [45], however the solution here considered only involves 200 terms, i.e. o = 200. In any case, the maximum relative error does not exceed 0.05%.

Experimental study

In this section we describe the experimental study performed to measure the forces on the indenter F at each time step n as a function of the indentation depth w.

Materials

p-(chloromethyl)phenyltrimethoxysilane (ABCR, Germany), acrylamide ( >99.5%, Fluka-Chemie AG, Switzerland), triethylamine (>99.5%, Sigma -Aldrich, Switzerland), sulfuric acid (95-97%, Sigma-Aldrich, Germany) and hydrogen peroxide (30 wt% in water, Merck, Germany) were all used as received. Tetrahydrofuran (99.5% extra dry, Acros, Germany) and toluene (>99.7%, Fluka-Chemie AG, Switzerland) were freshly distilled over sodium prior to use. Sodium N,N-diethyldithiocarbamate (97%, Fluka, Switzerland) was recrystallized from methanol. Water was deionized with a GenPure filtration system (18.2 MΩ cm, TKA, Switzerland). HEPES buffer was prepared by dissolving 1 mM 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES, BioChemika Ultra, Fluka, Switzerland) in Milli-Q water, and adjusting the solution pH to 7.4.

Brush synthesis

Surface-grafted PAAm brush films were fabricated by means of a "grafting-from" strategy on an iniferter-modified silicon substrate via UVLED surface-initiated polymerization (UVLED-SIP) described previously in [START_REF] Li | Surface-grafted, covalently cross-linked hydrogel brushes with tunable interfacial and bulk properties[END_REF][START_REF] Li | Poly(acrylamide) films at the solvent-induced glass transition: adhesion, tribology, and the influence of crosslinking[END_REF]. The reaction was initiated from dithiocarbamateterminated, silane-based self-assembled monolayers (SBDC) on piranha-cleaned silicon wafers (P/B <100>, Si-Mat Silicon Wafers, Germany). Prior to the polymerization, an aqueous acrylamide (AAm) solution (1 M) and the flask containing SBDC-modified wafer were degassed/flushed with Ar for 15 minutes. The degassed monomer solution was transferred to the substrate-containing flask under an argon atmosphere. The polymerization was started on the SBDC-modified substrate in an aqueous environment under the irradiation of the UVLED setup, which was built in-house and incorporated a LED (NCSU033A, NICHIA Corp., Japan) with a narrow emission spectrum at 365 ± 5 nm. The UV irradiation intensity (5.5 mW/cm 2 ) was measured at the substrate position using a radiometer (UVX radiometer with UVX-36 sensor, UVP, Upland, CA, USA). Substrates were exposed to UV light for 30 min. After the polymerization, the substrates with PAAm thin films were rinsed with Milli-Q water, then immersed in water for 24 hours, in order to remove any entangled or weakly bonded polymers.

The dry thickness of the surface-grafted PAAm films was determined with a spectroscopic ellipsometer (VASE, M-2000F, LOT Oriel GmbH, Darmstadt, Germany) at an incident angle of 70 o . The dry polymer was assumed to have a refractive index of 1.45. The thickness was determined via the analysis of a four-layer (Si/SiO2/PAAm/ambient) model with known thicknesses and refractive indices of the Si, SiO2 and ambient (software WVASE32, LOT Oriel GmbH, Darmstadt, Germany). The obtained dry thickness of the used polymer film was thus 122 nm ± 2 nm. The swelling ratio in HEPES aqueous buffer solution was measured to be 25.9 ± 1 [START_REF] Li | Poly(acrylamide) films at the solvent-induced glass transition: adhesion, tribology, and the influence of crosslinking[END_REF]. The swollen thickness used for the modelling was therefore assumed to be 3160 nm ± 170 nm. the Au-coated tipless cantilever (NSC-36, Mikromash, Estonia) was measured by the thermalnoise method before attaching the colloidal microsphere [39]. A silica microsphere of 12 µm in radius (EKA chemicals AB, Kromasil) was glued with UV-curable glue (Norland optical adhesive 63) to the end of the tipless cantilever by means of a home-built micromanipulator [START_REF] Ducker | Measurement of forces in liquids using a force microscope[END_REF][START_REF] Ducker | Direct measurement of colloidal forces using an atomic force microscope[END_REF]. Approach rates were set to 0.5, 2.0, 8.0, 40 and 80 µm/s. The exact approach rates ż(t) were measured using a Z-piezo sensor. True indentation depth w(t) as well as indentation rate ẇ(t) were then calculated by subtracting the deflection of the cantilever at the tip from the measured Z-sensor displacement. An acquisition rate of 2000 Hz was used, resulting in measured data points in time steps of ∆t = 0.5ms. The set-point for the indentation force was 18 nN, which corresponds to 500-700 nm of indentation depth. Several measurements were performed at each rate at at least four different locations of each sample. Representative curves used in the model are presented in Figure B. 16. The methodology used to define the exact point of contact at each rate of indentation is detailed in [START_REF] Mathis | Investigation of Fluid-Confinement in Lubricating Polymer-Brush Coatings[END_REF]. The indentation results shown in Figure B.16 will be used in Section 1.4 for comparison with the simulation results and to identify the permeability of the indented specimen.

Colloidal-probe AFM nanoindentation

Identification of permeability and comparison of the two adopted models Solving the inverse problem

In this section, we compare the experimental results with the simulated models. Before starting any identification process, we fitted F (w) and U(w) obtained from the experimental data to avoid numerical oscillations and divergence of the inverse problem. We illustrate some examples in Figure B.17 For the solid force computation, in the first approach we use the rod model previously described in Section 1.1. However, we should start by identifying the modulus of elasticity of the solid part. One way to identify the solid contribution to the reaction force is through indenting to a constant depth then allowing the system to reach a steady state at this indentation depth. Thus, the fluid velocity and consequently the pressure gradients, will tend to zero. While minimizing Equation (B.17) we overestimate the solid contribution, which becomes larger than the total spring force at some point. This indicates that very probably the rods are subjected to buckling and consequently the solid effective modulus cannot be computed by using such an approach.

A more accurate approach could therefore use the total reaction or spring force at a very low velocity ż ≈ 0 as the contribution of the solid part to the total reaction force. Using this approach, we illustrate the contribution of the solid force, the fluid force and the fitted spring force in Figure B.19 and B.20 for a penetration velocity ż = 8µm/s and ż = 80µm/s respectively. In Figure B.19, the blue curve labeled spring force is the experimentaly measured force in the spring at ż = 8µm/s, the solid contribution to the reaction force in green is the measured spring force at ż ≈ 0µm/s, and the fluid contribution to the reaction force in red is the difference between them. The fluid contribution to the reaction force is also the integral of the pressure over the contact area between the indenter and the tested specimen. One may notice that the fluid reaction force exceeds the solid one at high indentation rates, ż = 80µm/s for instance, as shown in Since deriving a function numerically at a node n requires knowing its value at nodes n+1 and n-1 (using Euler's centered derivation scheme), and since we also require the function's value at node n, in order to use Newton's algorithm, one iteration therfore requires therefore identifying the value of the function at 3 nodes. Thus, to obtain the fluid contribution, one has to solve Darcy's problem defined by Equation (B.8) a least 3 times for each iteration of Newton's algorithm, updating all the parameters and the value of the permeability. How- .20: Contribution of the solid and fluid force to the total force while displacing the cantilever base at a velocity of ż = 80µm/s, the fluid reaction force exceeds that of the solid ever, using the PGD, we have already computed the solution for all possible values of the parameters and we access the exact solution and sensitivity with respect to the permeability K in real-time. Therefore identifying the permeability (using the second approach to define the solid contribution to the reaction force) over 203 time steps using a Newton algorithm at every time step took only 34.38 seconds on a standard laptop. The identified permeability for different velocities ż is shown in Figure B.21. Due to the difficulties in determining the viscoelastic properties of the polymer brush in a swollen state it was not possible to distinguish the rate-dependent stiffening effects of the solid from the pure fluid contribution to the total stiffening. Therefore, the viscoelastic contribution of the polymer had to be lumped in the calculated permeability, which is addressed as the effective permeability. , the permeability as a function of indentation depth at each indentation rate shows an initial increase followed by a decreasing trend at larger indentation depths, i.e. larger strains. This type of behavior is observed for all indentation rates with permeability values increasing with increasing indentation rate. Moreover, the peak permeability for each indentation rate moves towards larger indentation depths when the indentation rate is increased.

Discussion

The obtained effective permeability, which includes pure fluid contribution as well as possible viscoelastic stiffening of the polymer brush, is plotted as a function of indentation depth in Figure B.21 and shows three different phenomena. One is the initial increase of permeability with indentation depth. This initial increase could be explained by a conformation of the brushes and their deformation in the lateral direction -the direction of the flow. Such lateral deformation could ease the flow through the mesh, resulting in increased permeability. Using a material with substantially lower Poisson's ratio should result in limited lateral deformation during compression and thus limited or even absent initial increase of permeability [START_REF] Mansour | The permeability of articular cartilage under compressive strain and at high pressures[END_REF].

The second obvious phenomenon in Figure B.21 is the monotonic decrease in permeability with increasing compression after reaching a peak permeability. Comparable behavior was also observed in confined compression of cartilage and polymer brushes [START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF][START_REF] Mansour | The permeability of articular cartilage under compressive strain and at high pressures[END_REF][START_REF] Mathis | Investigation of Fluid-Confinement in Lubricating Polymer-Brush Coatings[END_REF][START_REF] Mow | Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiments[END_REF]. In these cases, the authors showed that permeability decreases exponentially with compressive strain. As shown in a study by Abbott et al. the monomer density is lower far from the surface and increases with decreasing distance from the surface [2]. During indentation the increase in monomer density with increasing indentation depth is amplified by the compression of the structure. This substantially decreases the pore size (i.e. porosity), which consequently causes the decrease in permeability.

The third phenomenon seen in Figure B.21 is an increase in values the of the permeability at higher indentation rates. Apparently, the higher velocity of the interstitial fluid, which is caused by higher indentation rate, is accompanied by a relatively small increase in interstitial pressure. This is clearly shown in the Appendix B in Figures ?? and ??, which show the magnitude of the velocity v = -K µ • ∇P at point A(0, R/2, H/2) and the pressure at the same point as a function of the indentation depth w. We also show the amplitude of the pressure's gradient ∇P as a function of w in figure ??. We may notice that the pressure's gradient increases with increasing ż (8 times such as ∇P ż=80µm/s ≈ 8 × ∇P ż=2µm/s , the comparison is performed at the last timestep when indenting at ż = 80µm/s, at w = 0.514µm), while the velocity magnitude increased about 33 times with the same increase of ż. In Figure ?? we note a decreasing slope of the pressure's gradient as a function of w. However, as shown in Figure ??, the slope of the velocity v decreases more significantly as a function of w, which reflects the decrease of the permeability. We can also plot the velocity magnitude as a function of the pressure at point A, as shown in Figure B.22. We notice a large slope of v = f (P ) for low pressure due to increased permeability, and a stagnation at higher pressures as a result of decreased permeability. This behavior is different from what was observed for bovine cartilage, where permeability decreases substantially with increasing pressure over the entire range of compression [START_REF] Mansour | The permeability of articular cartilage under compressive strain and at high pressures[END_REF][START_REF] Mow | Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiments[END_REF]. The difference in behavior could be due to the complex structure of cartilage, where brushy charged proteoglycans are trapped in the stiffer collagen mesh, thus preventing the flow, reducing permeability and giving rise to interstitial pressure from the very beginning of compression. The difference in behavior could also be explained by relatively long reptation times of the polymer brush, due to its length. Increased indentation rates decrease the time available for rearrangement of the polymer.

The fact that the depth at which the maximum permeability is reached increases with indentation rate, could be due to competition between pore expansion at increased fluid pressure (at higher rates) and pore compression due to increased compressive strain. However, to prove this assumption, further experiments are needed.

Overall the poroelastic behavior observed is in line with characteristic properties of various biphasic materials, such as cartilage and hydrogels [38,[START_REF] Moore | Tribological rehydration of cartilage and its potential role in preserving joint health[END_REF][START_REF] Mow | Fluid transport and mechanical properties of articular cartilage: a review[END_REF]. Further studies, both experimental and via modeling, are required to extend the current understanding of fluid confinement and its role in defining rate-dependent materials properties.

Conclusion

In this work, we tackled the issue of identifying the properties of soft and thin poroelastic materials in a liquid environment. To efficiently identify the permeability, a high-dimensionality Darcy problem defined in a 7D space has been solved using the PGD. The potential of using PGD as a simulation tool to identify challenging materials parameters of soft biphasic matter is presented. A comparison with the finite-element approach confirmed the PGD results and revealed the potential of this technique, which allows the curse of dimensionality to be circumvented by solving the 7D problem as a sequence of problems in lower-dimensional spaces. Subsequently, calculation of permeability proceeds from this parametric solution by means of an inverse technique (Newton method), in order to fit the experimental data extracted from indentation tests performed at different indentation rates.

The calculated permeability seems to depend on the indentation rate as well as on the compressive strain. We attribute the initial increase in permeability with increasing strain to the conformation of the polymer brush, which eases the liquid flow through the mesh. The decaying permeability towards higher strains is most probably a consequence of the decreasing mesh size due to compressive deformation. The fact that permeability increases with deformation rate suggests that for PAAm brushes the flow velocity of the interstitial fluid increases more than the interstitial pressure when higher indentation rates are applied.

The combined approach using novel experimental and modeling techniques has enabled new insight into the poroviscoelastic properties of hydrated soft materials. The obtained results are encouraging, yet also show current concerns and limitations to be addressed to fully comprehend the physical and chemical interplay in both biological as well as synthetic soft matter.

Simulation of 3D printing 2.1 Introduction

In the past decade, additive manufacturing attracted interest because of its benefits and flexibility. It made possible the production of designs that were impossible to manufacture using traditional manufacturing techniques and opened the door to the advancement in topology optimization and design efficiency enhancement. Additive manufacturing is a wide technology with many variants. Among them, Selective Laser Melting -SLM-also known as Direct Metal Laser Sintering, that proceeds on a powder bed [22].

SLM consists of a concentrated laser that hits a metal powder layer and melts it according to the slicing of a CAD model. The powder is heated above the liquidus temperature of the metal hence producing parts with properties comparable to those of metal bulk, except in what concerns ductility.

However, due to the extreme process conditions some defects occur, such as cracking [START_REF] Gong | Analysis of defect generation in ti-6al-4v parts made using powder bed fusion additive manufacturing processes[END_REF] and part distortion [START_REF] Shiomi | Residual stress within metallic model made by selective laser melting process[END_REF] induced by installed residual stresses [START_REF] Heigel | In situ monitoring and characterization of distortion during laser cladding of inconel® 625[END_REF]. In order to account for these defects, scientists have tried to quantify them using either experimental methods as well as numerical simulation.

Experimentally, two types of measurements are used: with and without contact. In absence of contact, ultrasounds, thermal imaging such as infrared -IR-sensors [START_REF] Clijsters | In situ quality control of the selective laser melting process using a high-speed, real-time melt pool monitoring system[END_REF] or X-ray tomography to quantify the part defects such as porosity and cracks [START_REF] Du Plessis | Effects of defects on mechanical properties in metal additive manufacturing: A review focusing on x-ray tomography insights[END_REF], are widely considered, among many other techniques. This type of techniques can only inspect external exposed surfaces or its neighborhood inside the material depending on their penetration capability. Contact-based techniques use thermocouple joints to measure the temperature of the substrate (bottom surface or bed) [START_REF] Denlinger | Thermal modeling of inconel 718 processed with powder bed fusion and experimental validation using in situ measurements[END_REF] or other sensors enabling accessing to the installed stresses.

To circumvent the shortcomings of the experimental approach, researchers resorted to numerical techniques in an attempt to simulate the printing process. Among the numerical techniques, the finite elements analysis -FEM-is mostly used to simulate the temperature field in additive manufacturing processes [START_REF] Loh | Numerical investigation and an effective modelling on the selective laser melting (slm) process with aluminium alloy 6061[END_REF] [START_REF] Li | Thermal behavior during selective laser melting of commercially pure titanium powder: Numerical simulation and experimental study[END_REF]. A thermal simulation is then coupled to a mechanical simulation in order to predict the distortion of 3D printed parts. Simulation concerns three different scales: Micro, Meso and Macro.

Microscale analysis focuses on the rich local physics including the interaction between the heat source and the substrate, the phase-change and the associated metallurgical transformations, as well as the localized thermal gradients. The fast cooling and the induced stresses along the printing path is usually described at the mesoscopic scale. These two approaches are usually coupled [58] [181] for reproducing quite well experimental findings. However, the computational cost makes impossible envisaging its application at the part level (macroscopic scale).

In order to validate the numerical results, several approaches are used depending on the considered process. For SLM, the so-called "Twin Cantilever" test shown in Fig. B.23 is usually considered. In this experiment a T-shaped part is printed on a support and the tip deflection after removing the support serves to calibrate and/or validate.

However in most additive manufacturing processes, distortion is difficult to predict because of the change of the mechanical properties during the process [START_REF] Sun | The influence of as-built surface conditions on mechanical properties of ti-6al-4v additively manufactured by selective electron beam melting[END_REF]. The high temperatures cause the microstructure to be altered, and the cyclic heating may cause for some materials a thermal fatigue leading to the material embrittlement [21]. These changes are difficult to model and take into account since they depend on the printing parameters such as the laser power and speed, the material itself, the powder characteristics, the heating spot size, ... affecting not only the microstructure but also the interlayer cohesion due to the varying porosity [START_REF] Du Plessis | Effects of defects on mechanical properties in metal additive manufacturing: A review focusing on x-ray tomography insights[END_REF].

In these validation experiments, the printing parameters are varied and the distortion is studied. A quite good agreement between the experimental and numerical results are noticed when addressing the thermal modeling. Whereas in what concerns the distortion the error varies depending on the process parameters and materials employed, ranging from 10% to a total mismatch [START_REF] Li | Efficient predictive model of cantilever distortion in selective laser melting[END_REF] [START_REF] Erik | Residual stress and distortion modeling of electron beam direct manufacturing ti-6al-4v[END_REF].

The present work aims at proposing a thermal simulation of SLM while studying the reheating effect when the heat source re-passes close to a solidified element. In addition, it will present a novel approach for simulating the distortion of a 3D printed part. This approach is used to generate a calibration model from which the stresses will be extracted and used to perform a macroscale simulation that gives the distortion as a function of the temperature reached by the molten metal and of the printing trajectories appropriately parametrized, in real time by using the non-intrusive Proper Generalized Decomposition approach [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] [START_REF] Quesada | Computational vademecums for real-time simulation of surgical cutting in haptic environments[END_REF].

Transient Simulation of a Moving Heat Source

This section discusses the finite element thermal modeling associated with a moving heat source hitting a powder bed and moving at a velocity V, as encountered in SLM.

Heat Equation

The heat transfer mechanisms in SLM are conduction and natural convection as shown in Fig. B.24. Note that radiation can be neglected since metallic materials have an excellent heat conduction ability, hence most of the heat will be diffused by conduction to the surrounding material [START_REF] Roberts | A three-dimensional finite element analysis of the temperature field during laser melting of metal powders in additive layer manufacturing[END_REF]. In the reduced order modeling that we develops in the present work, the lasermatter interaction is not taken into account, only the effective heat transferred to the material is considered.

The transient heat equation reads

ρCp ∂T ∂t = k ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 + ∂ 2 T ∂z 2 + q latent , (B.18)
where q latent is the latent heat, i.e. the heat released by the phase change.

The heat source is considered having a gaussian distribution [START_REF] Roberts | A three-dimensional finite element analysis of the temperature field during laser melting of metal powders in additive layer manufacturing[END_REF] and hence q source is expressed by

q(r) = 2P πr 2 0 e -2r 2 r 2 0 (B.19)
in which P is the laser power (W ), r 0 is the laser spot radius and r is the distance from the spot center. The initial temperature is given by T (x, t = 0) = 293K and the free convection boundary condition reads q = h(T -T ∞ ), with h the exchange coefficient. In the heat equation (B.18), Cp refers to the heat capacity and k to the thermal conductivity.

The implicit semi-discrete formulation reads

T n+1 -T n ∆t = k ρCp ∂ 2 T n+1 ∂x 2 + ∂ 2 T n+1 ∂y 2 + ∂ 2 T n+1 ∂z 2 + q ′ , (B.20)
in which T n+1 is the temperature field at time step n+1, i.e. t = (n+1)∆t, with T 0 = 293K, and q ′ is the latent heat to heat capacity ratio assumed temperature dependent [START_REF] Radovic | Numerical simulation of steel ingot solidification process[END_REF]

q ′ =      q ′ = 0 if T > T l q ′ = ∆H Cp dfs dT ∂T ∂t if T s < T < T l q ′ = 0 if T < T s , (B.21)
where T l and T s are respectively the liquidus and the solidus temperatures, ∆H is the latent heat and f s the solid fraction. The value of ∆H is taken from the literature [START_REF] Radovic | Numerical simulation of steel ingot solidification process[END_REF] assuming that the solidification in additive manufacturing is similar to that in casting. Obviously this hypothesis can be relaxed and other expressions considered. The previous equation shows that there is no heat release when cooling in the liquid state or solid state, but only when there is a solid-liquid mixture. In Eq. (B.21), the term dfs dT reads [START_REF] Pinheiro | Mould heat transfer with mould flux lubrication[END_REF] df

s dT = - 1 (T l -T s )(1 -2/π) 1 -sin π(T -T l ) 2(T s -T l ) . (B.22)
In order to simulate the moving heat source, at each time step ∆t the laser moves a distance ∆x = V ∆t. Hence the FEM discrete problem needs at each time step n, the solution of the linear problem expressed from

[K] n [T ] n+1 = [Q] n,n+1 , (B.23)
where both the right hand member and the coefficient matrix change at each time step, the former because the thermal source displacement and the last moreover due to the nonlinearities.

Thermal Properties

The thermal conductivity is a temperature dependent function and hence it must adapt to the temperature change at each iteration. The used material in the simulation is the 316L stainless steel. In [3] authors found correlations for the thermal conductivity and the temperature

k = a T b e cT e d/T (B. 24
)
where T is the temperature in (K) and a, b, c and d are constants corresponding to iron, reported in Table B.1. The thermal conductivity is assumed to be the one of the bulk material expressed by Eq. (B.24) everywhere except in the powder layer in which it is reduced of one order of magnitude for taking into account the effects of porosity in the powder bed. The heat capacity can be expressed from the mixture rule

C p =      C p,solid if T < T s f s C p,solid + (1 -f s ) C p,liquid if T s < T < T l C p,liquid if T > T l , (B.25)
where f s is the solid fraction expressed by [START_REF] Pinheiro | Mould heat transfer with mould flux lubrication[END_REF]

f s = (T l -T ) + 2 π (T s -T l ) 1 -cos π(T -T l ) 2(Ts-T l ) (T -T s )(1 -2/π) . (B.26)

Results and Discussion

Numerical simulations performed and discussed in what follows consider the parameters reported in Table B.1. Using these values, the simulation is performed in a thin 3D plate of length 15 mm, width 5 mm and thickness 1 mm, where a temperature snapshot is depicted in

Heating

When studying the maximum temperature at each point along the laser path, it results the profile depicted in Fig. B.26. This result clearly proves that except at the beginning, all the points on the laser trajectory will reach almost the same maximum temperature, for a given heat source that advances with a given velocity.

Then, the cooling profile of every point in this region of uniform maximum temperature is depicted in is of about 44 ms, time that is referred as critical time t cr . Thus, if the time elapsed before coming back to the vicinity of a certain point is higher than that critical time t cr , this latter has cooled down to the almost ambient temperature T amb and it will experience the same reheating as any point in the neighbor path-line, as the marked point in Thus, if the elapsed time between two passes is higher than t cr , the previous analysis proves that all the points start their reheating from an almost ambient temperature and thus will reach almost the same maximum temperature when reheated, and will follow a similar cooling.

However, when the elapsed time between two passes is lower than t cr , reheating applies on a still hot initial state. This situation applies when a printing path starts close to the end of the previous one as illustrated The previous analysis allows neglecting the thermal heterogeneities induced by elapsed times lower than t cr as soon as the printed part are large enough.

Thermo-Mechanical Modelling Parametric Distortion

In what follows we assume a thermal load ∆T = T max -T amb , with T max discussed in the previous section, and the deposition trajectories parametrized as proposed in [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF] where the trajectories coincide with the curves Φ = cte, with Φ(x, y) solution of the problem

∆Φ(x, y) = µ 5 (B.28) with (x, y) ∈ ω = [0, L] × [0, H],
and where the printed part Ω ⊂ ω.

Eq. (B.28) is solved with the boundary condition expressed from

         Φ(x, y = 0) = µ 1 + (µ 2 -µ 1 ) x L Φ(x, y = H) = µ 4 + (µ 3 -µ 4 ) x L Φ(x = 0, y) = µ 1 + (µ 4 -µ 1 ) y H Φ(x = L, y) = µ 2 + (µ 3 -µ 2 ) y H . (B.29)
The parametric solution performed by using the so-called Proper Generalized Decomposition -PGD-allows computing Φ(x, y, µ 1 , . . . , µ 5 ) that as discussed in [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF] allows describing most usual printing trajectories.

In [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF] authors assume that the local residual tension follows the printing trajectory, the only non-zero component of the residual stress tensor, σ 0 | (x 0 ,y 0 ) , as shown in Fig. B.34 where x 0 and y 0 represents local frame while x and y represents the global one.

The virtual work principle reads

Ω ǫ * : (σ -σ 0 )dΩ = 0, (B.30)
or by using the elastic behavior

Ω ǫ * : K : ǫ dΩ = Ω ǫ * : σ 0 dΩ, (B.31)
with K the fourth order elasticity tensor, and the residual tension scaling with ∆T as discussed in [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF].

Thermo-Mechanical Coupled Model

In this section we will present a new macroscale approach to simulate the distortion of a part processed by additive manufacturing. In this approach the cooling of the deposited metal will be divided into two stages: T > T melting and T < T melting . In the first stage, the metal is in a liquid phase and no residual stresses are induced because the liquid metal adapts to the substrate deformation that will expand when heated from the ambient temperature. Then, in the second stage, the metal solidifies and hence sticks to the substrate surface. The whole system (deposited later and substrate) continues the cooling from the temperature existing when the deposited layer started its solidification, until reaching the final ambient temperature. The flowchart is presented in Fig. B. 35. Predictions were validated experimentally as described in the next section.

Experimental Validation

The performed experiment consists of evaluating the distortion of a cantilevered plate with one deposit line on its top surface. The geometry is shown in Fig

. B.36.
The experiment is performed using Laser Metal Deposition -LMD-, with the deposit dimensions being: width 0.8 mm, thickness 175 µm and the substrate thickness 2 mm with its left boundary clamped (all the other remaining free). In addition the material used is stainless steel 316L.

For the thermal simulation, the transient heat transfer equation is solved with a prescribed initial temperature in the deposit layer T (t = 0) = 2000K and T (t = 0) = 293K elsewhere. As discussed in the previous section, the thermal simulation provides the temperature field in the substrate when the deposit-substrate interface reaches the melting temperature T s . The substrate heating produces the substrate plate deflection depicted in The quantity of interest is the maximum deflection at the tip of the plate. When comparing the simulated results with the experimental measurements of the maximum tip deflection, an excellent agreement was observed. Note that the modulus of elasticity of the deposit is taken to be slightly lower than that of the substrate due to the porosities and phase change during solidification (E deposit = 175 GPa and E substrate = 200 GPa) [START_REF] Riede | Material characterization of aisi 316l flexure pivot bearings fabricated by additive manufacturing[END_REF]

Calibration Model for the Macroscale Simulation

In order to perform fast simulations on large scale models, a high-resolution calibration model is proposed from which the residual stresses to be used as local stresses in the whole part (to proceed as discussed in Section 2.3) are extracted.

For addressing the solution of the associated thermomechanical problems in the extremely degenerated domains in which thicknesses are orders of magnitude much smaller that the other in-plane characteristic dimensions, the in-plane-out-of-plane separated representation allowing using thousands of elements along the thickness is retained.

That separated representation expresses the displacement field according to u(x, y, z)

≈ N i=1 P i (x, y) • T i (z), (B.32)
where • refers to the Hadamard product and P i and T i are vectors involving respectively the in-plane and out-of-plane coordinates.

The calibration model geometry consists of three deposited layers with three printed lines per layer as shown in Fig. B.40, all of them along direction y. In this model we assume that the deposit of interest, highlighted in red in Fig. B.40 of 0.8 mm width and 200 µm thickness, is a representative one accounting for the effects of the neighboring deposits and hence its stress state can be used at any point in a simulated part.

The simulation is performed by adding a deposit after the other to accumulate the distortions and the installed stresses. In the deformed model shown in Fig. B. 41 we can see that the highest distortion is observed between the tracks and between the layers and mainly on the edges. In addition, in the middle of the representative domain the distortions reached an almost uniform value, meaning that the deposit length is sufficient to extract the mechanical state.

From the computed strains we can then calculate the stresses and also the Von Mises effective stress. Fig. B.42 proves that there is a stress concentration on the interface between the layers which is a problem widely encountered in additive manufacturing [32]. The plot shows a region in the middle of the domain in which stresses are almost uniform and all in-plane stresses vanish. These stresses will be used as local stresses to enrich the parametric model proposed in [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF] whose assumption (residual tension aligned with the printing path-line) seems too simplistic when compared with the stresses reported in equivalent mesh resolution. It is important to note that few minutes in a standard laptop allowed for a resolution equivalent to 2.5 millions of 3D finite elements.

Parametric Macroscopic Modeling of Part Distortion

In this section we will present the application of the proposed methodology for the simulation of the distortion of a barrel. The simulated part geometry is shown in Fig. B.46. The part is made of Titanium Grade 5 and printed on a Titanium substrate.

The cross sectional area being circular we calculate the solution of the Poisson equation (B.28) in a slightly bigger rectangular domain. Then, the solution will be projected from the rectangular embedding domain onto the part.

In order to obtain circular trajectories, Eq. (B.28) is solved with µ 5 = 1 enforcing Φ = 0 outside the part domain, by using a penalty formulation. The technique proposed in [START_REF] Quaranta | Parametric evaluation of part distortion in additive manufacturing processes[END_REF] is extended by incorporating the richer stress state that is parametrized by the deposition temperature T . Thus, a parametric solution representing the part distortion u(µ 1 , µ 2 , µ 3 , µ 4 , µ 5 , T ) is computed. The deposition temperature can be easily correlated with the process parameters: deposition velocity and laser power. For computing such a parametric solution several runs (about one thousand) for different values of the parameters were performed and the computed results combined for defining the solution for any other choice by interpolation [START_REF] Ghnatios | Towards parametric rtm processes: The interpolative mapping[END_REF].

Moreover, to avoid axial distortion occurring when all trajectories start at the same location, the starting positions will be randomly chosen in each layer. Using this parametric solution one can modify the trajectory or process parameters and evaluate their impact in almost real-time.

The methodology performs in larger parts without major difficulties, as the one illustrated in

Conclusion

This paper started by discussing transient thermal simulations of SLM. The reheating effect when the heat source re-passes next to an already printed point of a previous pass was quantified. The analysis proved that if the metal had enough time to cool down to reach the ambient temperature, the heat source will reheat it to an almost representative temperature ensuring a certain uniformity of the mechanical state all along the part.

Then a thermo-mechanical coupled approach for calculating the distortion in a 3D printed part was presented. This model was used to propose a calibration procedure in order to extract the characteristic stresses to be used for calculating the distortion at the part level. Finally, several simulations were done by varying the trajectory and the temperature of the deposited material to define the parametric part distortion with respect to the process parameters.

159 solution iteratively starting from an initial one built by some heuristic until a stopping criterion is met. The stopping criterion can be elapsed time, number of iterations, etc. The operation of meta-heuristic algorithms works in such a way to determine the final solution, while only some existing solutions are visited. Each meta-heuristic algorithm under its own specific process attempt to find good solutions with no guarantee to be the greatest [START_REF] De Freitas | Optimization in software testing using metaheuristics[END_REF]. The optimization procedure for any type of meta-heuristic algorithm can be described as follows [START_REF] Khajehzadeh | A survey on meta-heuristic global optimization algorithms[END_REF]:

• Initializing the population in the search domain F by seeding the population with random values.

• Evaluating the fitness for each individual of the population.

• Generating a new population by reproducing selected individuals through evolutionary operations, such as crossover, mutation, and so on.

• Looping to step 2 until stopping criteria are satisfied.

Optimization problems arise in various disciplines such as engineering design, manufacturing system, economics etc. Identification of the optimal machining parameters for example is very important for reduction of machining costs, product quality improvement and increased productivity and profit. Therefore, the meta-heuristic algorithms have been increasingly used to further improve the solution of machining optimization problems with complex nature in many applications [START_REF] Yildiz | A novel particle swarm optimization approach for product design and manufacturing[END_REF]. It is reported that the meta-heuristic algorithms have been applied in machining because of their ability to deal with highly complex, non-linear, and multi-dimensional machining optimization problems [START_REF] Mukherjee | A review of optimization techniques in metal cutting processes[END_REF].

A review paper dedicated to machining parameters optimization by means of the metaheuristic algorithms is presented by Yusup et al. [START_REF] Yusup | Evolutionary techniques in optimizing machining parameters: Review and recent applications[END_REF]. As has been reported in the literature, three types of meta-heuristic-based search algorithms e.g. GA, SA and PSO have been mostly applied in the domain of the machining parameters optimization. However, in recent years there is an increasing trend in the application of other meta-heuristic algorithms such as ACO, ABC, IHSA, and CSA for solving machining optimization problems. In [START_REF] Madic | Comparison of meta-heuristic algorithms for solving machining optimization problems[END_REF] the authors have chosen four meta-heuristic algorithms, namely, RCGA, SA, IHSA and CSA for optimal combinations of different machining parameters for five case studies taken from the literature. IHSA has proved to be the most efficient meta-heuristic algorithm in terms of computational time and solution quality. In [START_REF] Said | A comparative study of meta-heuristic algorithms for solving quadratic assignment problem[END_REF], the authors applied GA, TS, and SA as Meta-heuristic algorithms for solving the Real life Quadratic Assignment Problem (QAP). Results showed that GA, TS, and SA algorithms have effectively demonstrated the ability to solve QAP optimization problems. Computational results showed that genetic algorithm has a better solution quality than the other meta-heuristic algorithms for solving QAP problems. Tabu-search algorithm has a faster execution time than the other meta-heuristic algorithms for solving real-life QAP problems.

TS has proved remarkably powerful in finding high-quality solutions to computationally difficult combinatorial optimization problems drawn from a wide variety of applications [1,[START_REF] Drezner | Tabu search model selection in multiple regression analysis[END_REF][START_REF] Glover | Tabu Search[END_REF]. The authors in [13] compared six meta-heuristic optimization algorithms applied to solving the traveling salesman problem: three classical approaches: GA, SA and TS, and compare them with three recently developed ones: QA, PSO and HS. It was shown that simulated annealing and tabu-search outperform newly developed approaches in short simulation runs with respect to solution quality, standard deviation of results and time needed to reach the optimum. SA finds best solutions, yet tabu-search has lower variance of results and faster convergence.

In this work, Tabu-search meta-heuristic optimization algorithm possibilities are explored in the framework of kernel regression optimization. Thus, a pool of variables or a search domain F is defined, and the best least squares regression kernel would be automatically selected by the algorithm. It is not the first time this technique is used. In [START_REF] Hage | Statistically validated and optimized tabu search estimation of cutting tool life in turning[END_REF], the authors combined tabu-search to regression (dubbed TS-REG) which yielded to a robust methodology that was compared to regression alone, artificial neural network and genetic algorithm and found that TS-REG is more reliable when dealing with large number of parameters. While in [START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF], tool life in turning has been estimated using a novel combination of tabu-search optimization and regression analysis (named TS-REG) and compared against regression analysis and artificial neural networks. The required time to select among 16 independent variables only the significant ones with highest R-square was 3 minutes (max), a minimum amount of time compared to other deterministic, Bayesian, and neural networks approaches. Comparing regression standalone with tabu-search combined with regression TS-REG, the latter select only the significant variables with the lowest p-values, the highest R-square better than deterministic, Bayesian and neural networks approaches, Comparing alternative regression with t-search combined with regression, the R-square was the highest, the p-values were the lowest and the mean absolute percentage error for the latter was the lowest. It is shown that Tabu-search optimized regression outperforms stepwise, backward or forward algorithms used in classical statistical packages offered by different commercial software [START_REF] Pacheco | A variable selection method based on tabu search for logistic regression models[END_REF].

The article starts with a review of the used taby search meta-heuristic algorithm and the selected regression methods, then tests the approach on a classical buckling problem. Later on, the optimized regression is tested on different complicated target functions, while exploring the effect of design for experiment selection on the results.

Tabu-search optimized regression method

In this section we illustrate the optimized regression using the Tabu-search algorithm technique. For the sake of simplicity, linear multiple regression will be considered for this work even though any nonlinear kernel-based regression could be used. We assume a set of inputs x with known outputs y. A regression tends to find the weights of the inputs such as:

Xb + ǫ = y (C.1)
with X a set of vectors or variables, where each X i = f i (x) ∈ F , and X will be named as the kernel of functions depending on x. ǫ is the bias of the regression. In our work we are considering a list of possible kernel functions F , consisting of a list of monomial functions, which can take negative and positive exponents. The weights b are found using the least squares minimization technique:

b = X T X -1 X T y (C.2)
and an estimation of the fitted function ŷ can be found using:

ŷ = Xb (C.3)
The error of the fitting is therefore the bias ǫ = yŷ which can also be written as:

ǫ = (I -H) y, (C.4)
with I the identity matrix and H = X X T X -1

X T the influence matrix, which is also the derivative of the estimation with respect to the measurements.

The regression results highly depends on the choice of the kernel functions X. Thus, an optimization scheme is used to find the best combination of monomial functions in X starting from a predefined pool of functions F . The selected optimization scheme is a multi-start Tabusearch algorithm. The Tabu-search algorithm starts from a random position in the domain, and defines an action to take to a neighboring element in the domain. A neighbor of the currently selected solution is defined as:

• Adding one element f * (x) from the pool of monomial functions F such as:

X * = X ∪ f * (x)
• Substituting one element from X by another element from the pool F

• Deleting one element X i from the current kernel solution X

The coefficient of determination R 2 of the regression is computed for the regression at each step or at each change of kernel neighborhood X. The neighborhood i considered as a better one if the coefficient of determination increased. However, the coefficient of determination increase is not always an indicator of the convergence of the data fitting in the regression. Thus, a composite cost function C is created to optimize the data fitting:

C = 0.5 × R 2 + (1 -MAP E) , (C.5)
where MAP E is the mean absolute relative error defined by:

MAP E = y -ŷ T • y -1 • 1 N , (C.6)
N being the number of measurements in the y vector. The algorithm acts by maximizing C while only conserving the statistically relevant monomial terms in the found regression. Tstudent tests are thus performed after each iteration to test the relevance of each one of the found monomial terms.

Once a kernel is selected by the Tabu-serach algorithm, the outliers can be removed if having a high leverage. For example a measurement i is considered as an outlier if:

H ii > 3p N , (C.7)
p being the number of selected vectors in the kernel X. The outliers are therefore removed before attempting any regression fitting.

To avoid cycling and rechecking tested neighborhoods, the Tabu-search algorithm makes use of the Tabu list L. If a variable X j taken out of the pool F decreases the cost function C instead of increasing it, the variable X j is added to the Tabu list. If a variable is listed in the Tabu list, it is taken out of the pool of possible kernel functions F . The Tabu list has however a fixed length L. Once the Tabu list reaches its full length, the oldest variable in this list is deleted and enters back the list of possible kernel functions F . The Tabu list helps minimizing the cycling effects of testing some already tested neighborhoods. However, the variables reducing C are not penalized forever, since they may reenter the list F after a while, considering that the best currently selected combination of variables X might have changed, and the variable X j can potentially increase the cost function C with this new combination.

The presented algorithm is first applied on numerically generated data before attempting to fit the experimental data.

Numerically generated buckling data

Classical buckling of straight columns is a challenging issue for regression and fitting due to the high non linearities involved with respect to the parameters of the column. Classical mechanics relations gives the critical buckling load around the z axis P cr of straight isotropic columns with rectangular cross section as a function of the column parameters as:

P cr = π 2 EI z L 2 (C.8)
With E being the modulus of elasticity of the column, L the length of the column, I z its second moment of area around the z axis defined by:

I z = 1 12 bh 3 (C.9)
With b the dimension of the cross section parallel to the z axis and h the dimension normal to the cross section as illustrated in figure C.1. The relations clearly illustrates the column length L to the power -2 and its dimension b to the power 3. We compute the value of the buckling load for E = 200GP a and various values of L, b and h, obtaining therefore numerical data for about 2000 points.

The Tabu algorithm illustrated in section 1.2 is later on used to identify the best polynomial function that fits the obtained data. The initial pool of monomial F was initiated to all possible combinations of parameters having a power of -3 to +3, including the zeros, which leaves us with 194 potential monomials and the constant term. The algorithm converges to obtain therefore the monomial combination defined in table C.1, with R 2 = 1 and MAP E = 1.

One may note from table C.1 that the identified constant for the correct monomial function is exactly equal to π 2 × E. All the other terms are negligible with a weight ratio less that 10 17 with respect to the correct monomial. These weights can be attributed to numerical errors. Now considering the absolute critical buckling load in mode 1, the equation of buckling becomes even more complicated even using classical mechanics:

P crz = π 2 EIz L 2 for h ≤ b P cry = π 2 EIy L 2 for b ≤ h (C.10)
with the second moments of area I z and I y defined by:

I z = 1 12 bh 3 I z = 1 12 hb 3 (C.11)
Eventually the absolute minimum buckling load would lead to the first failure for P cra = min P crz ; P cry . Using the Tabu-search algorithm and attempting to fit P cra in one polynomial function is not a good idea since a step function needs theoretically an infinite number of monomials to correctly represent it. However such fitting may lead to a good fitting results at the expense of increasing the number of monomial terms and the computation time. However, fitting P crz and P cry in two separated polynomial functions and later on taking the solution as the minimum of the two polynomial results P cra = min P crz ; P cry would lead to the exact solution of the problem. This fact would inspire our next part of the work fitting therefore each mode of buckling with a different polynomial function.

Artificially generated data fitting

In this section we illustrate data fittings of three different numerically generated data set using 4 variables x 1 , x 2 , x 3 and x 4 . The data sets are first generated with a polynomial easy functions consisting of four monomials. Later on a more challenging data set is generated with 20 different monomial terms generated using random powers and monomial coefficients. Finally, a challenging function is used with monomial terms taking possibly non-integer powers. The results are shown for different meshing used for the input variables x i , all the considered input are included in [0; 10]. The considered pool of functions F consists of all possible combination of monomials with a degree lower than 4.

An easy monomial function

In this section we illustrate the usage of Tabu-search optimized regression to fit the polynomial target function:

y 1 = 5x 2 1 + 0.26x 1 x 3 3 + 0.53x 2 x 3 + 26x 4 (C.12)
The tested function y 1 consists apparently of an easy test for the Tabu-search algorithm. The test is performed on 501 generated data points and split as follows: 80% (401 points) are used for estimation/training while the remaining 20% (100 points) are used for prediction/evaluation of the regression. The test is performed on three different meshes for x i : 2. A mesh consisting of all possible combinations of selected values for each one of the x i , as usually performed in an experimental plan, named mesh 2.

3. A uniform mesh without repetition of nodal values, named mesh 3.

The regression results for mesh 1 using random values for each x i yields the following least-square fitted equation (C.13) while enforcing the presence of a constant in the polynomial regression: The difference between y 1 and ŷ1 is the presence of two terms, a constant and the interaction between x The regression results for mesh 2 using all possible combination of values from a selected pool for each x i yields the following regression when enforcing the presence of a constant in the polynomial regression: correlation may exist between the x i and these variables are therefore not totally independent. For instance, x i may be derived from x j using x i = a • x j + b with a and b are two constants. Thus, the Tabu-search have found the best regression to fit the data, taking into consideration the correlation between the x i , which can't be considered as independent variables, since we can get any x i form scaling and translating the others. The goodness of fitting shown in figure C. 4 shows therefore an excellent fitting results for both the training data set and the testing or evaluation data set.

ŷ1 = 5x
ŷ1 = 5x
The trial test in this section on the target function y 1 illustrates the ability of the Tabusearch optimized regression to identify the exact functions y 1 when using independent input variables x i . The use of correlated input variables would yield good fittings and evaluation, without having the exact initial target function.

A complicated polynomial function

The target test function in this section y 2 is polynomial, however made of 20 different monomial terms with random integer exponents between 0 and 4, while keeping the degree of each One may note that for mesh 1 using radom values for x i , the Tabu-search optimized regression with 20 multi-starts was able to find the exact 20 monomial terms of y 2 6 times out of 20 multi-starts from random initial solutions. The identified function corresponded to the 20 monomials with their exact coefficients a j , along with few other monomial terms with coefficients several order of magnitudes lower than a j . The composite cost function C is illustrated in figure C.8 for the 20 multi-starts. Using mesh 2, which is common in experimental testings, the Tabu-search was able to find 14 exact monomial terms only, even after 20 multistarts, while still finding excellent correlations for both training and testing data sets as shown 6. Finally, for a finite element mesh 3, the Tabu-search optimized regression was able to find excellent fittings for both training and evaluation datasets using only 3 monomial terms due to the correlation between the non-independent input variables.

A challenging target function

In this section, we use a challenging target function y 3 to fit using the Tabu-search optimized regression. y 3 consists of five terms with a mix of integer and non-integer exponents:

y 3 = a 1 x 0.5 1 x 2 3 x 4 + a 2 x 0.1 3 x 4 + a 3 x 0.3 1 + a 4 x 3 1 x 2 + a 5 x 1 x 0.2 3 , (C.17)
with a j random coefficients of y 3 terms ∈ [0; 10]. y 3 is a challenging target function since it contains only one term from the pool F , while all other terms are taken from outside the pool of functions F . The regression is performed on the three meshes as performed in the previous sections. The goodness of fitting curves for meshes 1, 2 and 3 are respectively illustrated in figures C.9, C.10 and C.11.

The illustrated results in figure C.9 shows accurate fittings even though the considered y 3 terms are not initially found inside the F pool of variables. The optimized regression was able 11 both shows excellent fitting results for both training and evaluation data sets. However, the regression using mesh 2 was able to identify the monomial term x 3 1 x 2 included in y 3 and F , but the one using mesh 3 never found that term even though figure C.11 show better results than C.10. This is explained by the non-independent nature of x i while using linear 1D finite element mesh to evaluate their values. The optimal regression using mesh 3 consisted of only 3 monomial terms, while the one found while using mesh 2 consisted of 21 monomial terms.

Conclusion

In this work we illustrate the use of Tabu-search optimization as a mean to find the best regression kernel, starting from a large pool of variables F . The ability to optimize the regression kernel improves the fitting abilities as well as the predictive ability of a regression as illustrated in sections 1.3 and 1.4. The application of the algorithm is illustrated on buckling of straight columns and on three different sets of artificially generated data. Four parameters were used as regression input with a large pool of monomials is initially considered. Different input parameters are also tested and the benefits of using of independent input parameters is shown. The Tabu-search optimized regression was also able to find the best representation of the data when using correlated input variables x i , with a lower number of terms than the ones used to create initially the target function y. The final fits illustrate low mean absolute percentage errors, high R 2 and a high predictive ability as illustrated in section 1.4. It is also important to highlight that only statistically relevant terms are derived by the designed algorithm.

Hybrid modeling combining the proper generalized decomposition and model learners 2.1 Introduction

Data-driven techniques are becoming increasingly popular in many field, engineering is not an exception. Although science always relied on data, it is only until the beginning of the 21th century that the eruption of data-driven modeling took place [47]. This phenomenon was motivated by both the increase of data availability and sharing, as well as the increase in computing power. Many algorithms are developed based on linear or non-linear regressions, decision trees, random forests, neural networks... All these models tends to uncover a relation between an input and an output, while neglecting the physical phenomenon involved. As a result, the modeling interpretation becomes too complicated if not impossible in many situations.

On the other hand, the development of simulation techniques for solving PDEs made a tremendous progress in the filed of developing real-time simulations on the fly, with minimal calculation time thanks to the domain decomposition techniques [45,[START_REF] Ghnatios | A reduced-order model manifold technique for automated structural defects judging using the pgd with analytical validation[END_REF][START_REF] Perez | Sensitivity thermal analysis in the laser-assisted tape placement process[END_REF]. Moreover, many model reduction techniques perform an a posteriori model reduction by compressing the available data using the most relevant information [11,20,[START_REF] Leonenko | On the resolution of the fokker-planck equation using a high-order reduced basis approximation[END_REF][START_REF] Patera | Reduced basis approximation and a posteriori error estimation for a boltzmann model[END_REF]. The available techniques allow, after a pre-possessing of f line calculation, to use the computed parametric solutions in real-time in what is called the online phase for real-time optimization and control [START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF][START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven of thermal processes[END_REF]. The model reduction techniques were successfully used in many optimization problems, since the gradients with respect to the process parameters are readily available in the treated prob-lems [START_REF] Ghnatios | Simulation avancée des problèmes thermiques rencontrés lors de la mise en forme des composites[END_REF][START_REF] Perez | Sensitivity thermal analysis in the laser-assisted tape placement process[END_REF].

Using model reduction techniques or reduced order basis methods lead in general to a good approximation, which tends to the finite element approximation when increasing the size of the projection basis [8,9,[START_REF] Porsching | Estimation of the error in reduced basis method solution of nonlinear equations[END_REF]. However, model reduction involves truncation errors in general, on top of the modeling error. Therefore, a room for improving model reduction techniques is possible. In this work, we attempt to use advanced machine learning algorithms to correct modeling and simulation errors. Moreover, and despite the impressive progress in computational power, many problems remain intractable using classical model reduction methods [START_REF] Ghnatios | Advanced separated spatial representations for hardly separable domains[END_REF], or remains too intrusive for an inclusion in an existing commercial software [START_REF] Ghnatios | Spuriousfree interpolations for non-intrusive pgd-based parametric solutions: Application to composites forming processes[END_REF]. Despite numerous developments in both intrusive and non-intrusive manners [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF], the calculation of solutions using model reduction techniques faces major challenges when changing a parameter induces drastic changes in the solution map [START_REF] Ghnatios | Towards parametric rtm process: The interpolative mapping[END_REF].

With the recent development of artificial intelligence and machine learning algorithms, some combinations of model reduction techniques with machine learning are addressed to account for error reduction and reliability improvements, either by error correction [4] or by improvement of modal selection using the quantities of interest [START_REF] Hartman | A deep learning framework for model reduction of dynamical systems[END_REF].

In this work, we aim to combine the physical models generated in the reduced basis, with machine learning algorithms, to create a physically informed machine learning algorithm. Eventually, with the multi-parametric optimization involved in any machine learning technique, involving many iterations and derivatives computations, any classical simulation algorithm will fail to create a physically informed neural network due to the multiple simulations required at each iterations. Therefore, we leverage the availability of a parametric simulation using the PGD model reduction technique, to learn the physical problem parameters along with the machine learning algorithms hyper parameters. The selected application for this work is a non-linear visco-elastic biphasic material. In fact, thanks to the recent progress in simulation technology and computing power, the mechanical behavior of biological tissues is nowadays one of the most active research topics. However, often biological tissues are biphasic by nature, which renders their effective modeling and simulation a challenging issue, even with the impressive progress achieved recently. Human cartilage, for example, is a biphasic material, where the fluid pressurization is believed to be the main load carrying phenomenon [START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]. The microstructure and the fluid-solid interactions in such materials is still not fully understood [START_REF] Moore | Tribological rehydration of cartilage and its potential role in preserving joint health[END_REF][START_REF] Pawaskar | Modelling of fluid support inside articular cartilage during sliding[END_REF]. Previous works attempted to model soft permeable matter as a superposition of a fluid force contribution to an elastic solid one [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]. The previous modeling do not take into consideration the possibility of having a non-elastic behavior of the solid counterpart of the biphasic material. In this work, we apply the previous model to the atomic force microscopy nanoindentation of a biphasic material with a large thickness, when compared to the indentation depth.

First of all, we will revisit the biphasic soft material modeling in section 2.2. Later on, we use the PGD to simulate the indenter behavior, and explain how the obtained reduced order model results can lead to an improvement of the model in section 2.3. Once the PGD results are available, we leverage a multidimensional parametric solution and the experimental results shown in section 2.4, to identify the permeability of a biphasic hydrogel using an elastic behavior of the solid component of the hygrogel in question, in section 2.5. The results being none fully satisfactory, a correction of the model is performed twice, once using an analytical modeling, and another time using a fully coupled, physically informed, machine learning algorithm, in sections 2.5 and 2.5 respectively. It is finally shown that the modeling coupling the simulation

PGD solution of the multidimensional Darcy equation

The problem is solved in a cylindrical domain of height H and radius R with the following boundary conditions:

     v = 0 at z = 0 v = (0, 0, -U) at z = H P = 0 at r = R (C.27)
which consists of a zero velocity on the bottom support, the indentation velocity U at the top of the domain and zero pressure at r = R. The shown modeling results in a steady-state simulation problem. Since we aim to identify the permeability of the domain, one will have to solve the problem several times during the optimization algorithm for each value of the permeability. Moreover, at each time step, the same problem should be solved again, at each change in the value of the parameters U, R or H. On top of this complication, computing the gradients numerically during an optimization problem will require several more solutions of the direct differential equation, when using a classical approach. Therefore, disposing of a parametric solution of the form P (x, y, z, H, U, R, K) becomes a leverage which we will use during the optimization problem.

In this simulation, we adopt a quasi-static approach while solving the Darcy equation, assuming the equilibrium is reached at each time step. Therefore the Darcy equation is written by:

v = -µ K ∇P (C. 28
)
where P is the pressure, µ the fluid's viscosity and K the domain permeability.

Another challenge is faced in this simulation, using a domain of initial height H 0 = 3mm, while the indentation depth has the order of magnitude of few micro meters. Therefore, all the pressure variation exists only at the surface of the domain, within a thickness of few micrometers. Such problem is complicated using classical simulation techniques since it requires an extremely refined mesh in the thickness direction near the indented surface. The result is a large number of required degrees of freedom leading to unrealistic calculation time using classical calculation techniques. However, using the PGD technique, we solve the problem using an in-plane-out-of-plane decomposition. Such decomposition results in solving the problem in (x, y) domain separately from the thickness z domain. Thus refining the mesh in the thickness direction as much as needed is possible without compromising the computation time [6,[START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF]. We may also use the PGD ability to solve problems using geometrical and configuration parameters as extra coordinates of the problem [START_REF] Ghnatios | Poroelastic properties identification through mico indentation modeled by using the proper generalized decomposition[END_REF]. This approach leads to a large problem dimensionality, but is tractable using the PGD since it allows circumventing the curse of dimensionality by using the separation of variables [45]. The solution is a parametric 7D solution of the pressure P defined in the following form:

P = N i=1 X i (x, y) • Z i (z) • U i (U) • L i (K) • R i (R) • H i (H) (C.29)
We use in this simulation 100 000 nodes in the out-of-plane thickness direction z, logarithmically distributed, with an extremely refined mesh on the top surface of the domain, to avoid oscillation near the top surface, whereas we use 363 nodes in the in-plane domain. The permeability domain consists of 10 000 nodes logarithmically distributed to capture the behavior at low permeabilities. Moreover, we use 10 000 nodes in the thickness dimension H, 5000 in the radius dimension R and 1000 nodes in the imposed velocity dimension U. The PGD simulation resulted in a 3D pressure field for any indentation velocity U, domain dimensions H and R and permeability K. Therefore the results can be used on the fly, online, at any time step t, once R, H and U are computed from measured variables, to identify the domain permeability K for example. The simulation is performed within 2 hours on a normal portable PC, core i7 with 8 Gb RAM. The result obtained is the equivalent of 5 × 10 14 3D simulations, each with 36 300 000 degrees of freedom. The solution for one combination of these parameters is illustrated in figure C.13. Readers unfamiliar with the construction of parametric solutions and domain decomposition using the PGD are directed to [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF][START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF] and the numerous references therein. What follows consists of a small introduction of the solution of the problem using the PGD. We first assume the function P as having a separated form known until an order n, and we attempt to define P n+1 using:

P n+1 = n i=1 X i (x, y) • Z i (z) • U i (U) • L i (K) • R i (R) • H i (H) + 6 j=1 R j , (C.30)
with R j , j = 1 • • • 6 being the new functions to be identified defining the unknown terms in P n+1 , for instance the functions of (x, y), z, U, K, R and H. Since the geometrical dimensions are extra parameters of the PGD solution, the (x, y) domain is mapped into an apparent domain (ξ, ζ), a circle of center (0, 0) and radius 1. The thickness domain z is mapped into φ ∈ [0; 1]. The mappings are defined by:

     x = R • ξ y = R • ζ z = H • φ (C.31)
The Jacobian of the transformation becomes therefore:

J =    R 0 0 0 R 0 0 0 H    (C.32)
Finally, the integral form of the problem becomes after replacing the (x, y, z) coordinates with the apparent domain ones:

ξ ζ φ KP * µ H ∂ 2 P ∂ξ 2 + H ∂ 2 P ∂ζ 2 + R 2 H ∂ 2 P ∂φ 2 dξdζdφ = 0 (C.33)
Therefore, the exact notation of P is now:

P n+1 = n i=1 X i (ξ, ζ) • Z i (φ) • U i (U) • L i (K) • R i (R) • H i (H) + 6 j=1 R j , (C.34)
Eventually the problem becomes non linear and its solution is obtained using a fixed point iterative algorithm. To compute R k , we first define the test function k noted P * k using:

P * k = R 1 × • • • × R a k st × • × R 6 (C.35)
Then we replace into the weighted integral for of the problem leads to a 2D problem in case k = 1, and a 1D problem for all other values of k. Once the fixed point iterative algorithm is converged for P n+1 , the newly found functions R j will be used to compute P n+2 and so on until the convergence of the residual of the governing differential equation of the problem defined in equation (C.23). One may note that using the PGD, the gradient of P with respect to the process parameters are found in a straightforward manner, simply by deriving the function of the permeability L i (K). For instance:

∂P ∂K = N i=1 X i (ξ, ζ) • Z i (φ) • U i (U) • ∂L i (K) ∂K • R i (R) • H i (h) (C.36)

Experimental study

The indentation of the hydrogel poroelastic material was performed at 5 different probe velocities ż. The domain consisted of a relatively thick hydrogel, of thickness H 0 = 6mm. The results are obtained using atomic force microscopy [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF], and illustrated in figure C.14. As illustrated in figure C.14, the indentation depth has the order of magnitude of few micrometers. Therefore no sensible variation in the material properties are expected to appear in the studied domain.

Identification of the permeability and comparison of two models

In this section, we solve the minimization problem defined by:

(K; p solid ) = argmin (F exp -F sim ) 2 (C.37)
where p solid is whatever property of the solid component to be identified from the minimization problem, F exp is the experimentally measured force in the spring and F sim is the simulated force, which is defined as the sum of the contribution of the solid component and the fluid component of the reaction force. Two scenarios are considered for this work:

Modeling using a fluid in an elastic solid mesh

In this section, we solve equation (C.37) considering only the contribution of the elastic part of the solid component as defined in equation (C.26). Therefore, equation C.37 is a function of only one parameter, the permeability of the domain K since the solid properties are known for instance: The minimization is performed using a classical Newton algorithm using the cost function to minimize C = (F exp -F sim ) 2 , and the unknown parameter to identify X = K are updated at each iteration using:

F sim = F elastic solid + F f luid = F spring@ ż≈0 + F f luid (C.38) Indentation depth w (m) Indentation force F (t) (N) 0 0.5 1 1.5 2 ×10 -6
X i+1 = X i -∇C(X i ) -1 C(X i ) (C.39)
The Newton algorithm leads to the permeabilities shown in figure C.15.

As one may note from the results illustrated in figure C.15, that a small transition phase in the first 0.6 µm exists but later on the results are steady as a function of w. These results are expected since the indentation depth is small with respect to the thickness and therefore w should not affect the permeability values. On the other hand, one may notice a dependency of the permeability on the indentation velocity ż. Such result is unexpected and is a possible sign of the presence of a viscoelastic solid and/or other factors that may affect the material reaction. These factors are clearly linear with respect to the indenter velocity ż as illustrated in figure C. 16. In fact, figure C.16 illustrates, as a function of ż, the average steady-state permeability value for each indentation velocity, from w = 0.6µm until the end of the indentation experiment. To explain the speed dependence of the apparent permeability shown in figure C.15, another term, proportional to the rate of deformation ǫ or ż, can be added to the equation C.38.

A nonlinear correcting term

In this section, we compare the experimental results to the simulation of a biphasic material with a non-elastic correcting term in the equilibrium equation. Using the results shown in figure C.16, one might assume that lumping the effects of material viscoelasticity into the permeability term generates inaccuracies that reflect as an apparent permeability that increases with speed. Assuming a constant permeability, we added another term that represents a rate to deformation dependent response of the solid polymer network into the equilibrium equation. Therefore, and since ż is proportional to ǫ, we suggest the following correction:

σ z = σ elastic + γ • ǫz (C.40)
Where γ is the correction factor representing the non-elastic response of the solid, σ z is the normal stress on the solid and ǫ z is the normal strain in the z direction. The dependency of γ • ǫz should be linear as a function of ż using this model, since w comparable to z [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF], and ǫ z is linear as a function of w. The cost function (C.41) to minimize is now written as:

(K; γ) = argmin F exp -F sim 2 (C.41)
The force in the solid component is now therefore defined as:

F solid = F spring@ ż≈0 + γ • ǫz • π • R 2 (C.42)
The cost function now has 2 unknowns to be identified. Therefore we use the Gauss-Newton algorithm to identify the two missing properties. Thus, we first define the following notation:

f v (K, γ) = -F spring@ ż=vµm/s + F spring@ ż≈0 + γ • ǫz • π • R 2 + A T op P • dA (C.43)
A T op being the indented area on the top surface of the indenter. We may redefine the cost function to minimize by C such as:

C = f 2 2 + f 2 8 + f 2 40 + f 2 80 , (C.44)
where f i is the error function defined by equation C.43 for the given indentation speed ż = i. Now the optimization problem is rewritten as:

(K; P solid ) = argmin [C] (C.45)
Problem (C.45) is solved using Guass-Newton algorithm. The expected result is a constant permeability K. Indeed, the identification yields a constant value of the permeability K ≈ 10 -16 but a negative viscosity γ as shown in figure C.17. This shows that the constant permeability is in fact generating higher contribution to the reaction force than the real experimental one, and a negative term shall be added to go back to the correct value.

An informed neural network

In section 2.5, the modeling error was corrected with a non-linear term using a visco-elastic approach in the modeling of the solid mesh behavior. However, the results were not physical and the final model error using a constant permeability was still large. The modeling therefore, although did enhance the end outcome, remained sub-optimal. With the presence of this human ignorance with respect to the behavior of the indented specimen, an informed artificial intelligence neural network fitting the modeling bias, becomes an appealing solution. In this section, we combine the fluid-reaction force model of section 2.5, to a shallow neural network defined in section ??, defining therefore a hybrid twins approach, coupling the simulation results to a data-driven modeling technique, enhancing therefore the numerical simulation with experimental data. Eventually one can choose any fitting technique with enough parameters to represent the "ignorance", as neural networks are one of the fitting techniques. The Figure C.17: The identified viscosity of the solid part of the biphasic material as a function of ż author choice is motivated only by the convenience and ease of use of the considered technique.

In this section we use a shallow neural network made of 6 neurons in the first layer with Sigmoid activation function and one neuron in the output layer with a linear activation function. The model representation is given by figure C. 18. In figure C.18, the first block to the left represents the model depicted in section 2.5, which output is the model estimation of the fluid force F f luid using as hyper parameter the permeability K only. We define this modeling error in the section by E using:

E = F exp -F spring@ ż≈0 -F f luid , (C.46)
where F f luid is computed using the PGD results as F f luid = A T op P • dA. Equation (C.46) represents the error the model defined in section 2.5, where the simulated reaction force was defined as F sim = F spring@ ż≈0 + F f luid . Therefore E is the error of the elastic model. We will therefore attempt to use a neural network to model E. We define Y 1 as the output of the neural network, which should consists an estimation of E. We also define Y as an estimation of the experimental reaction force such as:

Y = Y 1 + F spring@ ż≈0 + F f luid (C.47)
With this approach, eventually the fitting Y 1 of E with the neural network requires the knowledge of F f luid , the output of the PGD algorithm, and therefore a fitting of K. Therefore, a simultaneous fitting of both K and the hyper parameters Θ i of the defined neural network is the most suitable approach for this process.

To compute the problem's hyper parameters, the following cost function is defined: i consists of a regularization term to avoid over fitting, the parameters Θ i are the neural network hyper parameters and λ is the regularization term. Later on, the minimization algorithm identifies the hyper parameters using:

J = 1 2 E 2 + 1 2 (Y 1 -E) 2 + λ 1 2 n i=1 Θ 2
(K, Θ i ) = argmin J (C.49)
To achieve convergence, we first converge the fluid force using a regularized Newton algorithm, as the gradient with respect to the fluid permeability K is much larger than the gradient with respect to the neural network hyper parameter. For instance, we minimize the cost function J * defined by:

J * (K) = 1 2 E 2 , (C.50)
using an iterative regularized Newton update algorithm defined by:

K i+1 = K i - J * i ∂J * i ∂K + c 1 J * i (C.51)
The constant c 1 is a regularization constant to avoid numerical explosion of the algorithm, used equal to 1 in this work. The value of K converges within several iteration to 1.98 × 10 -16 P a.s, a quite comparable value to the previous modeling effort.

Once the value of K estimated, it will be used to the initiation of the minimization problem defined in equation C. 49. The problem is minimized using Levenberg-Marquardt gradient decent algorithm with a non-uniform learning rate. For instance, the learning rate of the hyper parameter K is lower than the learning rate of the neural network hyper parameters Θ i or biases b i , one bias per layer in the neural network. For training the neural network, the data is split into 70% for training set, 15% on the cross-validation set and 15% for the test set. The hyper parameter λ is fitted to the cross validation set. A batch gradient decent minimization is used to accelerate the convergence process.

The final estimated permeability K = 3.86 × 10 -16 , and the performance of the neural network fitting algorithm is illustrated in figure C. 19. The mean square error is MSE = 0.48% on the best performance on the cross validation set. As expected, the permeability is larger than the one identified in section 2.5, which was a bit overestimating the fluid reaction forces.

Discussion

Many factors may affect the modeling. Moreover, simulating all the local behavior around the indenter may become unacceptably expensive form a computational time point of view, taking into consideration the local effects and the thickness of the domain with respect to the indenter. Therefore the adopted suitable approach is an empirical modeling of the error with by including a correcting factor in the equilibrium equation proportional to the strain rate ǫ. In a first attempt, the parametric simulation of the Darcy equation is performed using the PGD and the Gauss-Newton optimization algorithm is used to identify the parameter γ, as identified in section 2.5. The obtained result γ is linear as a function of ż, which was expected from the reported error is figure C.16. The obtained permeability is constant therefore as a function of ż after the first transient region K ≈ 10 -16 , which is in fact less than the expected values shown in figure C.16. However, the final results estimates a single value of the permeability.

On a second approach, a physically informed neural network is used to correct the modeling error using a set of 25 hyper parameters (18 weights and 7 biases) to fit simultaneously the permeability using the model developed in section 2.5 and the neural network parameters. The algorithm converges within few minutes on a normal portable PC, leading to an excellent performance of the model. The estimation modeling of the error E is illustrated in figure C.20, with a much better performance than the visco-elastic solid modeling explained in section 2.5.

One may note that the availability of the PGD parametric solution in the separated form defined in equation C.29 made the fully coupled minimization of J possible, and therefore possible the identification of the permeability K along with the hyper parameters of the neural network. Otherwise, each iteration of the Levenberg-Marquardt algorithm would require at many resolution of the direct PDE defined in equation C.23, which least to prohibitive and unrealistic calculation times.

Conclusion

In this work, we derived an expression of the error with respect to the previous work modeling biphasic materials as a combination of an elastic solid and a pressurized fluid. The suggested error modeling first identifies a strain-rate term and a parameter γ linear with respect to the indentation speed ż. Later on, a neural network is used for a superior fitting of the modeling error. The identification is performed by comparing a parametric simulation to an experimental atomic force microscopy nanoindentation of a biphasic soft hydrogel with a thickness several order of magnitude larger than the indentation depth. The simulation is performed by using an in-plane-out-of-plane decomposition allowing an impressive refinement of the mesh in the thickness direction and therefore identifying with high fidelity the pressure distribution inside the studied part.

The used hybrid-twins approach is only possible when all the solutions of the direct problem are readily available, which is the case in a model reduction framework, for instance the PGD framework. The absence of the parametric multidimentional solutions leads to prohibitive calculation times during the optimization process.

3 Non Parametric probabilistic method for PGD solutions

Introduction

Nowadays, simulation-based decision making in engineering and sciences is widely accepted and used in predicting material and parts behavior. Many partial differential equation (PDE) models are widely accepted and used in the industry, while constantly being improved to increase their fidelity. On the other hand, many advanced methods for solving PDEs are being designed to improve solution speed, aiming for real-time applications [12,45,[START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF][START_REF] Grepl | efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Nguyen | An efficient reduced-order modeling approach for nonlinear parametrized partial differential equations[END_REF]. These techniques are commonly known by model order reduction techniques (or simply model reduction techniques), mainly aiming at constructing a reduced basis for the solution of PDEs [START_REF] Ghnatios | A reduced-order model manifold technique for automated structural defects judging using the pgd with analytical validation[END_REF][START_REF] Perez | Sensitivity thermal analysis in the laser-assisted tape placement process[END_REF].

Regarding model order reduction techniques, the methods are divided into two large categories: (i) "a priori" model reduction techniques, which aims at constructing a solution in a reduced basis separated form, before having any knowledge of the full order model or also the high definition model (HDM) solution [START_REF] Ghnatios | First steps in the space separated representation of models defined in complex domains[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]. (ii) "a posteriori" model reduction techniques, where the reduced order basis of projection is constructed from previous HDM solutions computed using efficient greedy algorithms to construct as accurate as possible reduced order basis [11,20,[START_REF] Leonenko | On the resolution of the fokker-planck equation using a high-order reduced basis approximation[END_REF][START_REF] Patera | Reduced basis approximation and a posteriori error estimation for a boltzmann model[END_REF]. The PGD method appears is known to be an "a priori" model reduction technique which is widely used nowadays to compute solutions in a fairly separable domain [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF][START_REF] Ghnatios | Advanced separated spatial representations for hardly separable domains[END_REF][START_REF] Ghnatios | Spuriousfree interpolations for non-intrusive pgd-based parametric solutions: Application to composites forming processes[END_REF]. Using model reduction techniques or reduced order basis methods lead in general to a good approximation, which tends to the finite element approximation when increasing the size of the projection basis [8,9,[START_REF] Porsching | Estimation of the error in reduced basis method solution of nonlinear equations[END_REF]. However, it is known that model reduction techniques in general involve a truncation error (either being an "a priori" or "a posteriori" techniques) as well as modeling errors, inherited from their full-order models counterparts, also known as high definition models (HDM) [START_REF] Ghnatios | Simulation avancée des problèmes thermiques rencontrés lors de la mise en forme des composites[END_REF]. Therefore, a margin of improvement of the solution do exist.

On the other hand, with the current development of machine learning and data-driven techniques, fitting data using surrogate models based on regressions, neural network is becoming widely available and increasingly popular [START_REF] De Weg | Neural network-based surrogate model for a bifurcating structural fracture response[END_REF][START_REF] Ghnatios | An efficient tabu-search optimized regression for data-driven modeling[END_REF][START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF][START_REF] Xu | Solution of diffusivity equations with local sources/sinks and surrogate modeling using weak form theory-guided neural network[END_REF]. However, these surrogate models do not impose naturally physical conditions like the conservation of mass or energy. Therefore, some other PDE fitting techniques are being developed, based on satisfying thermodynamic constraints, using methods like the GENERIC formalism [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF][START_REF] Gonzalez | Thermodynamically consistent datadriven computational mechanics[END_REF] or the General Equilibrium for Non-Equilibrium Reversible-Irreversible Coupling [START_REF] Grmela | Multiscale thermodynamics[END_REF][START_REF] Ottinger | Beyond equilibrium thermodynamics[END_REF]. However, often engineering applications data is scarce and expensive to generate in sufficient quantities to solely rely on data models. Therefore, to leverage recent developments of artificial intelligence and machine learning algorithms, combinations of model reduction techniques with machine learning are built to address the modeling and truncation errors, either by error correction [4] or by improvement of modal selection using the quantities of interest [START_REF] Hartman | A deep learning framework for model reduction of dynamical systems[END_REF].

Another appealing approach for creating data-driven surrogate models while leveraging previously well-established models is the "Digital Twins" approach [4,47]. Digital twins are built based on the correction of previously established PDEs using surrogate models, to learn new physics or establish new models quantifying the error [START_REF] Ghnatios | A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials[END_REF][START_REF] González | Learning corrections for hyperelastic models from data[END_REF], but also to incorporate the variability of the experimental results in the HDM or reduced order models solutions [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF][START_REF] Soize | Probabilistic learning for modeling and quantifying modelform uncertainties in nonlinear computational mechanics[END_REF]. Multiple other works aimed to quantify the errors induced by modeling, either for the HDM solutions [START_REF] Ghanem | Stochastic Finite Elements: A spectral approach[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A spectral approach[END_REF][START_REF] Maitre | Spectral Methods for Uncerainty Quantification with Applications to Computational Fluid Dynamics[END_REF][START_REF] Schueller | Computational methods in stochastic mechanics and reliability analysis[END_REF], or for using model reduction techniques relying on reduced order basis (ROB) [36,[START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamicsrecent theoretical advances[END_REF][START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. Moreover, a non-parametric probabilistic method (NPM) for error quantification was introduced in [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF] and extended to µ-parametric ROB methods in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. The NPM method has the advantage of quantifying uncertainties independently of their sources, either modeling errors, nonlinear modeling, truncation errors... This is performed by replacing the deterministic reduced basis by a stochastic counterpart [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. However, to the best knowledge of the authors, non of these techniques where adapted to be used for the Proper Generalized Decomposition (PGD) framework, but are rather used for the Finite elements HDM approach, or for a µ-parametric ROB solution [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF].

In this work, we aim to extend the NPM method presented in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF] to the PGD framework. The extension results in a stochastic solution based on a stochastic reduced basis, who's variations are restricted to the Stiefel manifold of the original PGD solution. The variations are controlled by several hyper parameters identified through solving an optimization problem on the fly, as the method do not require any direct solution of the PDE modeling the problem during optimization.

First of all, we will revisit NPM approach in section 3.2 of this work. Later on, we introduce the PGD simulation and experimental measurements of the automated tape placement, the application selected for this work, in section 3.3. In section 3.4 we illustrate the NPM results applied on the PGD solution obtained in section 3.3. Finally the article is wrapped up with discussions and conclusions in section 3.5.

NPM method for PGD solutions

Short review of NPM for reduced order basis

Let's consider for illustration purpose a PDE having the following semi-discrete formulation: The variable [q(t)] is the solution of the PDE equation depicted in equation (C.52) in the reduced coordinate system, after projection on the reduced order basis, and solving as a function of the time. The NPM method consists of replacing a reduced order basis [V ] ∈ R (N,n) by a stochastic counterpart [V] ∈ R (N,n) , where N is the number of degrees of freedom in the HDM model, and n is the dimension of the built reduced basis. Considering a set of constraints or Dirichlet boundary conditions defined by a matrix [B] ∈ R (N,N BC ) , the reduced basis [V ] should satisfy the following orthogonality and boundary conditions constraints:

[V ] T [M][V ] = [I n ] [B] T [V ] = [0] (C.54)
where [I n ] being the identity matrix in a domain having the dimension of the reduced basis n, and [M] ∈ [N, N] the mass matrix of the problem in question. Therefore, the stochastic reduced basis should satisfy the same constraints as the initial reduced basis:

[V] T [M][V] = [I n ] [B] T [V] = [0] (C.55)
Using NPM, the stochastic basis [V] is constructed using the maximum entropy principle, while variation of [V ] to construct [V] should be performed on the tangent vector T v space defined by [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]:

T v = [Z] ∈ R (N,n) : [V] T [M][Z] + [Z] T [M][V] = [0] (C.56)
with [Z] being therefore the first derivative of [V ] using the inner product

<< V 1 , V 2 >>= tr [V 1 ] T [M][V 2 ]
. Therefore, [Z] can be built in the following form:

[Z] = [V][a] + [V ⊥ ][b] (C.57)
of degrees of freedom in q j domain, and N the number of products of functions required to converge the PGD solution of the problem. Using NPM, we choose to enrich some (or all) of the PGD reduced coordinates solution using the same algorithm depicted in equation (C.58).

The results for a reduced basis of a domain [q j ] results in the following:

                       [U] = [G(β)][σ] [A] = [U] -[B]([B] T [U]) [D] = [Q] T j [M][A] + [A] T [M][Q] j /2 [Z] = [A] -[Q] j [D] [H] = [I n ] + s 2 [Z] T [M][Z] -1/2 [Q j ] = ([Q] j + s[Z]) [H] (C.62)
The algorithm illustrated in equation (C.62) represents an enrichment of one dimension in the PGD solution domain. The same algorithm can be designed to enrich as much dimensions as required, by generating the random basis [G(β)] as a higher dimensionality array. The number of hyper-parameters of the problem with increase accordingly.

In the optimization process, evaluating the quantities of interests O([Q], t) do not require any solution of the PDE problem, as the PGD defines a surrogate model replacing the PDE with a product of functions defined in every dimension of the domain. For instance, the evaluation of the quantities of interest for any solution requires only the knowledge of the solution y defined by: y(q 1 ; q

2 ; • • • ; q D ) = i=N i=1 j∈D d Q ij (q j ) j∈Dp Q ij (q j ), (C.63)
where D d are the deterministic, non-enriched basis, while D p are the probabilistic basis, enriched using NPM. Therefore, one solution and thus one evaluation of the quantities of interests, is computed only by using the product and sum of matrices.

The automated tape placement process, simulation and measurements

The automated tape placement process (ATP) is a one step, out-of-autoclave, composite manufacturing process [START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Levy | Simulation and optimization of the thermoplastic automated tape placement (atp) process[END_REF]. It aims a continuous deposit of prepreg tapes, while heating the deposit region and compressing the incoming tape to achieve consolidation. The process has the potential to achieve in a unique step the manufacturing of composite parts in the final shape.

Experimental process

The studied process is the deposition of prepregs on a heated cylindrical shape substrate, while heating with two external infrared heat source. The process is illustrated in figure C.21. The system is also equipped with two infrared camera to measure the temperature fields during deposition. The measurement is performed on 6 selected points as marked in figure C.22. The used material is prepreg tapes of carbon fiber impregnated with a thermoplastic matrix.

The measurements are performed at a constant time step ∆t = 0.5s. After measurements, the data is formatted and filtered using a moving least squares approach [33,[START_REF] Ghnatios | Knee joint injury risk assessment bymeans of experimental measurements and proper generalized decomposition[END_REF]

PGD simulation

In this part we explain the modeling and simulation performed in the context of the ATP process. The modeling is performed in the cylindrical (r, θ, z) coordinate system. When considering the reference as the IR 1 for example, one can work in a steady state formulation considering the convection-diffusion equation. In cylindrical coordinates, the governing equation becomes:

- ρC P r w ∂U ∂θ + K r 1 r ∂ ∂r r ∂U ∂r + K θ 1 r 2 ∂ 2 U ∂θ 2 + K z ∂ 2 U ∂z 2 = 0 (C.64)
with U being the temperature, ρ being the density, C p the thermal capacity at constant pressure, w the angular velocity, (K r , K θ , K z ) the diagonal components of the conductivity tensor. One may note that the deposition in the studied process orients the fibers in the tangential direction only, therefore the out-of-diagonal components of the conductivity tensor are all set to zero.

U(r, θ, z, T in , T b , w, Q 1 , Q 2 ) = i=N i=1 R i (r)Z i (θ, z)T i (T in )U i (T b )Q i (Q 1 )P i (Q 2 ) (C.66)
Readers unfamiliar with the PGD solution in separated form are refereed to the following references and their references therein [45,[START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]. Moreover, a thermal contact resistance is introduced between deposited ply, as described in [START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Ghnatios | Simulation avancée des problèmes thermiques rencontrés lors de la mise en forme des composites[END_REF], with a thermal conductance h = 6000W/m 2 .K. The involved problem parameters are depicted in table C.

NPM method applied for the PGD simulation of the automated tape placement process

Once the experimental and simulation results available, the NPM for PGD algorithm can be used to update and correct the simulated results with the experimental data. The experimental data are considered ground truth in this work, despite the large variability shown in some measurements. The cost function to optimize in this problem is given by: J = 0.9

Conclusion

In this work, we illustrate the possibility of using the Non-parametric Probabilistic Method (NPM) in the Proper Generalized Decomposition (PGD) framework. The formulation for the PGD framework shows the possibility to enrich the solution and correct it using the NPM, without the need to perform any extra PDE solutions. The selected application was the ATP, which deposits unidirectional fiber reinforced plastics 194 on a cylindrical shape using a continuous deposition, and infrared heating. The process is modeled using a novel approach and the simulation results are compared to the experimental ones, before and after enrichment. Before enrichment, the results are in good agreement with the experimental ones. The enrichment process induced a large improvement the results and illustrates the experimental variability as reflected on the presented results.

195

Résumé : La modélisation de la mise en forme des matériaux composites et métalliques est confrontée à de nombreux verrous scientifiques malgré les avancées récentes en matière de modélisation mécanique, analyse numérique, stratégies de discrétisation et capacité de calcul. En effet, la modélisation bien qu'elle soit avancée, est toujours incapable de reproduire tous les phénomènes physiques. Ceci est évidement dû à l'ignorance humaine de plusieurs phénomènes naturelles existants, mais pas découvert à présent. De plus, la capacité de calcul est en perpétuelle augmentation depuis des années. C'est ainsi que la modélisation par la voie des données prend de plus en plus de popularité parmi les chercheurs, promettant ainsi une modélisation "générique", pouvant potentiellement dépasser l'erreur humaine existante dans les modèles classiques.

Cepedant, il s'avère que la modélisation par les voie des données et l'intelligence artificielle n'est actuellement pas la meilleur façon de modélisation. En effet, en mécanique, plusieurs complications existent en partant sur une voie de modélisation par données. D'abord, les données sont chères. Chaque expérience doit-être contrôlée énormément pour éviter des artefacts expérimentaux qui influenceront les modèles générés. De plus, les moèles classiques en mécanique sont pertinents et reproduisent la réalité de façon fidèle dans la plupart des cas. Enfin, actuellement nous avons pas de méthodes de certification des modèles générés par des données et l'intelligence artificielle. Ainsi, ces modèles ne pourront pas être utilisées dans l'aéronautique par exemple.

En vue de ces complications, la voie la plus plausible pour faire des avancées technologique d'un point de vue modélisation sera l'utilisation des modéles hybrides, aussi connus sous le non de jumeaux digitaux. Utilisant les jumeaux digitaux, les résultats des modèles classiques sont conservés mais corrigés ou améliorés par des modèles à base de données. Ceci bien en considérant que les corrections auront une amplitude réduite par rapport aux solutions obtenus par la mécanique classique. De plus, il n'y aura pas nécessités d'avoir un nombre extensive d'expériences, comme la correction aura normalement une amplitude réduite. Finalement, la certification sera faite uniquement sur des modèles classiques et la correction ne sera qu'un surplus améliorant les résultats.

C'est ainsi que ce travail est divisée en trois parties majeurs:

• L'amélioration des méthodes de simulation avancée, notament par réduction de modèle, créant aussi des possibilités d'extension et d'innovation dans de nouveaux domaines, non traités actuellement par ces méthodes.

• Le développement des simulations avancées et des abaques numériques par des méthodes de réduction de modèles. Ces applications sétendent sur des modèles fluides, déformation mécaniques, thermiques, éctromagnétiques et couplages des différentes physiques.

• Le développement des méthodes permettant de construire des modèles à partir des données, mais aussi couplées aux méthodes de simulations existantes pour construire des jumeaux digitaux. Dans cette partie on valorise aussi les méthodes de reduction de modèles par la construction de jumeaux digitaux en temps réel.

Dans la première partie du travail, une étude de développements des méthodes de réduction de modèles est presentée dans plusieurs cas. D'abord, une extension de la méthode de stabilisation de la PGD ou Proper Generalized Decomposition pour satisfaire les condition LBB dans la simulation des écoulements des fluides est expliquée. De plus, la simulation par PGD est étendu pour adresser le problème de la non-intrusivité de la PGD, par une solution garantissant la création des abaques numériques par des méthodes non intrusives. Finalement, dans cette partie une présentation de l'accélération de la PGD par des méthodes d'intelligneces artificielles est introduites. Cette dernière methode possède aussi le potentiel d'être non intrusive.

Dans la deuxième partie de ce mémoire, une étude sur les différentes applications récentes de la simulation par des modèles réduits est présentée. Plusieurs applications industrielles sont conçues et implémentées à travers l'utilisation des modèles réduits. Quelques unes sont présentées en détail dans la deuxième partie de ce manuscrit. Plusieurs autres applications sont publiés dans les dernières années et sont cité dans ce manuscrit.

La dernière partie de ce travail porte sur la modélisation par la voie des données, et les nouvelles méthodes concu dans ce domaine. Ces méthodes sont validées expérimentalement et comparées a des méthodes de l'état de l'art et l'intérêt y est tiré pour les cas traités. De plus, dans un deuxième temps, les méthodes de sciences des données et d'intelligence artificielles sont couplés aux méthodes de réduction de modèles pour en consevoir des jumeaux numériques en temps réels. Quelques applications sont détaillés alors que d'autres sont cités. Finalement, ce mémoire présente aussi une nouvelle méthode de création de jumeaux digitaux probabilistes couplant la PGD avec la NPM ou Non-parametric Probabilistic Method.

Abstract : Modeling composite and metallic forming processes faces many scientific challenges nowadays, despite the recent advances in terms of mechanical modeling accuracy, numerical analysis, discretization techniques and computing power. In fact, recent advanced models are still unable to reproduce all physical phenomenons. This is obviously the fruit of human ignorance of different physical phenomenons yet to be discovered. Considering the recent advancement in computing power, the data-driven modeling techniques are becoming more popular among researchers, with the promise to deliver generic modeling, over-passing classical ones and alleviating the human ignorance of underlying physical phenomenons.

However, it appears that the data-driven and artificial intelligence modeling may not be the best solution in the mechanical engineering field. In fact, multiple complications are faced when it comes to mechanical modeling through data. First of all, data is expensive in mechanical engineering. Each experiment should be extremely well controlled to avoid experimental artifacts, which would consequently influence the generated models. Moreover, classical existing models perform well in reproducing reality, in many situations. Finally, there is currently no certification mechanism for models solely created by the data-driven and artificial intelligence route. Therefore, these models can't be implemented in aeronautical transport applications for example.

In view of the aforementioned complications, the most plausible solution to leverage the possibilities of machine learning in mechanics appears currently to be the modeling by using hybrid models, also known as digital twins. Using digital twins, the results of the classical, well established models, are conserved, but corrected or updated by data-driven modeling techniques. Eventually one would expect that such corrections have a small amplitude with respect to the results of the well established models. Moreover, a reduced or limited number of experiments is enough to train the digital twins, as their aim is not to create a new model, but rather to update an existing one in a certain instance for example. Finally, certification of the model can be performed solely on the classical models, while the correction is nothing but a plausible amelioration.

Following the aforementioned route, this work in divided into 3 main parts:

• Improving existing advanced simulation methods, mainly model reduction techniques, creating therefore the possibilities of novels applications and innovations, not tackled by these techniques until now.

• Developing novel advanced simulation and numerical abacus using model reduction techniques in different applications. These applications range from fluid mechanics applications, deformation and elasticity, thermal modelings, electromagnetism and coupled multiphysics.

• Developing advanced methods facilitating the construction of data-driven models, and also coupling them to the existing simulations and constructing digital twins. In this part we leverage the capabilities of model reduction techniques to construct real-time digital twins.

In the first part of this work, novel developments performed in the model reduction techniques are presented, for different physics. First of all, the developments in the PGD formulation to perform a stabilization of the fluid simulation and therefore solve fluid flow problems while satisfying the LBB condition is presented. Moreover, the PGD intrusiveness issue is tackled, creating therefore a non-intrusive PGD formulation able to create computational vdemecums, using commercial software and a fully non-intrusive approach. Finally, in this part, wee accelerate the PGD through using artificial intelligence techniques. This last method has also the potential to be fully non-intrusive.

In the second part of this thesis, a study of different model reduction applications in recent industrial problems is shown. Multiple industrial models are created and simulated using model reduction techniques. Some applications are presented in this work, while others are cited. Novel applications are being addressed frequently as mentioned in the manuscript.

The last part of this work shows novel data-driven modeling techniques created mainly for low-data applications. These methods are validated experimentally and compared to classical artificial intelligence and data-driven modeling techniques, with their interest highlighted for the treated cases. On a later stage, these methods as well as classical artificial intelligence techniques are coupled with the simulation, mainly using model reduction algorithms, to create real-time digital twins applications. Some applications are detailed in this manuscript while others are cited. Finally, this document also presents a novel method to create real-time probabilistic digital twin applications by coupling the PGD to the NPM or the Non-parametric Probabilistic Method.

3. 1

 1 The microstructure interface defined by a sinusoidal form, separated by mapping the domain into a prismatic one . . . . . . . . . . . . . . . . . . . . . . 3.2 Simulation of a 3D electric engine using a 1D/1D/1D decomposition. The simulation includes 4 extra process parameters with a total number of degrees of freedom exceeding the Avogadro number 6.02214076 × 10 23 . . . . . . . . 3.3 Interpolation of RTM filling using permeability interpolation. To the top left classical interpolation of the pressure field, to the top right a classical interpolation of the filling. The bottom line includes to the left an interpolation of the pressure field using a smart morphing of space and to the bottom right an interpolation of the filling using a morphing of space . . . . . . . . . . . . . . 3.4 Linear variation of temperature on the top surface . . . . . . . . . . . . . . . 3.5 Summary of the neural network used to predict generic functions . . . . . . . 3.6 Plot of generic function with maximum occurring residual in the domain . . . 3.7 Solution of problem with largest mean error at different elevations . . . . . . . 3.8 Maximum error at each in-plane node for solution with the largest mean error . 4.1 Parametric solution of the pressure field as a function of different considered parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Solution of the pressure field in the parametric space for a given node in the physical domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Thermal field during the deposition of a single path . . . . . . . . . . . . . . 4.4 Maximum temperature achieved in 3D printing as a function of the process parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Construction to exctract a representative volume element, which is the middle element of the deposited 9 paths . . . . . . . . . . . . . . . . . . . . . . . . 4.6 A zoom showing the quality of the mesh . . . . . . . . . . . . . . . . . . . . 4.7 A 3D printed wall with a top ring . . . . . . . . . . . . . . . . . . . . . . . . 4.8 A 3D printed wall with a top ring, a side view . . . . . . . . . . . . . . . . . 4.9 A zoom showing the quality of the mesh . . . . . . . . . . . . . . . . . . . . 4.10 A zoom showing the quality of the mesh . . . . . . . . . . . . . . . . . . . . 4.11 The modeling of the ATP process . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Vertical section in the deposited layers during the deposition of the 3 rd and 5 th layers. Temperatures in o C. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Dumbbell representation of a fabric patch. . . . . . . . . . . . . . . . . . . . 4.14 Compressible fabrics when subjected to a pressure through the squeeze flow. . 4.15 Considered domain, boundary conditions and in-plane mesh (red nodes indicated the matrix region whereas the white regions are occupied by the conductor). Three cross sections, a, 1 and 2 are considered for visualizing the results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 E y on section 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 E y on section 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v 4.18 E y on section 1 for a domain thickness of 0.2m. . . . . . . . . . . . . . . . . 4.19 E y on section a for a domain thickness of 0.2m. . . . . . . . . . . . . . . . 4.20 Composite material cross section, obtained by using micro-tomography of a composite tape. In blue the fibers and in yellow the matrix . . . . . . . . . . 4.21 Propagation of the electric field in the direction normal to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Through thickness propagation of the electric field in the direction normal to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Propagation of the electric field in the direction normal to the fibers for ω 2 = 24.5GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 Through thickness propagation of the electric field in the direction normal to the fibers for ω 2 = 24.5GHz. . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Propagation of the electric field in the direction parallel to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Through thickness propagation of the electric field in the direction parallel to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Propagation of the electric field in the direction parallel to the fibers for ω 2 = 24.5GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Through thickness propagation of the electric field in the direction parallel to the fibers for ω 2 = 24.5GHz. . . . . . . . . . . . . . . . . . . . . . . . . . 4.29 The scanned 3D part with the red constituting the fibres, the dark blue is the matrix and the light blue is the air phase . . . . . . . . . . . . . . . . . . . . 4.30 Propagation of the electric field in the direction normal to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.31 Propagation of the electric field in the parallel normal to the fibers for ω 1 = 2.45GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.32 Heat generation term (in W/m 3 ) for ω 1 = 2.45GHz and a wave propagation normal to the fiber orientation. . . . . . . . . . . . . . . . . . . . . . . . . . 4.33 Heat generation term (in W/m 3 ) for ω 1 = 2.45GHz and a wave propagation parallel to the fiber orientation. . . . . . . . . . . . . . . . . . . . . . . . . . 4.34 Temperature field at t = 50s term for ω 1 = 2.45GHz and a wave propagation normal to the fiber orientation. . . . . . . . . . . . . . . . . . . . . . . . . . 4.35 Temperature field at t = 50s term for ω 1 = 2.45GHz and a wave propagation parallel to the fiber orientation. . . . . . . . . . . . . . . . . . . . . . . . . . B.30 Maximum temperature along x (m) on a parallel close trajectory to the one followed by the laser spot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.31 Cooling profile of points on the neighboring line. . . . . . . . . . . . . . . . . B.32 Zig-Zag printing sequence exacerbating neighboring effects. . . . . . . . . . . B.33 Maximum temperature when t < t cr . . . . . . . . . . . . . . . . . . . . . . . B.34 Local and global frames. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.35 Simulation flowchart for the calculation of both displacement components d 1 and d 2 contribution to the total one d. . . . . . . . . . . . . . . . . . . . . . B.36 Validation experiment geometry. . . . . . . . . . . . . . . . . . . . . . . . . . B.37 Part distortion (m) caused by the substrate heating induced by the deposited molten metal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.38 Part distortion (m) caused by the substrate and the deposit cooling. . . . . . . B.39 Final part distortion obtained by adding the configurations shown in Figs. B.37 and B.38 (in m). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.40 Calibration model geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . B.41 Calibration model distortion (m). . . . . . . . . . . . . . . . . . . . . . . . . B.42 Von Mises stress (Pa) distribution. . . . . . . . . . . . . . . . . . . . . . . . B.43 Stresses in the deposit of interest along the deposition direction. . . . . . . . . B.44 Deformation zoon-in. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.45 Equivalent 3D finite element mesh associated with the 2D-1D (in-plane-outof-plane) separated representation resolution. . . . . . . . . . . . . . . . . . . B.46 Case study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.47 Penalized domain (black) with Φ = 0 (left) and the Poisson equation solution in the extended domain (right). . . . . . . . . . . . . . . . . . . . . . . . . . B.48 Resulting concentric circular trajectories. . . . . . . . . . . . . . . . . . . . . B.49 Axial distortion (in m) when all trajectories start at the same position (left) and when that starting point is randomly choosen (right). . . . . . . . . . . . B.50 Part distortion (m) for (a) T=2226 K, (b) T=2905 K and (c) T=4000 K, when considering circular trajectories. . . . . . . . . . . . . . . . . . . . . . . B.51 Part distortion: printing trajectories (left) and the distortion(m) (right). . . . . B.52 Part distortion: printing trajectories (left) and the distortion (m) (right). . . . B.53 Cylinder distortion (m) printed from concentric circular trajectories, random starting positions and a deposition temperature of 2710 K. . . . . . . . . . . . C.1 The column studied for buckling to the left with its cross section to the right. . C.2 Goodness of fitting of y 1 as a function of ŷ1 for x i mesh using random values . C.3 Goodness of fitting of y 1 as a function of ŷ1 for x i mesh using all possible combination of selected pool values . . . . . . . . . . . . . . . . . . . . . . . C.4 Goodness of fitting of y 1 as a function of ŷ1 using a uniform finite element 1D mesh for x i values. Excellent fitting is found despite not finding the exact y 1 fucntion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.5 Goodness of fitting of y 2 as a function of ŷ2 for x i taking random values for i = 1, • • • , 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.6 Goodness of fitting of y 2 as a function of ŷ2 for x i taking all possible combinations of a pool of values, for i = 1, • • • , 4 . . . . . . . . . . . . . . . . . . . C.7 Goodness of fitting of y 2 as a function of ŷ2 for x i taking uniform meshes values170 C.8 Cost C for different starts for the Tabu-searched regression trying to identify y 2 using random values in x i . . . . . . . . . . . . . . . . . . . . . . . . . . . C.9 Goodness of fitting of y 3 as a function of ŷ3 for x i taking random values for i = 1, • • • , 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Developping Computational Vamedecums • Research projects manager (projects with the collaboration of General Electric, EADS, AEROLIA ...) • Head of ICARE project with AEROLIA: Deliver a real time simulation app for the pilot cockpit beams as a function of the external load. An optimization was also performed. • PhD thesis guidance • Teaching courses and lab for engineering students November 2012-June 2013 Project Manager : New Composite Processes Institut de recherche technologique Jules Verne, Nantes, France. • Development of existing out-of-autoclave composite forming processes. Delivered a real time simulation application for the tape placement process as a function of the process parameters • Use of the PGD to simulate thermoplastic consolidation. Achieved a real time app prediction consolidation and degradation. • Coordinate projects with collaboration of industrials and academics in order to bridge the gap in the composite domain • PhD thesis guidance November 2009-October 2012 Research Engineer : Simulation of composite forming processes Ecole Centrale Nantes -Technocampus EMC2 • Development of advanced numerical solutions • Modeling and simulation of composite structures while forming • Optimization of industrial processes • Contributions to Ecole Centrale research projects (ITER, AMORE, ...) as part of my PhD thesis • A patent of invention and a registered mark resulted of my work Ecole des mines de Paris -Sophia Antipolis • Subject: Multiphasic simulation of complex fluid flows in a regulator. • Year of thesis beginning: 2018 • Year of thesis defense: 2021 • Supervisor: Elie Hachem and Thierry Coupez • PhD thesis supervision contribution: 20%

Figure 3 . 1 :

 31 Figure 3.1: The microstructure interface defined by a sinusoidal form, separated by mapping the domain into a prismatic one

Figure 3 . 2 :

 32 Figure 3.2: Simulation of a 3D electric engine using a 1D/1D/1D decomposition. The simulation includes 4 extra process parameters with a total number of degrees of freedom exceeding the Avogadro number 6.02214076 × 10 23

Figure 3 . 3 :

 33 Figure 3.3: Interpolation of RTM filling using permeability interpolation. To the top left classical interpolation of the pressure field, to the top right a classical interpolation of the filling. The bottom line includes to the left an interpolation of the pressure field using a smart morphing of space and to the bottom right an interpolation of the filling using a morphing of space

Figure 3 . 4 :

 34 Figure 3.4: Linear variation of temperature on the top surface

Figure 3 . 5 :

 35 Figure 3.5: Summary of the neural network used to predict generic functions

Figure 3 . 6 :

 36 Figure 3.6: Plot of generic function with maximum occurring residual in the domain

1 TemperatureFigure 3 . 7 :Figure 3 . 8 :

 13738 Figure 3.7: Solution of problem with largest mean error at different elevations

Figure 4 . 1 :

 41 Figure 4.1: Parametric solution of the pressure field as a function of different considered parameters

Figure 4 . 2 :

 42 Figure 4.2: Solution of the pressure field in the parametric space for a given node in the physical domain

Figure 4 .

 4 Figure 4.1 shows the abacus explored in real time on Paraview, an open source parallel visualization software. Figure4.2 illustrates the possibility of exploring the parametric domain in real-time, for a given node in the physical space (x, y, z).

Figure 4 . 3 :

 43 Figure 4.3: Thermal field during the deposition of a single path

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Maximum temperature achieved in 3D printing as a function of the process parameters

Figure 4 . 6 :

 46 Figure 4.6: A zoom showing the quality of the mesh

Figure 4 . 9 :

 49 Figure 4.9: A zoom showing the quality of the mesh

Figure 4 . 11 :

 411 Figure 4.11: The modeling of the ATP process

  (a) Deposition of the third layer (b) Deposition of the fifth layer

Figure 4 . 12 :

 412 Figure 4.12: Vertical section in the deposited layers during the deposition of the 3 rd and 5 th layers. Temperatures in o C.

Figure 4 . 13 :

 413 Figure 4.13: Dumbbell representation of a fabric patch.

  Figures 4.15 shows a selected configuration of preforms, and figures 4.16, 4.17, 4.18 and 4.19 shows the electrical wave propagation for different sections and different wave frequency in the domain.

  (a) Compressible fabric inside fluid squeeze flow (b) A top view of the compressible fabric

Figure 4 . 14 :Figure 4 . 15 :

 414415 Figure 4.14: Compressible fabrics when subjected to a pressure through the squeeze flow.

Figure 4 . 16 :

 416 Figure 4.16: E y on section 1.

Figure 4 . 17 :

 417 Figure 4.17: E y on section 2.

Figure 4 . 18 :Figure 4 . 19 :

 418419 Figure 4.18: E y on section 1 for a domain thickness of 0.2m.

Figure 4 . 20 :

 420 Figure 4.20: Composite material cross section, obtained by using micro-tomography of a composite tape. In blue the fibers and in yellow the matrix

Figure 4 . 21 :

 421 Figure 4.21: Propagation of the electric field in the direction normal to the fibers for ω 1 = 2.45GHz.

Figure 4 . 22 :

 422 Figure 4.22: Through thickness propagation of the electric field in the direction normal to the fibers for ω 1 = 2.45GHz.

Figure 4 . 23 :

 423 Figure 4.23: Propagation of the electric field in the direction normal to the fibers for ω 2 = 24.5GHz.

Figure 4 . 24 :

 424 Figure 4.24: Through thickness propagation of the electric field in the direction normal to the fibers for ω 2 = 24.5GHz.

Figure 4 . 25 :

 425 Figure 4.25: Propagation of the electric field in the direction parallel to the fibers for ω 1 = 2.45GHz.

Figure 4 . 26 :

 426 Figure 4.26: Through thickness propagation of the electric field in the direction parallel to the fibers for ω 1 = 2.45GHz.

(4. 18 )

 18 Equation(4.18) clearly reflects the difference between the results expected from the set

Figure 4 . 27 :

 427 Figure 4.27: Propagation of the electric field in the direction parallel to the fibers for ω 2 = 24.5GHz.

Figure 4 . 28 :

 428 Figure 4.28: Through thickness propagation of the electric field in the direction parallel to the fibers for ω 2 = 24.5GHz.

Figure 4 . 29 : 59 Figure 4 . 30 :

 42959430 Figure 4.29: The scanned 3D part with the red constituting the fibres, the dark blue is the matrix and the light blue is the air phase

Figure 4 . 31 :

 431 Figure 4.31: Propagation of the electric field in the parallel normal to the fibers for ω 1 = 2.45GHz.

Figure 4 . 32 :

 432 Figure 4.32: Heat generation term (in W/m 3 ) for ω 1 = 2.45GHz and a wave propagation normal to the fiber orientation.

Figure 4 . 33 :

 433 Figure 4.33: Heat generation term (in W/m 3 ) for ω 1 = 2.45GHz and a wave propagation parallel to the fiber orientation.

Figure 4 . 34 :

 434 Figure 4.34: Temperature field at t = 50s term for ω 1 = 2.45GHz and a wave propagation normal to the fiber orientation.

Figure 4 . 35 :

 435 Figure 4.35: Temperature field at t = 50s term for ω 1 = 2.45GHz and a wave propagation parallel to the fiber orientation.

Figure 5 . 1 :

 51 Figure 5.1: The Tabu-search optimized regression algorithm.

Figure 5 . 2 :

 52 Figure 5.2: A digital twin algorithm correcting the modeling error with a neural network.

Figure A. 2 :Figure A. 6 :

 26 Figure A.2: The x component u in m/s of the velocity vector v in the simulated lid-driven cavity

Figure A. 12 :Figure A. 13 :Figure A. 14 : 2 -Figure A. 15 : 5 Figure A. 16 :Figure A. 17 :Figure A. 18 :

 1213142155161718 Figure A.12: The 3D velocity vectors generated by the squeeze flow of the considered domain

Figure A. 17

 17 Figure A.17 shows the x component u of the solution of the problem in the in-plane (x, y) domain at a section normal to z, at z = h/2, while figure A.18 shows the y component v of the solution at the same section. The z component w and the pressure P are shown in figures A.19 and A.20 respectively. The slight oscillations appearing in figure A.19 are generated by the small number of PGD product of functions n calculated. Such oscillation is decreased at the expense of computing more products of functions and therefore spending more time to compute the solution.

Figure A. 19 :Figure A. 20 :μFigure A. 21 :

 192021 Figure A.19: w component of the velocity field at z = h/2, µ = 1P a.s and U = 1m/s, using h = 0.1m. We remind the reader to mind the axis

Figure A. 22

 22 Figure A.22 illustrates a section in the 7D domain as z = h/2, µ = 0.55P a.s, U = 0.5m/s, h = 0.01m, t = 0.3m, and shows the x component u of the velocity field on this section, while figure A.23 illustrates the y component v of the same velocity field. Figures A.24 and A.25 illustrates respectively the z component w of the velocity field and the pressure P on the same section. Figure A.26 illustrates the 3D solution at the same combination of the parameters, while figure A.27 illustrates the solution in the parametric domain µ, h, t) at x = t/2, y = t/2, z = h/2 when U = 0.5m/s.

Figure A. 22 :Figure A. 23 :Figure A. 25 :μFigure A. 27 :

 22232527 Figure A.22: u component of the velocity field at z = h/2, when µ = 1P a.s, U = 1m/s, h = 0.01m and t = 0.1m

Figure A. 28 :

 28 Figure A.28: Domain containing a nonplanar internal boundary Because of the small domain thickness, i.e., H ≪ L, standard techniques usually fail in the discretization of the weak form. Those based on the use of separated representations

1 .

 1 Mapping Ω b into R b . The first mapping with r = (r, s) ∈ R b = (0, L) × (0, 1) reads x = r, y = s h(r), (A.47) or, equivalently, r = x, s = y h(x) .

  the differential operators involved in the weak form (A.40),∇ = ( ∂ ∂x , ∂ ∂y ) T as ∂r -∂• ∂s s h ′ (r) h(r) , with dx = det(J b )dr = h(r)dr, (A.51)allows rewriting the weak form (A.40). For this purpose, we define matrix B b

64 ) 4 Figure A. 29 :

 64429 Figure A.29: Interface in Ω defining the top and bottom domains Ω u and Ω b

Figure A. 31 :

 31 Figure A.31: Solution in the real domain Ω = Ω u ∪ Ω b with the interface location highlighted.

Figure A. 32 :

 32 Figure A.32: Solution in the reference domain R = R u ∪ R b .

Figure A. 33 :

 33 Figure A.33: Solution in the real domain Ω = Ω u ∪ Ω b with the interface location highlighted.

Figure A. 34 :

 34 Figure A.34: 3D case study involving a non-planar interface.

Figure A. 35

 35 depicts on the reference domain R the solution on the domain cross section t = L z /2, while Fig. A.36 shows the solution at the corresponding cross section of the real domain Ω. Figure A.37 depicts the 3D plot of the solution on the real 3D domain Ω.

Figure A. 35 :

 35 Figure A.35: Temperature u(r, s, L z /2).

Figure A. 36 :

 36 Figure A.36: Temperature u(x, y, L z /2).

Figure A. 37 :

 37 Figure A.37: 3D solution representation. The domain dimensions are 1m × 10 -3 m × 1m.

h 2 (Figure A. 38 :

 238 Figure A.38: The geometrical transformation depicted in eEq. (A.69)

  Figure A.39 depicts the PGD solution of the problem while Fig. A.40 illustrates the one related to the solution of the 100 1D problems, both compared in Fig. A.41 that reveals a maximum relative error of 0.7%. The maximum error appears at the bottom left and top right regions of the domain, places where the mesh is excessively deformed after transformation.

Figure A. 39 :

 39 Figure A.39: First numerical example: Parametric solution u(x, X) computed by using the separated representation constructor operating on mapped domain R.

Figure A. 40 :

 40 Figure A.40: First numerical example: Solution for the different inclusion positions, u(x; X i ), i = 1, . . . , 100.

Figure A. 41 :

 41 Figure A.41: First numerical example: Relative error between between solutions depicted in Figs. A.39 and A.40.

93 .

 93 Figure A.42: Second numerical example: Parametric solution u(x, X) computed by using the separated representation constructor operating on mapped domain R.

Figure A. 43 :

 43 Figure A.43: Second numerical example: Solution for the different inclusion positions, u(x; X i ), i = 1, . . . , 100.

Figure A. 44 : 2 Figure A. 45 :

 44245 Figure A.44: Second numerical example: Relative error between between solutions depicted in Figs. A.39 and A.40.

Figure A. 46 :

 46 Figure A.46: The heat source in the real domain (X, s)

  Fig. A.50 lead to the one illustrated in Fig. A.51. In such case, a simulation using standard non-intrusive-PGD fails. RTM front propagation for a permeability k 2 = 5 × 10 -11

Figure A. 50 :

 50 Figure A.50: RTM front propagation solution for two different permeabilities k 1 and k 2 (Metric system units) at the same filling time

Figure A. 51 :

 51 Figure A.51: Classical front interpolation of two RTM fluid front at a permeability k = 10 -11 from two snapshots at k 1 = 5 × 10 -12 and k 2 = 5 × 10 -11

Figure A. 52 :enforced at position P 2 Figure A. 53 :Figure A. 54 :

 5225354 Figure A.52: An elastic rope supporting a load at a position x = X

  Figure A.55: Solution in the 2D real domain (x, X) and transformed domain (s, X)

  Figure A.56(a) shows the solution obtained using a classical interpolation in the x domain at X = L/2, while Fig. A.56(b) shows the result of the interpolation when performed in the s domain and then using an inverse transformation.

Figure A. 56 :Figure A. 57 :

 5657 Figure A.56: Solution at X = L/2 using classical interpolation and the proposed mapping technique

  obtained, capturing different temperature fields T i (x) for different locations of the heat source G i . The location of the different thermal sources is depicted in Fig. A.58. Snapshot solutions (two of them depicted in Fig. A.59

  Fig. A.60. Transformations between the two coordinate systems X = (x, y) and R = (r, s) for each triangle in Fig. A.60 read:

Figure A. 58 : 5 Figure A. 59 :Figure A. 60 :

 5855960 Figure A.58: Locations of different thermal sources G i in the different computed snapshots

Figure A. 61 :

 61 Figure A.61: Solution using classical interpolation and the one based on the proposed mapping

Figure A. 62 :Figure A. 63 :-Figure A. 64 :

 626364 Figure A.62: Interpolated solution in the reference domain

  Figure A.67 shows a classical interpolation of 2 solutions for different k x and k y , while Fig. A.68 shows the fronts R 1 (θ) and R 2 (θ). The interpolation of the front at every θ yields the result illustrated in Fig. A.69.

Figure A. 65 :-Figure A. 66 :

 6566 Figure A.65: Solution using flow front radius interpolation in the (r; θ) coordinates for k x = k y = 0.5k

(A. 84 )Figure A. 67 :Figure A. 68 :Figure A. 69 :Figure A. 70 :Figure A. 71 :

 846768697071 Figure A.67: Solution using classical interpolations for different permeability tensors

Figure A. 72 :

 72 Figure A.72: Different interpolated pressure fields at the same time step

5 ( 5 ( 5 ( 120 -

 555120 Figure A.73: Interpolated flow front positions for different combination of permeabilities

Figure A. 75 :

 75 Figure A.75: Reference solution obtained using finite element approach, calculated solution with the proposed mapping and relative errors between the pressure fields in the two solutions

Figure A. 76 :Figure A. 77 :O 3 Figure A. 78 :Figure A. 79 : 122 -Figure A. 80 :

 76773787912280 Figure A.76: The shell part considered for the SMC simulation. The coordinate s is a curvilinear coordinate used to map the 3D shell part into a 2D domain.

Figure B. 1 :

 1 Figure B.1: Model order reduction techniques are built on domain decomposition.

Figure B. 4 :

 4 Figure B.4: The complete model of the indenter

Figure

  Figure B.5:The "rods" approach to modeling the solid part

Figure B. 6 :

 6 Figure B.6: For a given combination of the 4 parameters, we illustrate the pressure field in the (ξ, γ, φ) domain

Figure B. 6

 6 Figure B.8: The velocity fields to a constant difference K µ • v in the 3D domain

  and B.14 respectively. Both solutions are almost identical, the relative error between the two solutions being shown in Figure B.15.

6 Figure B. 10 :Figure B. 11 :KFigure B. 12 :Figure B. 13 : 2 xFigure B. 14 :

 610111213214 Figure B.10: The pressure field in the thickness direction at x = 0 and y = 0 for a given combination of parameters

Figure B. 15 :

 15 Figure B.15: Relative error between the PGD solution shown in Figure B.13 and the finiteelement reference solution, shown in Figure B.14

Figure B. 16 :

 16 Figure B.16: Indentation force F as a function of the indentation depth w

Figure B. 17 :Figure B. 18 :

 1718 Figure B.17: Fitting of the indentation velocity U as a function of the indentation depth w

FigureFigure B. 19 :

 19 Figure B.19: Contribution of the solid and fluid force to the total force while displacing the cantilever base at a velocity of ż = 8µm/s

  Figure B.20: Contribution of the solid and fluid force to the total force while displacing the cantilever base at a velocity of ż = 80µm/s, the fluid reaction force exceeds that of the solid

2 )Figure B. 21 :

 221 Figure B.21: Identified poroelastic material permeability K for different ż

PressureFigure B. 22 :

 22 Figure B.22: Velocity magnitude at A(0, R/2, H/2) as a function of the pressure for different ż values (indicated ż values are in µm/s)

Figure B. 23 :

 23 Figure B.23: The Twin Cantilever experiment setup.

Figure B. 24 :

 24 Figure B.24: Heat transfer mechanisms involved in SLM

Figure B. 25 :

 25 Figure B.25: Temperature field at a given time t, with the source moving from the left boundary towards the right one

  Fig. B.25.

Figure B. 26 :

 26 Figure B.26: Maximum temperature at each position x (in m)

  Fig. B.27. Reheating Printing implies considering close trajectories as depicted in Fig. B.28. When the heat source comes close to a certain point during the next printing path, this point will be reheated, effect known as neighboring effect. As shown in Fig. B.27, the characteristic cooling time in order to reach a temperature verifying the condition T -T amb T amb 100 < 0.5%, (B.27)

  Fig. B.29.

Figure B. 28 :

 28 Figure B.28: Printing strategy

  Figure B.30: Maximum temperature along x (m) on a parallel close trajectory to the one followed by the laser spot.

Figure B. 31 :

 31 Figure B.31: Cooling profile of points on the neighboring line.

Figure B. 32 :

 32 Figure B.32: Zig-Zag printing sequence exacerbating neighboring effects.

Figure B. 33 :

 33 Figure B.33: Maximum temperature when t < t cr .

Figure B. 34 :

 34 Figure B.34: Local and global frames.

Figure B. 35 :

 35 Figure B.35: Simulation flowchart for the calculation of both displacement components d 1 and d 2 contribution to the total one d.

Figure B. 36 :

 36 Figure B.36: Validation experiment geometry.

  Fig. B.37.

Figure B. 37 :

 37 Figure B.37: Part distortion (m) caused by the substrate heating induced by the deposited molten metal.

Figure B. 38 :

 38 Figure B.38: Part distortion (m) caused by the substrate and the deposit cooling.

  [21].

Figure B. 39 :

 39 Figure B.39: Final part distortion obtained by adding the configurations shown in Figs. B.37 and B.38 (in m).

Figure B. 40 :

 40 Figure B.40: Calibration model geometry.

Fig. B. 43

 43 extracts the stresses in the deposit of interest.

  Figure B.41: Calibration model distortion (m).

Figure B. 42 :

 42 Figure B.42: Von Mises stress (Pa) distribution.

Figure B. 43 :

 43 Figure B.43: Stresses in the deposit of interest along the deposition direction.

Figure B. 44 :

 44 Figure B.44: Deformation zoon-in.

  Fig. B.47 shows the penalized domain and the Poisson equation solution. Hence, after calculating the solution, curves Φ = cte define the trajectories shown in Fig. B.48.

  Figure B.45: Equivalent 3D finite element mesh associated with the 2D-1D (in-plane-out-ofplane) separated representation resolution.

Figure B. 46 :

 46 Figure B.46: Case study.

Figure B. 47 :

 47 Figure B.47: Penalized domain (black) with Φ = 0 (left) and the Poisson equation solution in the extended domain (right).

Figure B. 48 :

 48 Figure B.48: Resulting concentric circular trajectories.

Figure B. 49 :

 49 Figure B.49: Axial distortion (in m) when all trajectories start at the same position (left) and when that starting point is randomly choosen (right).

Figure B. 50 :

 50 Figure B.50: Part distortion (m) for (a) T=2226 K, (b) T=2905 K and (c) T=4000 K, when considering circular trajectories.

  Fig. B.53 consisting of a cylinder of inner diameter 35 mm, thickness 10 mm and height 65 mm.

Figure B. 51 :

 51 Figure B.51: Part distortion: printing trajectories (left) and the distortion(m) (right).

Figure B. 52 :

 52 Figure B.52: Part distortion: printing trajectories (left) and the distortion (m) (right).

Figure B. 53 :

 53 Figure B.53: Cylinder distortion (m) printed from concentric circular trajectories, random starting positions and a deposition temperature of 2710 K.

Figure C. 1 :

 1 Figure C.1: The column studied for buckling to the left with its cross section to the right.

R 2 7 Figure C. 2 :

 272 Figure C.2: Goodness of fitting of y 1 as a function of ŷ1 for x i mesh using random values

2 2 and x 2 3 .

 3 Both extra terms are weighted many orders of magnitude lower than the other terms originally included in y 1 . The resulted regression yielded to an accurate fitting and prediction as shown in the goodness of fitting curve illustrated in figure C.2. Figure C.2 shows an excellent fitting with a minimal error on both the fitting and the training data.The two extra terms can be associated to numerical errors in the least square polynomial fitting.

R 2 = 1 ,Figure C. 3 : 1 ŷ1R 2 = 1 ,Figure C. 4 :

 2131214 Figure C.3: Goodness of fitting of y 1 as a function of ŷ1 for x i mesh using all possible combination of selected pool values

2 ŷ2R 2 = 1 , 8 Figure C. 5 : 4 Figure C. 6 :

 2218546 Figure C.5: Goodness of fitting of y 2 as a function of ŷ2 for x i taking random values for i = 1, • • • , 4

2 ŷ2R 2 = 1 ,Figure C. 7 :Figure C. 8 :

 22178 Figure C.7: Goodness of fitting of y 2 as a function of ŷ2 for x i taking uniform meshes values

3 ŷ3R 2 = 1 , 7 Figure C. 9 :Figure C. 10 :

 3217910 Figure C.9: Goodness of fitting of y 3 as a function of ŷ3 for x i taking random values for i = 1, • • • , 4

Figure C. 10

 10 Figure C.10 as well as C.11 both shows excellent fitting results for both training and evaluation data sets. However, the regression using mesh 2 was able to identify the monomial term x3 1 x 2 included in y 3 and F , but the one using mesh 3 never found that term even though figure C.11 show better results than C.10. This is explained by the non-independent nature of x i while using linear 1D finite element mesh to evaluate their values. The optimal regression using mesh 3 consisted of only 3 monomial terms, while the one found while using mesh 2 consisted of 21 monomial terms.

R 2 = 1 ,Figure C. 11 :

 2111 Figure C.11: Goodness of fitting of y 3 as a function of ŷ3 for x i taking uniform meshed values

Figure C. 13 :

 13 Figure C.13: The pressure (in Pa) of the 3D solution found using the PGD framework for a given combination of parameters

Figure C. 14 :

 14 Figure C.14: The experimental AFM nanoindentation reaction force in the spring F (t) as a function of the indentation depth w

Figure C. 15 : 16 PaFigure C. 16 :

 151616 Figure C.15: The identified paermeability as a function of the indentation depth w for different indentation velocities ż

Figure C. 18 : 2 n i=1 Θ 2

 1822 Figure C.18: An informed neural network to fit the ignorance of the model or the modeling error, Y being the experimentally measured force, F exp , while Y 1 , the output of the neural network consists of a fitting of the modeling error E

Figure C. 19 :

 19 Figure C.19: Neural network fitting performance on the three sets: training, validation and test sets

Figure

  Figure C.20: E fitted by the neural network

  [M][ÿ(t)] + [C][ ẏ(t)] + [K][y(t)] + [F (t)] = 0, (C.52)with t being the time domain, and the solution [y] can be projected on a reduced basis [V ] such as:[y(t)] = [V ][q(t)](C.53)

  . The experiment was repeated 4 times, depositing 6 plies in each experiment. The final filtered results from camera 1, for test one, one of the four measurements, are illustrated in figure C.23.

  (a) Picture of the ATP winding process

Figure C. 21 :

 21 Figure C.21: The ATP winding process with the two infrared heat sources marked as IR 1 and IR 2 , the incoming tape temperature is T in , the cylindrical deposition base has a radius R c = 54mm and a controlled temperature T b .

  (a) Thermal camera placed above the incoming tape, point 3 is an input for the simulation (b) Thermal camera placed below the incoming tape, point 5 is an input for the simulation

Figure C. 22 :

 22 Figure C.22: The thermal camera with the control points 1 and 2 from figure C.22(a) and 1 to 4 from figure C.22(b).

Figure C. 23 :

 23 Figure C.23: The measurement data performed by thermal camera 1 at control points 1 and 2.

K

  Multiplying equation (C.64) by r first and by a test function U * , then integrating by part, equation (C.64) will lead to the weak form of the problem: being the conductivity tensor, n the outward normal and Γ the boundary of the domain Ω. Equation (C.65) is therefore the governing equation to be solved, coupled with the boundary conditions depicted in figure C.24. A boundary condition of symmetry impose the continuity of the temperature between ply l at θ = 2π and ply l -1 at θ = 0. The tip of the very first layer is considered adiabatic. In the experimental results, T b and T in are constantly changing, therefore they are considered as extra parameters of the problem, varying inside a predefined interval: T in ∈ [25; 350] o C and T b ∈ [25; 250] o C. Moreover, the deposit angular velocity w as well as the thermal heat flux generated by the infrared IR 1 and IR 2 , Q 1 and Q 2 respectively, are also considered as extra parameters of the problem. The solution is therefore obtained in a separated form as:

2 .

 2 The simulation results are depicted in figure C.25 for the deposition of layers 3 and 5. The simulated results are compared to the experimental ones in figures C.26 and C.27. The experimental results show large variability. Only 4 experiments were available to compute the mean and the standard deviation. The results are shown for w = 0.0533 rad/s, Q 1 = 9250W/m 2 and Q 2 = 4500W/m 2 . T in and T b are experimental inputs.

Figure C. 24 :

 24 Figure C.24: The modeling of the ATP process

  Figure C.25: Vertical section in the deposited layers during the deposition of the 3 rd and 5 th layers. Temperatures in o C.

5 Figure C. 27 :Figure C. 28 : 3 .

 527283 Figure C.26: Comparison of the simulated temperature in the middle of the top layer, with experimental measurements, for the 3 rd and 5 th layer deposition. Temperatures in o C.

(

  E exp -E(U sim )) . 2 + 0.1 i=6 i=1(σ expσ(U sim )) .2 (C.[START_REF] Ducker | Measurement of forces in liquids using a force microscope[END_REF])With U sim = U (r exp , θ exp , z = Z max /2, T exp in , T exp b , w exp , Q exp 1 , Q exp 2 ),the evaluation of the simulation abacus at the points having the same input parametric values as the performed experiments. In equation (C.67), σ refers to the standard deviation. One may note that evaluating the cost function (C.67) would require multiple evaluation of the PDE solution, however no extra solution is required in the PGD framework. Only multiplications of vectors and matrices are involved. The enriched PGD solution with NPM for PGD slgorithm are illustrated in figures C.28 and C.29 for the temperature fields in the 3 rd and 5 th deposited layer respectively. The results are computed after performing m = 100 calculation of the enrichment basis. The PGD enrichment is limited to the 5 most energetic PGD product of vectors. Figures C.28 and C.29 show the results for NPM-PGD enhancement. The results shown an adaptation of the solution to overlap as much as possible with the experimental results, while translating the experimental results' variability into a simulation confidence interval.

Figure C. 29 : 5 .

 295 Figure C.29: Enrichement of the deposition of the fifth layer using the NPM for PGD algorithm, l = 5. Temperatures in C.

EDUCATION October 2012 Ph.D. Mechanical Engineer, under

  the Chair EADS-Ecole Centrale Nantes Title :Advanced Simulation of Thermal Problems Encountered in Composites Forming Processes Thesis supervisor : Prof. Francisco Chinesta Grade : Very honorable with the congratulations of the jury's president

	Professional address:
	Mechanics of materials Research Group
	Mechanical Engineering Department
	Notre Dame University -Louaizé
	P.O.Box: 72, Zouk Mikhael
	Zouk Mosbeh -Liban
	Email: cghnatios@ndu.edu.lb
	Tel.: +961 9 218 950 -Ext: 2015 / Fax: +961 9 218 771
	Cell.: +1 650 334 8941
	2 September 2009
	Masters in Mechanical Engineering-Applied Mechanics, Ecole Centrale Nantes
	Conception and Simulation of Complex Products and Systems
	Grade : Good

1 Identity GHNATIOS, Chady Born on Octobre 20, 1986 at Beirut -Lebanon, married. September 2009 Masters in Mechanical Engineering, Lebanese University Faculty of Engineering Last engineering year performed in Ecole Centrale Nantes, under a double degree program with Ecole Centrale Nantes Grade : Very Good • Accomplishment of more than 15 industrial research projects and granting agencies research projects as principal investigator or co-investigator, with some budgets exceeding a million EUR. • Research associate to ESI-Group, innovation department.

  Advani."3D Modeling of Squeeze Flow of Multiaxial Laminates". Journal of Non-Newtonian Fluid Mechanics, vol 234, 2016, pages 188-200. 19. C. Ghnatios, E. Abisset-Chavane, F. Chinesta, R. Keunings. "3D modelling of squeeze flow of unidirectional and fabric composite inserts". American Institute of Physics extended proceedings, 1769, 2016, pages 1 -10.

	20. J.-V. Aguado, D. Borzacchiello, C. Ghnatios, F. Lebel, R. Upadhyay, C. Binetruy, 7 PUBLICATIONS F. Chinesta. "A Simulation App based on reduced order modeling for manufacturing
	optimization of composite outlet guide vanes". Advanced Modeling and Simulation in Books Engineering Sciences, vol 4-1, 2017, pages 1-26.
	21. C. Ghnatios, C. Mathis, R. Simik, N.D. Spencer, F. Chinesta. "Modeling soft, per-1. F. Chinesta, N. Hascoet, A. Ammar, E. Cueto, A. Falco, C. Ghnatios, J.-L. Duval, meable matter with the Proper Generalized Decomposition (PGD) approach, and veri-A. Chambard, "Methodologies on Data reduction, Data-driven modeling, Model Order fication by means of nanoindentation". Soft Matter, 13, 2017, pages 4482-4493. DOI: Reduction & Twining for engineers in small bytes", ESI Group, San Jose (USA), 2020. 10.1039/c7sm00246g.
	search
	estimation of wear characteristics and cutting forces in compact core drilling of basalt "Composites Manufacturing Processes: Towards an Advanced Simulation". Revue de rock using PCD tool inserts", computers and Industrial Engineering, 136, 2019, pages Composites et Matériaux Avancés, 8, 2011, 1-15. 477-493.

Book Chapters

1. F. Chinesta, A.

Leygue

, E. Cueto, A. Ammar, C. Ghnatios, B. Bognet, "PGD Based Advanced Numerical Modeling of Materials and Processes", Centrale Nantes edition, Nantes (France), 2012. 2. F. Chinesta, E. Abisset-Chavanne, J.-V. Aguado, D. Borzacchiello, E. Lopez, A. Barasinski, D. Gonzalez, E. Cueto, C. Ghnatios, J.-L. Duval, "Big-data, machine learning, data-based models and data-driven simulations, avatars and internet of things, boarding the fourth industrial revolution". In "Modeling, successes and limites", Academy of technology-CNRS, Paris (France), pp. 131-146, 2017. 3. C. Ghnatios, B. Bognet, A. Barasinski, F. Chinesta, "Proper Generalized Decomposition Based Model Order Reduction Applied to Problems Defined in Degenerated Domains. In "ASME-ACIER book series", American society of mechanical engineers (USA), accepted. Refereed journal publications 1. C. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, P. Breitkopf, P. Villon. "Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: Application to pultrusion". Composites Part A, 42, 2011, 1169-1178. 2. B. Bognet, C. Ghnatios, F. Barasinski, F. Masson, F. Chinesta, A. Poitou, A. Leygue. 3. C. Ghnatios, F. Masson, A. Huerta, A. Leygue, A. Cueto, F. Chinesta. "Proper Generalized Decomposition Based Dynamic Data-Driven Control of Thermal Processes". Computer Methods in Applied Mechanics and Engineering, 213-216 (2012), 29-41. 4. C. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, A. Poitou."Optimizing Composites Forming Processes by Applying the Proper Generalized Decomposition". American Insitute of Physics proceedings, vol 1315, 2010, 197-202. 18. C. Ghnatios, E. Abisset-Chavane, C. Binetruy, F. Chinesta, S. 22. H. Tertrais, R. Ibanez, A. Barasinski, C. Ghnatios, F. Chinesta. "On the Proper Generalized Decomposition applied to microwave processes involving multilayered components". Mathematics and Computers in Simulation, vol 165, 2019, pages 347-363. 23. C. Ghnatios, E. Hachem. "A stabilized mixed formulation using the proper generalized decomposition for fluid problems". Computer Methods in Applied Mechanics and engineering, 346, 2019, pages 769-787. 24. M. Perez, A. Barasinski, B. Courtemanche, C. Ghnatios, F. Chinesta. "Sensitivity thermal analysis in the laser-assisted tape placement process". AIMS Materials Science, vol 5(6), 2018, pages 1053-1072. 25. C. Ghnatios, G. Asmar, E. Chakar, C. Bou Mosleh. "A reduced order model based manifold technique for automated defect judging with analytical validation". Comptes rendus de l'academie des sciences-mecanique, 347, 2019, pages 101-113. 26. C. Ghnatios, P. Simacek, F. Chinesta, S. Advani. "A non-local void dynamics modeling and simulation through the use of the Proper Generalized Decomposition". International journal of material forming, 2019, 1-14, doi: https://doi.org/10.1007/s12289-019-01490-7. 27. C. Ghnatios, E. Abisset-Chavanne, A. Ammar, E. Cueto, J.-L. Duval, F. Chinesta. "Advanced separated spatial representations for hardly separable domains". Computer Methods in applied mechanics and engineering, 354, 2019, pages 802-819. 28. R.-M. Hage, I. Hage, C. Ghnatios, I.-S. Jawahir, R. Hamade. "Optimized tabu 29. M. Perez, D. Prono, C. Ghnatios, E. Abisset, J.L. Duval, F. Chinesta. "Advanced modeling and simulation of SMC processes". International journal of material forming, doi: 10.1007/s12289-019-01506-2. 30. C. Ghnatios, R. Hage, I. Hage. "An efficient Tabu-search optimized regression for data-driven modeling", Comptes Rendus -Mécanique, 347(11), November 2019, pages 806-816.

  .5.

	Terms	Values	Normalized error values
	Mean of Residuals	0.1181 • C	2.95 × 10 -4
	Maximum of Residuals 1.8754 • C	4.69 × 10 -3

Table 3 . 3 :

 33 Quality of predicted generic function

	(z)
	i
	S
	0	0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09	0.1
		z

Table B . 1 :

 B1 Values of the parameters used in the simulation

	Parameter	Value	S.I unit
	E	200	Gpa
	Cp l	970	J/Kg.K
	Cp s	460	J/Kg.K
	h	0.5	W/m 2 .K
	T l	1663	K
	T s	1371	K
	∆H	270.10 3	J/Kg
	P	25	W
	r 0	50.10 -6	m
	V	0.1	m/s
	a	2227.664	-
	b	-0.627271	-
	c	2.09554.10 -4	-
	d	22.35452	-

  2 1 + 0.26x 1 x 3 3 + 0.53x 2 x 3 + 26x 4 + -7.73 -37.5x 2 1 x 2 + 1.4x 2 × 10 -11 (C.14) Again, the optimized regression was able to find the terms included in y 1 , with an error having 11 orders of magnitude lower than the tested function. The goodness of fitting is illustrated in figure C.3.Finally, the regression was also tested for a polynomial regression with the input variables x i taking uniform mesh values, parallel as what is performed in a unidirectional finite element mesh for example. The regression yields:

	ŷ1 = 125.68 + 0.13x 2 + 35.93x 4 3 + 1.49x 2

Table C . 2 :

 C2 Process and material parameters.

(C.15) Eventually the regression shown in equation (C.15) is not comparable a priori to the function y 1 shown in equation (C.12). However, since the meshes are uniform without repetition, a

The domain mesh consists of 1000 nodes in the x domain and 1000 nodes in the y domain, an equivalent of 3 000 000 degrees of freedom in a 2D space. The solution of the problem is achieved on a normal portable PC within only 11.5 seconds for µ = 1P a.s and u 0 = 1m/s. The solution of the x and y components of the velocity field, u and v, are illustrated in figures A.9 and A.10 respectively, while the pressure P is illustrated in figure A.11. We can clearly see the parabolic distribution of the x component of the velocity field at the domain exit.

The squeeze flow in 3D

Another benchmark problem previously solved using the PGD algorithm along with the penalty formulation is the squeeze flow problem. The problem consists of a prismatic domain of height h = 1cm in the z direction, and the in-plane dimensions are L x = L y = 1m. The boundary conditions are given by:

The problem is solved using an in-plane-out-of-plane decomposition. Thus, the 3D problem is solved as a sequence of 2D × 1D solutions. The problem variables are written as: By defining the characteristic functions of both subdomains, χ u (r) and χ b (r),

and

and using the notations

The problem can therefore be solved in a separated representation as a sequence of problems in the r and s coordinate domains, by using the standard PGD procedure (see [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF] for more details).

The considered domain mapping should always conserve the separability of the mapped coordinates as well the Jacobian of the transformation. Such approach is trivially satisfied when working in stratified domains.

Numerical examples

This section addresses different case studies: 2D and 3D domains with non-separable interfaces, a parametric problem involving an inclusion and finally the problem of a thermal source moving along a one-dimensional domain. In what follows, we use a uniform mesh to discretize the domains.

2D domains with non separable interface

In this section the applied boundary conditions (in the sequel all units are assumed to be in the metric system) read: [START_REF] Dohi | A note on portfolio optimization with path-dependent utility[END_REF] In what follows, we analyze two kinds of non-separable boundaries.

Appendix B

Advanced simulation of materials and processes 1 3D modeling of soft biphasic materials

The mechanical behavior of biological materials is currently a very active research topic. The lack of effective replacements for some human tissues motivates a better understanding of these tissues as well as their replacements. For example, lubricating mechanisms of articular cartilage are not yet fully understood. Cartilage, as well as all other biphasic materials, can be modeled as a solid mesh impregnated by a fluid [START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF]. However, such biphasic materials are challenging to model [START_REF] Harris | Finite element prediction of cartilage contact stresses in normal human hips[END_REF][START_REF] Pawaskar | Modelling of fluid support inside articular cartilage during sliding[END_REF]. This is best demonstrated in a recent study that aimed to model the Atomic Force Microscopy (AFM) nanoindentation of hydrogels -also considered a mimic for cartilage -by means of finite element and analytical methods [START_REF] Sherstova | Nanoindentation and finite element modeling of chitosan-alignate multilayer coated hydrogels[END_REF][START_REF] Shi | Microstructure and friction properties of pva/pvp hydrogels for articular cartilage repair as function of polymerization degree and polymer concentration[END_REF]. Interstitial fluid pressurization is believed to be one of the main load-carrying phenomena and was shown to play a major role in frictional properties [29,[START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF][START_REF] Moore | An analytical model to predict interstitial lubrication of cartilage in migration contact areas[END_REF]. However, an effective modeling approach, or a methodological simulation-based approach for prediction of the properties of such materials remians elusive.

In this work, we propose an effective method for the modeling of poroelastic materials in liquid environments. The experimental part of the work consists of colloidal AFM nanoindentation of a thin layer of poly acrylamide (PAAm) brushes in an aqueous buffer solution, forming a biphasic soft material [START_REF] Li | Tuning surface and mechanical properties by amplified polyelectrolyte selfassembly: where grafting-from meets grafting to[END_REF]. The suggested model is based on decoupling the fluid and solid behaviors -a frequently used approach to describe poroviscoelastic materials, both synthetic and natural [14,29,[START_REF] Espinosa-Marzal | Understanding the role of viscous solvent confinement in the tribological behavior of polymer brushes: a bioinspired approach[END_REF][START_REF] Han | Time dependant nanomechanics of cartilage[END_REF]. The fluid behavior is modeled using the Darcy equation, which is commonly applied to compute the interstitial fluid pressure in a porous medium [35,[START_REF] Tlupova | Boundary integral solutions of coupled stokes and darcy flows[END_REF].

To compare the simulated 3D model to experimental AFM nanoindentation results, we have combined an innovative numerical technique called Proper Generalized Decomposition (PGD) with classical non-linear solving techniques such as a Newton algorithm. However, the Newton algorithm explicitly uses the sensitivities of the solution of the 3D Darcy solution. This requirement leads to multiple solutions of the differential equation at each time-step, which increases the computation time tremendously for usual simulation techniques. At this point, the PGD seems to be an appealing technique to be used in inverse parameters identification. PGD is a reduced order model technique based on the use of separated representations of space-time, space-coordinates or space-parameters [43,[START_REF] Cueto | Improving computational efficiency in lcm by using computational geometry and model reduction techniques[END_REF][START_REF] Ghnatios | Modélisation avancée des procédés thermiques rencontrés lors de la mise en forme des composites[END_REF]]. An illustration of space decomposition is shown in Figure B.1. Using the PGD we can consider all unknown partial differential equation parameters as extra coordinates of the problem [43]. Thus, the resulting Darcy model addressed in this work is defined in a space of seven dimensions. Using the PGD

Appendix C

Digital twins and data-driven techniques 1 A Tabu-search optimized regressions 1.1 Introduction

Artificial neural networks are indeed a good and easy tool to use for fitting, classification, and deriving new models out of data, however their accurate training requires a large number of data points. Regression on the other side requires a lower number of data points to efficiently find a good fitting [41]. Moreover, regressions have the possibility to illustrate the most prominent parameters of the buckling load, allowing therefore the user to plan an optimization scheme for the considered panel. However, a regression should give a form to the fitting model, unlike the artificial neural network. The form of the regression is often named by its "kernel". Different regression optimization schemes exist in the literature, either by using support vector regression to correct any bias in the resulting data [START_REF] Zhang | Bi-sparse optimization-based least squares regression[END_REF], or by clustering the data and generating a "neighborhood based" regression [42]. Other optimization techniques tend to identify regression parameters by using unconventional optimization algorithms like the bee colony algorithm using in [START_REF] Zhu | Optimized support vector regression algorithm-based modeling of ship dynamics[END_REF]. However, previously mentioned type of optimization applies a constant regression kernels. In general kernel optimization or parameters learning are treated after a selection of a learning kernel, although kernel optimization may lead to better results [START_REF] Strazar | Approximate multiple kernel learning with least-angle regression[END_REF]. Some other works aim to optimize the selected basis functions for regression, or the selected kernel parameters for the regression, using meta-heuristic algorithms in a large search space [START_REF] Hage | Optimized tabu search estimation of wear characteristics and cutting forces in compact core drilling of basalt rock using pcd tool inserts[END_REF][START_REF] Pacheco | A variable selection method based on tabu search for logistic regression models[END_REF].

Glover [START_REF] Glover | Future paths for integer programming and links to artificial intelligence[END_REF] was the first to introduce the term meta-heuristic which represents a class of promising algorithms for solving hard optimization problems. Aiming to find near optimal solutions meta-heuristic algorithms work at efficient and comprehensive exploration of the search space, using the governing mechanisms which imitate certain strategies taken from nature, social behavior, physical laws, etc. . . Some popular global optimization algorithms include: genetic algorithm (GA), simulated annealing (SA), particle swarm optimization (PSO), ant colony optimization (ACO), artificial bee colony (ABC), taboo search (TS), quantum annealing, artificial immune system (AIS), improved harmony search algorithm (IHSA), real coded genetic algorithm (RCGA) and many more... The past 20 years have witnessed the development of numerous meta-heuristic algorithms in several fields, including artificial intelligence, computational intelligence and soft computing [START_REF] Gholizadeh | Comprative study of three metaheuristics for optimum design of trusses[END_REF]. The structure of meta-heuristic algorithms starts with an initial set of independent variables and then evolving to obtain the global minimum/maximum of the objective (fitness) function. Meta-heuristics seeks on improving the and artificial intelligence neural network perform better than both other models.

Modeling of AFM nanoindentation of soft matter

Classicaly, the AFM nanoindentation of a poroelastic brush is performed using a rigid spherical probe [START_REF] Ghnatios | Datadriven generic modeling of poroviscoelastic materials[END_REF][START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF], rigidly attached to a cantilever, as illustrated in Figure C.12. The indentation process is controlled by the motion of the base of the cantilever, z(t), as shown in figure C.12. An optical sensor measures z(t), the displacement of the cantilever, right above the colloidal probe (the indenter). In this work, the cantilever is modeled as a linear spring with a constant resulting stiffness k s = 2.88N/m. An increasing force F (t) is applied on the spring to indent the specimen. This force is measured during the indentation.

Specimen Indenter w(t)

.12: The reference for z = 0 is taken at the point of contact of the indenter with the specimen, w(t) is the penetration depth of the probe into the specimen

The indentation domain is cylindrical of initial height H 0 , the indenter is a semi-spherical probe of radius ρ = 12µm. The equilibrium equation is written as:

F f luid being the fluid contribution to the reaction force on the indenter, overcoming the spring compression force, while F solid is the solid phase contribution to the reaction, and F spring is the measured force F (t) imposed on the specimen. Knowing k s , the spring force is written as:

where w(t) represents the indentation depth into the specimen. The reference of z(t) = 0 is considered as the position of the spring base at the initial contact of the indentation probe with the specimen.

Note that the force of the spring is time dependent during the penetration phase. In the following, for the sake of notation simplicity, the force in the spring is simply denoted F (t). Next, we model the individual contributions of the fluid and the solid in the biphasic material.

The fluid contribution to the reaction force

The pressure in the fluid can be found using Darcy's law, modeling the fluid behavior in a porous media. One can easily demonstrate by a simple projection that the integral of the pressure over a spherical probe is equal to the integral of the same applied pressure over a disk of the same radius R(t). If ρ is the radius of the indenter, the radius of the contact disk is written as [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]:

R(t) changes as a function of time since the probe penetration depth w(t) is changing. Using Darcy's law we may write:

v being the fluid velocity, µ the fluid viscosity and K the effective permeability, which also takes into account the material porosity. In this work, we are using a highly porous media, with a fluid percentage above 90%, and therefore the conservation of mass is written as [START_REF] Ghnatios | Modeling soft permeable matter with the proper generalized decomposition (pgd) approach, and verification by means of nanoindentation[END_REF]:

Replacing (C.21) into (C.22) and assuming an homogeneous permeability K and fluid viscosity µ in the domain, we find Darcy's pressure equation to solve:

along with the boundary conditions specified in section 2.3. One may note that the permeability of the domain K, the radius R, the height of the domain H and the compression velocity U are time dependent. The indentation velocity can be determined by deriving the indentation depth w(t) at a given generic time step n:

where ∆t is the time step. Moreover, one may define the actual domain height H(t) under the indenter at a generic time step by:

During indentation, the height H n under the indenter is changing, as well as the indenter radius R n of contact with the specimen and the indentation velocity U n . While the permeability K is an unknown parameter to identify, we will need to solve the problem for different values of H n , U n , R n the radius at a time step n and permeability K. Therefore, having a parametric solution of the problem P (x, y, z, H, U, R, K) becomes an appealing approach.

Solid contribution to the reaction force

In this section we define the solid contribution to the reaction force by considering the polymer network as an equivalent homogeneous solid. We define σ z as the normal stress in the equivalent solid specimen along the vertical direction z, ǫ z is the strain on the same direction.

In this work, we first indent the biphasic material at a very low indentation speed ( ż = 0.5µm/s ≈ 0 for instance), which is considered slow enough to neglect any contributions of the fluid or solid viscoelastic effects. Therefore, the measured force in the spring is considered as the elastic contribution of the solid component. Thus, we define the elastic contribution of the solid to the reaction force by:

In addition to the force defined in equation (C.26), there may be another contribution coming from the viscoelastic component for example, as well as other modeling inaccuracy that will be discussed in section 2.5.

Defining [a], [b] and [V ⊥ ] leads to the following construction algorithm [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]:

with [σ] being an n × n upper diagonal matrix of hyper-parameters ∈ R to be identified, [G(β)] is a second order centered random matrix ∈ R (N,n) built from two uniform random distributions as well as the hyper-parameter β ∈ [0; 1], as described in the appendix D in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. s ∈ [0; 1] is also a hyper-parameter to define. The problem hyper-parameters are found using the following minimization problem:

We define w ∈ [0; 1] as a weight to insure as much as possible a convexity in the cost function. E(•) defines the mean of a distribution, V ar(•) define the standard deviation of a given distribution. E exp and V ar exp are the experimental values of the target values of the quantities of interest O([V], t), t being the time domain. The quantities of interest are defined as a function of the solution of the PDE equations of the modeled problem:

with [p(t)] ∈ R[n, m] the solution of the PDE in the reduced order basis coordinate system, m being the length of the mesh in the time domain. Computing [p(t)] would require a direct solution of the PDE (C.52) in the reduced coordinate system, for every variation of the reduced basis [V]. This gets even more time consuming when the derivatives are to be computed numerically during the optimization process. Interested readers are referred to the [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF][START_REF] Soize | Probabilistic learning for modeling and quantifying modelform uncertainties in nonlinear computational mechanics[END_REF] and their references therein.

To evaluate stochastic quantities of interests of the generated stochastic solution y, derived from the expression of NPM, one may need to repeat the algorithm (C.58) m times, to construct several solutions for different selection values of [G(β)]. In fact, [G(β)] is the only random distribution, generated out of two uniform random distributions as explained in the appendix of [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and high-dimensional nonlinear models[END_REF]. Therefore, one evaluation of the cost function depicted in (C.59) requires m solutions of the problem depicted in the reduced basis [V]. Such approach can lead to prohibitive calculation times when solving complex PDEs in a large computation domain expressing multiple number of reduced coordinates.

NPM method for PGD solutions

The PGD solutions are computed in a separated form such as [45]: y(q 1 ; q 2 ; • • • ; q D ) = i=N i=1 j=D j=1 Q ij (q j ), (C.61) with q j the separated coordinates of the problem, and Q ij (q j ) the recurrent solutions in the domain q j . In a discrete form, Q ij becomes a matrix [Q] j ∈ R[N j , N ], with N j the number