Full text
Finite element based model order reduction for parametrized one-way coupled steady state linear thermo-mechanical problems Nirav Vasant Shaha, Michele Girfoglioa, Peregrina Quintelab,c, Gianluigi Rozzaa,∗, Alejandro Lengomind, Francesco Ballarine, Patricia Barralb,c aScuola Internazionale Superiore di Studi Avanzati (SISSA), via Bonomea 265, Trieste, 34136, Italy bInstituto Tecnolóxico de Matemática Industrial (ITMATI), currently integrated in CITMAGA, s/n, Campus Vida, Rúa de Constantino Candeira, Santiago de Compostela, 15705, Spain cDepartamento de Matemática Aplicada, Universidade de Santiago de Compostela, Santiago de Compostela, 15782, Spain dPrimary & By-Products Department, ArcelorMittal, Global R&D Asturias, P.O. Box 90, Avilés, 33400, Spain eDepartment of Mathematics and Physics, Catholic University of the Sacred Heart, via Musei 41, Brescia, 25121, Italy Abstract This contribution focuses on the development of Model Order Reduction (MOR) for one-way coupled steady state linear thermo-mechanical problems in a finite element setting. We apply Proper Orthogonal Decomposition (POD) for the computation of reduced basis space. On the other hand, for the evaluation of the modal coefficients, we use two different methodologies: the one based on the Galerkin projection (G) and the other one based on Artificial Neural Network (ANN). We aim to compare POD-G and POD-ANN in terms of relevant features including errors and computational efficiency. In this context, both physical and geometrical parametrization are considered. We also carry out a validation of the Full Order Model (FOM) based on customized benchmarks in order to provide a complete computational pipeline. The framework proposed is applied to a relevant industrial problem related to the investigation of thermo-mechanical phenomena arising in blast furnace hearth walls. Keywords: Thermo-mechanical problems, Finite element method, Geometric and physical parametrization, Proper orthogonal decomposition, Galerkin projection, Artificial neural network, Blast furnace. 1. Introduction Due to technological developments occurred in recent years, the high-fidelity numerical computations, based on the so-called Full Order Models (FOM) (e.g., finite element or finite volume methods), are required to be performed for many configurations. This puts the computational resources under considerable stress. In this context, Model Order Reduction (MOR) has been introduced as an efficient tool to accelerate the computations with “affordable” and “controllable” loss of accuracy. The faster computations obtained by MOR helped in many query contexts, e.g. quick transfer of computational results to industrial problems. The basic idea on which MOR is based is related to the fact that often the parametric dependence of the problem at hand has an intrinsic dimension much lower than the number of degrees of freedom associated to the governing FOM. The development of a MOR consists in two main steps. The first one is the so-called offline stage when a database of several high-fidelity solutions is collected by solving the FOM for different values of physical and/or geometrical parameters. Then all the solutions are combined and compressed to extract a set of basis functions that approximate the low-dimensional manifold on which the solution lies. The second step is the online stage when the information obtained in the offline stage is used to efficiently compute the solutions for new parameters instances. For a comprehensive review on MOR, the reader is referred to, e.g., [3, 8, 10, 11, 24, 28, 43, 50]. In this work, we address the development of a MOR framework in a Finite Element (FE) environment for one-way coupled steady state linear thermo-mechanical problems. In the literature different MOR techniques have been proposed ∗Corresponding author Email addresses: [email protected] (Nirav Vasant Shah), [email protected] (Michele Girfoglio), [email protected] (Peregrina Quintela), [email protected] (Gianluigi Rozza ), [email protected] (Alejandro Lengomin), [email protected] (Francesco Ballarin), [email protected] (Patricia Barral) Preprint submitted to Finite Elements in Analysis and Design June 13, 2022 arXiv:2111.08534v4 [math.NA] 10 Jun 2022
in the context of thermo-mechanical problems. Guérin et al. [23] developed Rational Craig-Hale methodology for the investigation of thermo-mechanical coupling effects in turbomachinery. Benner et al. [9] compared the performance of Proper Orthogonal Decomposition (POD), Balanced Truncation, Padé approach and iterative rational Krylov algorithm for the approximation of the transient thermal field concerning an optimal sensor placement problem for a thermo-elastic solid body model. Zhang et al. [60] introduced reduced order variational multiscale enrichment method and tested their approach on proper benchmark tests related to the thermo-mechanical loading applied to a 2D composite beam and a functionally graded composite beam. More recently, Hernández-Beccero et al. [26] used Krylov Modal Subspace method for thermo-mechanical models as applicable to machine tools. We highlight that all these works are focused on MOR for the efficient reconstruction of the time evolution of the thermo-mechanical field. Regarding steady state problems, such as those dealt with in this work, Hoang et al. [31] used a two-field reduced basis algorithm based on the greedy algorithm in a physical parametrization setting. However, here we consider not only physical parameters but also geometrical ones. Concerning the methodology adopted, we use POD for the construction of reduced basis space. On the other hand, for the computation of the reduced coefficients, we consider two different approaches: a standard Galerkin projection (G) and Artificial Neural Network (ANN). POD-G aims to generate a MOR by a projection of the governing equations onto the POD space. So, the reduced coefficients associated to the POD bases are obtained by solving a set of algebraic equations. On the other hand, POD-ANN belongs to the category of machine learning methods based on systems that learn from data. In this framework, the reduced coefficients are computed by the employment of a properly trained neural network. There is a broad range of strategies for the development of new deep learning architectures for non-intrusive MOR approaches, i.e. without the need to access to FOM implementation: see, e.g., [27, 49, 56, 16, 41, 18, 58, 2, 17, 38]. We compare POD-G and POD-ANN in terms of relevant computational features including errors and speed up. From this viewpoint, the present contribution draws inspiration by [27] in which it has been shown that POD-ANN performs better than POD-G both in terms of efficiency and accuracy for steady state problems in heat transfer and fluid dynamics. We underline that whilst in [27] nonlinear problems are addressed, here we deal with a linear modeling framework by obtaining significantly different results. Our approach is applied within an industrial framework related to the investigation of thermo-mechanical phenomena arising in blast furnace hearth walls. Neural networks can in theory represent any functional relationship between inputs and output. However, many applications are still unexplored, and we retain that the use of POD-ANN method in the context of thermomechanical problems of industrial interest could open the door towards the application of deep learning techniques to new multiphysics scenarios. The workflow of this paper is organized as follows. In Sec. 2 the physical problem related to the blast furnace hearth walls is introduced. The full order model in strong formulation is described in Sec. 3. The corresponding weak formulation is then derived in Sec. 4 and the finite element analysis is introduced in Sec. 5. Subsequently, in Sec. 6, the MOR approach is described, and results obtained are shown and discussed. Conclusions and perspectives are drawn in Sec. 7. Finally, we dedicate a wide Appendix to the FOM validation in order to show the functionality of the developed finite element code. 2. Physical problem Steelmaking is a very old process that has contributed to the development of technological societies since ancient times. The previous stage is the ironmaking process, which is performed inside a blast furnace, whose general layout is shown in Figure 1. It is a metallurgical reactor used to produce hot metal from iron ore. For further details the reader is referred, e.g., to [15, 20]. The blast furnace operates at a high temperature (up to 1500 °C). The associated thermal stresses significantly limit the overall blast furnace campaign period. In this context, thermo-mechanical modeling has been used extensively either to support experimental campaign or to design various components. Vázquez-Fernández et al. [57] simulated the stationary heat transfer in a trough of a blast furnace to ensure durability based on the location of critical isotherm. In Barral et al [5, 6], the transient behaviour of the temperature during a tapping and a full campaign cycle was analyzed. Numerical modeling of heat flows in the blast furnace hearth lining was used by Swartling M. et al. [55] for improving experimental assessment. Thermo-mechanical modeling of blast furnace hearth was also developed by Brulin et al. [14]: they used micro-macro approach with homogenization method for replacing bricks and mortars by an equivalent material. The blast furnace operates under different conditions, each of which is governed by a different mathematical model. Considering the objectives of the present work, the following simplifications are considered: 2
Figure 1: Blast furnace [Courtesy: ArcelorMittal, Spain]. •Taphole operation is not part of this study. The perforation action of the taphole and the important pressures in the draining of the hot metal and slag produce important mechanical stresses located in the area that requires a deeper analysis and that is out of the scope of this work. •Since the objective is to be able to calculate in real time the effects of wall design on blast furnace operation, we focus on the steady state operations. •We assume that the hearth is made up of a single homogeneous, elastic, and isotropic material with temperatureindependent material properties. •Heat transfer only by conduction within hearth walls will be considered. The temperature of the molten metal inside the hearth is assumed constant and known. Therefore, the fluid region will not be part of the problem. 3. Full order model In this section, we discuss the mathematical formulation corresponding to the physical problem described in Sec. 2. We present the thermo-mechanical model in cylindrical coordinates endowed with suitable boundary conditions in Sec. 3.1. Then, we introduce the axisymmetric hypothesis and derive the axisymmetric thermo-mechanical model in Sec. 3.2. (a) Section of hearth geometry [Courtesy : ArcelorMittal]. Γ𝑠 𝑓 Γ+ Γ𝑜𝑢𝑡 Γ− (b) Simplified domain Ωand its boundaries. Figure 2: Hearth geometry: three dimensional domains, the real one (a) and its simplification as well as its boundaries (b). 3
𝑟 𝑦 𝛾𝑠 𝑓 𝛾𝑜𝑢𝑡 𝛾𝑠 𝛾− 𝛾+ 𝜔 𝑟𝑚𝑎𝑥 𝑦𝑚𝑎𝑥 Fluid region Figure 3: Vertical section of the hearth geometry 𝜔and its boundaries. 3.1. Thermo-mechanical model in cylindrical coordinates We consider the three dimensional domain Ωas in Figure 2 (b), corresponding to the simplified hearth geometry. We represent the displacement vector field as −→ 𝑢and the temperature scalar field as 𝑇. Based on the simplifications listed in Sec. 2, the energy and momentum conservation equations [12, 21, 22] can be written in a cylindrical coordinates system, (𝑟, 𝑦, 𝜃), with (𝑟, 𝑦) ∈ 𝜔, the vertical cross section of Ωin 𝑟−𝑦plane, (see Figure 3), and 𝜃∈ [0,2𝜋), as, −1 𝑟 𝜕 𝜕𝑟 𝑟𝑘 𝜕𝑇 𝜕𝑟 −𝜕 𝜕𝑦 𝑘𝜕𝑇 𝜕𝑦 −1 𝑟 𝜕 𝜕𝜃 𝑘 𝑟 𝜕𝑇 𝜕𝜃 =𝑄 , in Ω,(1) 𝜕𝜎𝑟𝑟 𝜕𝑟 +𝜕𝜎𝑟 𝑦 𝜕𝑦 +1 𝑟 𝜕𝜎𝑟 𝜃 𝜕𝜃 +𝜎𝑟𝑟 −𝜎𝜃 𝜃 𝑟+𝑓0,𝑟 =0,in Ω, 𝜕𝜎𝑟 𝜃 𝜕𝑟 +𝜕𝜎𝜃 𝑦 𝜕𝑦 +1 𝑟 𝜕𝜎𝜃 𝜃 𝜕𝜃 +2𝜎𝑟 𝜃 𝑟+𝑓0, 𝜃 =0,in Ω, 𝜕𝜎𝑟 𝑦 𝜕𝑟 +1 𝑟 𝜕𝜎𝜃 𝑦 𝜕𝜃 +𝜕𝜎𝑦𝑦 𝜕𝑦 +𝜎𝑟 𝑦 𝑟+𝑓0,𝑦 =0,in Ω. (2) Here −→ 𝑓0is the body force density term, 𝑄is the heat source term and 𝑘 > 0is the thermal conductivity tensor. The thermo-mechanical stress tensor 𝝈is related to the strain tensor 𝜺through the Hooke’s law: 𝝈(−→ 𝑢)[𝑇]=𝜆Tr(𝜺(−→ 𝑢))𝑰+2𝜇𝜺(−→ 𝑢)−(2𝜇+3𝜆)𝛼(𝑇−𝑇0)𝑰,(3) where 𝑇0is the reference temperature, 𝜺is defined as, 𝜺(−→ 𝑢)=1 2(∇−→ 𝑢+∇−→ 𝑢𝑇),(4) 𝛼is the thermal expansion coefficient, and 𝜆and 𝜇are the Lamé parameters of the material. These latter can be expressed in terms of Young modulus, 𝐸, and Poisson ratio, 𝜈, as: 𝜇=𝐸 2(1+𝜈), 𝜆 =𝐸𝜈 (1−2𝜈)(1+𝜈).(5) If 𝑨denotes the matrix, 𝑨=𝐸 (1−2𝜈)(1+𝜈) 1−𝜈 𝜈 𝜈 000 𝜈1−𝜈 𝜈 000 𝜈 𝜈 1−𝜈000 0001−2𝜈 20 0 0 0 0 0 1−2𝜈 20 0 0 0 0 0 1−2𝜈 2 ,(6) 4
the stress-strain relationship (3) can be expressed in vector formulation as, {𝝈(−→ 𝑢)[𝑇]} =𝑨{𝜺(−→ 𝑢)}− (2𝜇+3𝜆)𝛼(𝑇−𝑇0){𝑰},(7) where the following column vectors have been considered: {𝝈}={𝜎𝑟𝑟 𝜎𝑦𝑦 𝜎𝜃 𝜃 𝜎𝑦 𝜃 𝜎𝑟 𝜃 𝜎𝑟 𝑦 }𝑇, {𝜺}={𝜀𝑟𝑟 𝜀𝑦𝑦 𝜀𝜃 𝜃 2𝜀𝑦 𝜃 2𝜀𝑟 𝜃 2𝜀𝑟 𝑦 }𝑇, {𝑰}={111000}𝑇. (8) We introduce the notations for the boundaries of the domain Ω(Figure 2), and its vertical cross section in 𝑟−𝑦 plane, 𝜔(Figure 3): Γ𝑜𝑢𝑡 =𝜕Ω∩ (𝑟≡𝑟𝑚𝑎𝑥)=𝛾𝑜𝑢𝑡 × [0,2𝜋), Γ+=𝜕Ω∩ (𝑦≡𝑦𝑚𝑎𝑥 )=𝛾+× [0,2𝜋), Γ−=𝜕Ω∩ (𝑦≡0)=𝛾−× [0,2𝜋), Γ𝑠 𝑓 =𝜕Ω\(Γ𝑜𝑢𝑡 ∪Γ+∪Γ−)=𝛾𝑠 𝑓 × [0,2𝜋), 𝛾𝑠=𝜕𝜔 ∩ (𝑟≡0), (9) where 𝑟𝑚𝑎𝑥 ∈R+and 𝑦𝑚𝑎𝑥 ∈R+. Moreover, on the boundary, we define normal force 𝜎𝑛and tangential force −→ 𝜎𝑡as follows: 𝜎𝑛=(𝝈−→ 𝑛) ·−→ 𝑛 , −→ 𝜎𝑡=𝝈−→ 𝑛−𝜎𝑛−→ 𝑛 , (10) where −→ 𝑛is the unit normal vector that is directed outwards from Ω. •On the upper boundary, Γ+, the applied force, −→ 𝑔+, and the density of heat flux, 𝑞+, are known. Therefore, the following boundary conditions are considered: (−𝑘∇𝑇) ·−→ 𝑛=𝑞+,𝝈−→ 𝑛=−→ 𝑔+.(11) •On the bottom boundary, Γ−, the convection heat transfer occurs with the heat exchanger at temperature 𝑇−and heat transfer coefficient ℎ𝑐,−. The normal displacement is null, and we denote with −→ 𝑔−the shear force. Therefore, it is verified: (−𝑘∇𝑇) ·−→ 𝑛=ℎ𝑐,−(𝑇−𝑇−),−→ 𝑢·−→ 𝑛=0,−→ 𝜎𝑡=−→ 𝑔−.(12) •On the inner boundary, Γ𝑠 𝑓 , convection heat transfer with the fluid phase occurs and hydrostatic pressure 𝑔𝑠 𝑓 is acting. So, the following boundary conditions are considered: (−𝑘∇𝑇) ·−→ 𝑛=ℎ𝑐, 𝑓 (𝑇−𝑇𝑓),𝝈−→ 𝑛=−→ 𝑔𝑠 𝑓 ,(13) where 𝑇𝑓is the fluid temperature, assumed to be known and constant at the steady state, and ℎ𝑐, 𝑓 is the convective heat transfer coefficient. In addition, −→ 𝑔𝑠 𝑓 =−𝑔𝑠 𝑓 −→ 𝑛. •On the outer boundary, Γ𝑜𝑢𝑡 , a convective heat flux and known applied force −→ 𝑔𝑜𝑢𝑡 are assumed: (−𝑘∇𝑇) ·−→ 𝑛=ℎ𝑐,𝑜𝑢𝑡 (𝑇−𝑇𝑜𝑢𝑡 ),𝝈−→ 𝑛=−→ 𝑔𝑜𝑢𝑡 ,(14) ℎ𝑐,𝑜𝑢𝑡 being the convective heat transfer coefficient and 𝑇𝑜𝑢𝑡 the ambient temperature. 5
3.2. Axisymmetric thermo-mechanical model In the context of blast furnace application, the body force density term −→ 𝑓0, as well as surface forces, −→ 𝑔+,−→ 𝑔−,−→ 𝑔𝑠 𝑓 , −→ 𝑔𝑜𝑢𝑡 , have zero component in −→ 𝑒𝜃direction and they do not depend on 𝜃. Besides, the heat source term, 𝑄, the heat flux density, 𝑞+, the heat transfer coefficients, ℎ𝑐,−, ℎ𝑐, 𝑓 , ℎ𝑐,𝑜𝑢𝑡 , and temperatures 𝑇−, 𝑇𝑓, 𝑇𝑜𝑢𝑡 are assumed to be only dependent on (𝑟, 𝑦)coordinates. Therefore, a symmetry hypothesis is applicable that leads to significant computational savings. The associated axisymmetric model is reduced to consider conservation equations (1), (2) defined in 𝜔and the boundary conditions (11) - (14), where the Γboundaries are replaced by 𝛾such that Γ = (𝛾\𝛾𝑠) × [0,2𝜋)(see Figures 2 and 3, and definitions (9)). Moreover, the usual symmetry conditions on 𝛾𝑠are added: (−𝑘∇𝑇) ·−→ 𝑛=0,−→ 𝑢·−→ 𝑛=0,−→ 𝜎𝑡=−→ 0.(15) Therefore, the axisymmetric thermo-mechanical model can be summarized as: •Thermal model: −1 𝑟 𝜕 𝜕𝑟 𝑟𝑘 𝜕𝑇 𝜕𝑟 −𝜕 𝜕𝑦 𝑘𝜕𝑇 𝜕𝑦 =𝑄 , in 𝜔 , (16) with boundary conditions: on 𝛾+:−𝑘𝜕𝑇 𝜕𝑦 =𝑞+, on 𝛾−:𝑘𝜕𝑇 𝜕𝑦 =ℎ𝑐,−(𝑇−𝑇−), on 𝛾𝑠 𝑓 :−𝑘𝜕𝑇 𝜕𝑟 𝑛𝑟−𝑘𝜕𝑇 𝜕𝑦 𝑛𝑦=ℎ𝑐, 𝑓 (𝑇−𝑇𝑓), on 𝛾𝑜𝑢𝑡 :−𝑘𝜕𝑇 𝜕𝑟 =ℎ𝑐,𝑜𝑢𝑡 (𝑇−𝑇𝑜𝑢𝑡 ), on 𝛾𝑠:𝜕𝑇 𝜕𝑟 =0. (17) •Mechanical model: 𝜕𝜎𝑟𝑟 𝜕𝑟 +𝜕𝜎𝑟 𝑦 𝜕𝑦 +𝜎𝑟𝑟 −𝜎𝜃 𝜃 𝑟+𝑓0,𝑟 =0,in 𝜔 , 𝜕𝜎𝑟 𝑦 𝜕𝑟 +𝜕𝜎𝑦𝑦 𝜕𝑦 +𝜎𝑟 𝑦 𝑟+𝑓0,𝑦 =0,in 𝜔 , (18) In vector notation, axisymmetric stress-strain relationship can be expressed as, {𝝈(−→ 𝑢)[𝑇]} =𝑨{𝜺(−→ 𝑢)}− (2𝜇+3𝜆)𝛼(𝑇−𝑇0){𝑰}, 𝑨=𝐸 (1−2𝜈)(1+𝜈) 1−𝜈 𝜈 𝜈 0 𝜈1−𝜈 𝜈 0 𝜈 𝜈 1−𝜈0 0001−2𝜈 2 , {𝝈}={𝜎𝑟𝑟 𝜎𝑦𝑦 𝜎𝜃 𝜃 𝜎𝑟 𝑦 }𝑇, {𝜺}={𝜀𝑟𝑟 𝜀𝑦𝑦 𝜀𝜃 𝜃 2𝜀𝑟 𝑦 }𝑇, {𝑰}={1110}𝑇, (19) with the boundary conditions : 6
on 𝛾+:𝜎𝑟 𝑦 =𝑔+,𝑟 , 𝜎𝑦𝑦 =𝑔+,𝑦 , on 𝛾−:𝑢𝑦=0, 𝜎𝑟 𝑦 =−𝑔−,𝑟 , on 𝛾𝑠 𝑓 :𝜎𝑟𝑟 𝑛𝑟+𝜎𝑟 𝑦 𝑛𝑦=𝑔𝑠 𝑓 ,𝑟 , 𝜎𝑟 𝑦𝑛𝑟+𝜎𝑦𝑦𝑛𝑦=𝑔𝑠 𝑓 ,𝑦 , on 𝛾𝑜𝑢𝑡 :𝜎𝑟𝑟 =𝑔𝑜𝑢𝑡,𝑟 , 𝜎𝑟 𝑦 =𝑔𝑜𝑢𝑡,𝑦 , on 𝛾𝑠:𝑢𝑟=0, 𝜎𝑟 𝑦 =0. (20) 4. Weak formulation of axisymmetric thermo-mechanical model In this section, we derive the weak formulation related to the axisymmetric thermo-mechanical model (16)-(20). First, in Sec. 4.1, we introduce the relevant function spaces for temperature and displacement fields as well as for the model data, including boundary conditions, physical properties, and source terms, such that the problem is well defined. Next, the weak formulations for the thermal and mechanical problems are reported in Sec. 4.2 and 4.3, respectively. 4.1. Functional spaces For data of both thermal and mechanical problem, we introduce the weighted Sobolev spaces, 𝐿2 𝑟(𝜔), with norm || · ||𝐿2 𝑟(𝜔)as, 𝐿2 𝑟(𝜔)=𝑓:𝜔↦→ R,∫𝜔 𝑓2𝑟𝑑𝑟𝑑𝑦 < ∞, ||𝑓||2 𝐿2 𝑟(𝜔)=∫𝜔 𝑓2𝑟𝑑𝑟𝑑𝑦 . (21) Analogously, given 𝛾a subset of 𝜕𝜔, the boundary of 𝜔, 𝐿2 𝑟(𝛾)=𝑔:𝛾↦→ R,∫𝛾 𝑔2𝑟𝑑𝛾 < ∞.(22) Let 𝐿∞(𝜔)be the space 𝐿∞(𝜔)={𝑓:𝜔↦→ R,sup 𝜔|𝑓| ≤ 𝐶 , 𝐶 ≥0}, ||𝑓||𝐿∞(𝜔)=sup 𝜔|𝑓|.(23) Analogously, 𝐿∞(𝛾)is defined. For the temperature, we introduce the weighted Sobolev space, 𝐻1 𝑟(𝜔), with norm || · ||𝐻1 𝑟(𝜔)as, 𝐻1 𝑟(𝜔)=𝜓:𝜔↦→ R,∫𝜔 𝜓2+𝜕𝜓 𝜕𝑟 2 +𝜕𝜓 𝜕𝑦 2!𝑟𝑑𝑟𝑑𝑦 < ∞, ||𝜓||2 𝐻1 𝑟(𝜔)=∫𝜔 𝜓2+𝜕𝜓 𝜕𝑟 2 +𝜕𝜓 𝜕𝑦 2!𝑟𝑑𝑟𝑑𝑦 . (24) On the other hand, the following space Vfor the displacement is considered: V=(𝐻1 𝑟(𝜔) ∩ 𝐿2 1/𝑟(𝜔)) × 𝐻1 𝑟(𝜔).(25) It will be equipped with the inner product, <−→ 𝑢 , −→ 𝜙 >V=∫𝜔𝜙𝑟𝑢𝑟+𝜙𝑦𝑢𝑦+𝜕𝑢𝑟 𝜕𝑟 𝜕𝜙𝑟 𝜕𝑟 +𝜕𝑢𝑟 𝜕𝑦 𝜕𝜙𝑟 𝜕𝑦 +𝑢𝑟 𝑟 𝜙𝑟 𝑟+𝜕𝑢𝑦 𝜕𝑟 𝜕𝜙𝑦 𝜕𝑟 +𝜕𝑢𝑦 𝜕𝑦 𝜕𝜙𝑦 𝜕𝑦 +𝜕𝑢𝑟 𝜕𝑦 𝜕𝜙𝑦 𝜕𝑟 +𝜕𝑢𝑦 𝜕𝑟 𝜕𝜙𝑟 𝜕𝑦 𝑟𝑑𝑟𝑑𝑦 , (26) 7
and with the norm, ||−→ 𝜙||2 V=<−→ 𝜙 , −→ 𝜙 >V.(27) Its closed and convex subspace U, U={−→ 𝜙=𝜙𝑟𝜙𝑦∈V, 𝜙𝑦=0on 𝛾−, 𝜙𝑟=0on 𝛾𝑠},(28) will be the set of admissible displacements. The subspace Uis equipped with the same norm as spave Vi.e. ||−→ 𝜙||U=||−→ 𝜙||V,∀−→ 𝜙∈U. Finally, the function space for stress tensor is defined as, S={𝝈=[𝜎𝑖 𝑗 ]∈[𝐿2 𝑟(𝜔)]3×3, 𝜎𝑖 𝑗 =𝜎𝑗𝑖, 𝜎𝛼3=0, 𝛼 =1,2}.(29) For more details the reader is referred, e.g., to [29, 30, 35]. 4.2. Weak formulation for thermal model We assume the following hypotheses on the data: (TH1) The heat source term, 𝑄, is such that 𝑄∈𝐿2 𝑟(𝜔). (TH2) The convection temperatures, 𝑇𝑠 𝑓 ,𝑇−and 𝑇𝑜𝑢𝑡 , as well as the boundary heat flux 𝑞+are such that 𝑇𝑠 𝑓 ∈𝐿2 𝑟(𝛾𝑠 𝑓 ), 𝑇−∈𝐿2 𝑟(𝛾−), 𝑇𝑜𝑢𝑡 ∈𝐿2 𝑟(𝛾𝑜𝑢𝑡 ), 𝑞+∈𝐿2 𝑟(𝛾+). (TH3) The thermal conductivity 𝑘(𝑟, 𝑦)and the convective heat transfer coefficients, ℎ𝑐, 𝑓 (𝑟, 𝑦),ℎ𝑐,𝑜𝑢𝑡 (𝑟, 𝑦)and ℎ𝑐,−(𝑟, 𝑦)are such that 𝑘(𝑟, 𝑦) ∈ 𝐿∞(𝜔), 𝑘 (𝑟, 𝑦)> 𝑘0>0, ℎ𝑐, 𝑓 (𝑟, 𝑦) ∈ 𝐿∞(𝛾𝑠 𝑓 ), ℎ𝑐, 𝑓 (𝑟, 𝑦)> ℎ𝑐, 𝑓 ,0>0, ℎ𝑐,𝑜𝑢𝑡 (𝑟, 𝑦) ∈ 𝐿∞(𝛾𝑜𝑢𝑡 ), ℎ𝑐,𝑜𝑢𝑡 (𝑟, 𝑦)> ℎ𝑐,𝑜𝑢𝑡,0>0, ℎ𝑐,−(𝑟, 𝑦) ∈ 𝐿∞(𝛾−), ℎ𝑐,−(𝑟, 𝑦)> ℎ𝑐,−,0>0, where 𝑘0,ℎ𝑐, 𝑓 ,0,ℎ𝑐,𝑜𝑢𝑡,0and ℎ𝑐,−,0are suitable constants. In order to propose a weak formulation for the thermal model (16) - (17), in the following we assume sufficient regularity to perform the calculations. We multiply the energy equation (16) by 𝑟𝜓(𝑟, 𝑦), integrate over the domain 𝜔 with respect to (𝑟, 𝑦)variables, apply Gauss divergence theorem and use boundary conditions (17) to obtain: ∫𝜔 𝑟𝑘 𝜕𝑇 𝜕𝑦 𝜕𝜓 𝜕𝑦 +𝜕𝑇 𝜕𝑟 𝜕𝜓 𝜕𝑟 𝑑𝑟𝑑𝑦 +∫𝛾𝑠 𝑓 𝜓ℎ𝑐, 𝑓 𝑇𝑟𝑑𝛾 +∫𝛾𝑜𝑢𝑡 𝜓ℎ𝑐,𝑜𝑢𝑡𝑇𝑟𝑑𝛾+ ∫𝛾− 𝜓ℎ𝑐,−𝑇𝑟𝑑𝛾 =∫𝜔 𝜓𝑄𝑟𝑑𝑟𝑑𝑦 +∫𝛾𝑠 𝑓 𝜓ℎ𝑐, 𝑓 𝑇𝑓𝑟𝑑𝛾+ ∫𝛾𝑜𝑢𝑡 𝜓ℎ𝑐,𝑜𝑢𝑡𝑇𝑜𝑢𝑡𝑟𝑑𝛾 +∫𝛾− 𝜓ℎ𝑐,−𝑇−𝑟𝑑𝛾 −∫𝛾+ 𝜓𝑞+𝑟𝑑𝛾 . (30) It is to be noted that under assumptions (TH1)-(TH3) all integrals of the proposed weak formulation are well defined for all 𝑇, 𝜓 ∈𝐻1 𝑟(𝜔). The left hand side of equation (30) is bilinear and symmetric. So, we define in 𝐻1 𝑟(𝜔) × 𝐻1 𝑟(𝜔) the operator: 𝑎𝑇(𝑇, 𝜓)=∫𝜔 𝑟𝑘 𝜕𝑇 𝜕𝑦 𝜕𝜓 𝜕𝑦 +𝜕𝑇 𝜕𝑟 𝜕𝜓 𝜕𝑟 𝑑𝑟𝑑𝑦 +∫𝛾𝑠 𝑓 𝜓ℎ𝑐, 𝑓 𝑇𝑟𝑑𝛾 +∫𝛾𝑜𝑢𝑡 𝜓ℎ𝑐,𝑜𝑢𝑡𝑇𝑟𝑑𝛾+ ∫𝛾− 𝜓ℎ𝑐,−𝑇𝑟𝑑𝛾 . (31) 8
The right hand side of equation (30) is linear and the following operator defined on 𝐻1 𝑟(𝜔)is introduced: 𝑙𝑇(𝜓)=∫𝜔 𝜓𝑄𝑟𝑑𝑟𝑑𝑦 +∫𝛾𝑠 𝑓 𝜓ℎ𝑐, 𝑓 𝑇𝑓𝑟𝑑𝛾 +∫𝛾𝑜𝑢𝑡 𝜓ℎ𝑐,𝑜𝑢𝑡𝑇𝑜𝑢𝑡𝑟𝑑𝛾+ ∫𝛾− 𝜓ℎ𝑐,−𝑇−𝑟𝑑𝛾 −∫𝛾+ 𝜓𝑞+𝑟𝑑𝛾 . (32) Then, we can define the following problem: •Weak thermal model (WT) : Under the assumptions (TH1)-(TH3), find 𝑇∈𝐻1 𝑟(𝜔)such that, 𝑎𝑇(𝑇, 𝜓)=𝑙𝑇(𝜓),∀𝜓∈𝐻1 𝑟(𝜔).(33) By using Cauchy-Schwarz inequality, the trace operator properties, and Friedrich’s inequality [13, 39], it can be shown that, under the assumptions (TH1)-(TH3), 𝑎𝑇(𝑇, 𝜓)and 𝑙𝑇(𝜓)are continuous on 𝐻1 𝑟(𝜔) × 𝐻1 𝑟(𝜔)and 𝐻1 𝑟(𝜔), respectively, and 𝑎𝑇(𝜓, 𝜓)is coercive on 𝐻1 𝑟(𝜔) ×𝐻1 𝑟(𝜔). Hence the conditions of the Lax-Milgram theorem [13] are satisfied and accordingly the weak thermal model (𝑊𝑇)has a unique solution. 4.3. Weak formulation of the mechanical model We assume the following hypotheses on the data: (MH1) The body force density, −→ 𝑓0, is such that −→ 𝑓0∈ [𝐿2 𝑟(𝜔)]2. (MH2) The boundary forces −→ 𝑔+,−→ 𝑔𝑠 𝑓 ,−→ 𝑔𝑜𝑢𝑡 and −→ 𝑔−are such that −→ 𝑔+∈ [𝐿2 𝑟(𝛾+)]2,−→ 𝑔𝑠 𝑓 ∈ [𝐿2 𝑟(𝛾𝑠 𝑓 )]2,−→ 𝑔𝑜𝑢𝑡 ∈ [𝐿2 𝑟(𝛾𝑜𝑢𝑡 )]2,−→ 𝑔−∈ [𝐿2 𝑟(𝛾−)]2. (MH3) The Young’s modulus 𝐸(𝑟, 𝑦), the coefficient of thermal expansion 𝛼(𝑟, 𝑦), and the Poisson’s ratio 𝜈(𝑟, 𝑦)are such that 𝐸(𝑟, 𝑦) ∈ 𝐿∞(𝜔), 𝐸 > 𝐸0>0, 𝛼(𝑟, 𝑦) ∈ 𝐿∞(𝜔), 𝛼 > 𝛼0>0, 𝜈(𝑟, 𝑦) ∈ 𝐿∞(𝜔), 𝜈0< 𝜈 < 𝜈1, 𝜈0>0, where 𝐸0,𝛼0,𝜈0and 𝜈1are suitable constants. Analogously to what has been done for the thermal model, to propose a weak formulation of the mechanical model (18) - (20), in the following we assume sufficient regularity to perform the calculations. Given a function −→ 𝜙=(𝜙𝑟, 𝜙𝑦), we multiply the first equation of (18) by 𝑟𝜙𝑟(𝑟, 𝑦), the second one by 𝑟𝜙𝑦(𝑟, 𝑦), we sum both, integrate over 𝜔, apply Green formula and use equation (3) and (20) to obtain ∫𝜔 𝑨{𝜺(−→ 𝑢)} · {𝜺(−→ 𝜙)}𝑟𝑑𝑟𝑑𝑦 =∫𝜔(2𝜇+3𝜆)𝛼(𝑇−𝑇0){𝑰}·{𝜺(−→ 𝜙)}𝑟𝑑𝑟𝑑𝑦+ ∫𝜔𝜙𝑟𝑓0,𝑟 +𝜙𝑦𝑓0,𝑦 𝑟𝑑𝑟𝑑𝑦 +∫𝛾𝑠 𝑓 −→ 𝜙·−→ 𝑔𝑠 𝑓 𝑟𝑑𝛾 +∫𝛾𝑜𝑢𝑡 −→ 𝜙·−→ 𝑔𝑜𝑢𝑡𝑟𝑑𝛾+ ∫𝛾− −→ 𝜙·−→ 𝑔−𝑟𝑑𝛾 +∫𝛾+ −→ 𝜙·−→ 𝑔+𝑟𝑑𝛾 , ∀−→ 𝑢 , −→ 𝜙∈U, (34) where 𝑇is assumed to be the solution of the weak thermal model (𝑊𝑇). Notice that under assumptions (MH1)-(MH3), and since 𝑇∈𝐻1 𝑟(𝜔), all integrals in (34) are well defined for all −→ 𝑢 , −→ 𝜙∈U. The left hand side of equation (34), 𝑎𝑀(−→ 𝑢 , −→ 𝜙)=∫𝜔 𝑨{𝜺(−→ 𝑢)} · {𝜺(−→ 𝜙)}𝑟𝑑𝑟𝑑𝑦 , (35) 9
(a) Subdomains decomposition. (b) Close up of the mesh. Figure 6: Discretization of the domain ˆ𝜔. (a) Subdomains decomposition. (b) Close up of the mesh at bottom right. Figure 7: Improper domain decomposition: poor mesh quality at bottom right. (a) Subdomains decomposition. (b) Close up of the region affected by poor mesh quality. Figure 8: Improper domain decomposition: poor mesh quality under variation of the diameter 𝐷4. 16
where 𝐴is the area of the element, and 𝑙1,𝑙2and 𝑙3are the lengths of its three edges. The minimum value of 𝑞𝑒was 0.25, that is sufficiently far from zero. Notice that we use a coarser mesh with respect to the one used for the FOM benchmark tests in Appendix A. Such a choice is justified by the fact that the FOM solution is required to be solved at many parameters values, so using a fine mesh can be very costly and make prohibitive the collection of the high-fidelity database. The ranges of physical and geometrical parameters for training and testing are reported in Table 2. The sampling is carried out by using a Latin Hypercube Sampling (LHS) approach [33] which is a statistical method for generating near-random samples of parameter values from a multidimensional distribution. LHS divides the parameter space into equal partitions and samples parameters from each partition. In this manner, it is ensured that the patterns from entire parameter space are represented. The process has been repeated multiple times in order to ensure that random nature of samplings do not affect the final result. ANN has been trained by using the 70% of the total data provided by the full order model whilst the remaining 30% is used for the validation. Parameter Minimum value Maximum value 𝑡02.3 2.4 𝑡10.5 0.7 𝑡20.5 0.7 𝑡30.4 0.6 𝑡43.05 3.35 𝐷013.5 14.5 𝐷18.3 8.7 𝐷28.8 9.2 𝐷39.8 10.2 𝐷410.4 10.8 𝑘9.8 10.2 𝜇1.9e9 2.5e9 𝜆1.2e9 1.8e9 𝛼0.8e-6 1.2e-6 Table 2: Parameters ranges used for MOR training and testing. The accuracy of our MOR approach is quantified by the relative error defined as follows 𝜖𝑟𝑒𝑙,𝑋ℎ=||𝑋ℎ−𝑋𝑟𝑏 ℎ|| ||𝑋ℎ|| ,(69) where 𝑋ℎand 𝑋𝑟𝑏 ℎare the finite element solution and the corresponding reduced basis solution, respectively. We consider the projection error between the finite element solution 𝑋ℎand its projection on the reduced basis space 𝑋𝜋 ℎ, 𝜖𝑝𝑟𝑜 𝑗,𝑋ℎ=||𝑋ℎ−𝑋𝜋 ℎ|| ||𝑋ℎ|| ,(70) as benchmark for the relative error. || · || is the relevant norm (|| · ||𝐻1 𝑟,ℎ (𝜔)and || · ||Uℎ). 6.3.1. Thermal model We consider four numerical experiments that differ in terms of kind (physical and/or geometrical) and number of the parameters considered: •Numerical experiment (i): 1 physical parameter: Ξ = {𝑘}. •Numerical experiment (ii): 1 physical parameter and 3 geometric parameters: Ξ = {𝑘, 𝑡0, 𝐷2, 𝐷4}. •Numerical experiment (iii): 1 physical parameter and 6 geometric parameters: Ξ = {𝑘, 𝑡0, 𝑡2, 𝑡4, 𝐷0, 𝐷2, 𝐷4}. •Numerical experiment (iv): 1 physical parameter and all (10) geometric parameters: Ξ = {𝑘, 𝑡0, 𝑡1, 𝑡2, 𝑡3, 𝑡4, 𝐷0, 𝐷1, 𝐷2, 𝐷3, 𝐷4}. 17
Table 3 shows the number of samples provided by the full order model, 𝑛𝑡𝑟 , as well as the number of samples used for training and testing of ANN. Regarding the computation of POD space, for numerical experiment (i),50 FOM snapshots were considered while for the other ones 1000. The eigenvalues decay is shown in Figure 9. We see that the decay related to the numerical experiment (iv) is the slowest. This is due to the fact that in the numerical experiment (iv) we consider a larger number of parameters, so the system exhibits a greater complexity, and the modal content is more wide. Figure 10 shows the relative error (69) both for POD-ANN, related to different values 𝑛𝑡𝑟 and depth of hidden layers 𝐻, and POD-G. We also report the projection error (70). We observe that the performance of the POD-ANN method crucially depends on the values of 𝑛𝑡𝑟 and 𝐻. As expected, if we expand the training set and increase the depth of hidden layers, we obtain more accurate predictions when the number of parameters considered starts to get significative (numerical experiments (iii) and (iv)). Unlike [27], we observe that the POD-G method results to be in general more accurate than the POD-ANN method. This could be justified by considering that in the nonlinear framework, investigated in [27], the affine expansion could not be enforced and an Empirical Interpolation Method (EIM) [7] is used within the POD-G approach. Its implementation introduces interpolation error during the assembling of the reduced equations system by significantly affecting the accuracy of the POD-G method. Illustrative representations of the computed FOM and MOR are displayed in Figure 11 related to the numerical experiment (iv) for the parameters tuple Ξ = {2.365,0.6,0.6,0.5,3.2,14.10,8.50,9.2,9.9,10.6,10}. We use 4 POD basis. The POD-ANN solution was computed with 𝑛𝑡𝑟 =4500 and 𝐻=70. As we can see from Figure 11, both MOR approaches are able to provide a good reconstruction of the temperature field. We conclude by proving some information about the efficiency of our MOR approach. We report in Table 4 some estimations related to the offline time for all the numerical experiments carried out. We observe that the time taken by POD for numerical experiments (ii)-(iv) is much larger than the numerical experiment (i). This is fully justified by the fact that for the numerical experiments (ii)-(iv) we consider a larger number of snapshots (1000 instead of 50 as discussed above) for the computation of the reduced space. On the other hand, it should also be noted that the ANN training is faster for the numerical experiment (i) where only physical parameters are involved. This could be attributed to the fact that the introduction of geometric parameters increases the complexity of the input-output map that ANN is expected to learn. If on one hand the offline cost of the POD-G method is most composed of time taken by the computation of the snapshots from which the reduced space is extracted and the time taken by the computation of the POD modes, on the other hand the one related to the POD-ANN method is mainly associated to the computation of training data. So, when the parameter space is large (as for the numerical experiments (ii)-(iv)), the total offline cost of the POD-ANN method could be importantly greater than the one related to the POD-G method. We report in Table 5 the online time related to the POD-G and POD-ANN methods for all the numerical experiments carried out. As can be seen, the online time of POD-G method increases significantly in presence of geometric parameters by moving from 7𝑒−4s (numerical experiment (i)) to 1.3/1.5𝑒−2s (numerical experiments (ii)-(iv)). On the other hand, the time taken by POD-ANN online stage remains relatively constant for all the numerical experiments under investigation, around 5𝑒−4. So the computational efficiency of POD-ANN is much higher, of almost two order of magnitude, than POD-G when geometrical parametrization is considered. 𝑛𝑡𝑟 Training Testing Numerical experiment (i) 100 70 30 Numerical experiment (ii) 500 350 150 1500 1050 450 Numerical experiment (iii) 2000 1400 600 2500 1750 750 Numerical experiment (iv) 3500 2450 1050 4500 3150 1350 Table 3: Thermal model: number of total samples 𝑛𝑡𝑟 by FOM and number of samples used for training and testing of ANN. 18
Figure 9: Thermal model: plot of the eigenvalues {𝜃𝑖 𝑇}50 𝑖=1sorted in descending order for all the numerical experiments considered. (a) Numerical experiment (i). (b) Numerical experiment (ii). (c) Numerical experiment (iii). (d) Numerical experiment (iv). Figure 10: Thermal model: error analysis for POD-G and POD-ANN for all the numerical experiments carried out. 𝑡𝑃𝑂𝐷−𝐺 𝑜 𝑓 𝑓 𝑡𝑃𝑂𝐷−𝐴𝑁 𝑁 𝑜 𝑓 𝑓 𝑡𝑃𝑂𝐷 𝑡𝑡𝑟 𝑡𝐹𝑂𝑀 𝑡𝑝𝑟 𝑜 𝑗 Numerical experiment (i) ≈5≈1.2e1 8.4e-1 2.6e-1 8e-2 6.0e-4 Numerical experiment (ii) ≈9e1 ≈1.3e2 1.2e1 5.9e-1 7.5e-4 Numerical experiment (iii) ≈9e1 ≈3.7e2 1.2e1 7.4e1 7.1e-4 Numerical experiment (iv) ≈9e1 ≈5.0e2 1.2e1 4.5e1 7.3e-4 Table 4: Thermal model: time (in s) taken by (i) the entire offline stage (𝑡𝑜 𝑓 𝑓 ), (ii) the computation of the POD modes (𝑡𝑃𝑂𝐷), (iii) the training of ANN (𝑡𝑡𝑟 ), (iv) the computation of a FOM solution (𝑡𝐹𝑂𝑀 ) and (v) the projection of a FOM solution on the POD space (𝑡𝑝𝑟𝑜 𝑗 ). Concerning POD-ANN, we use 𝑛𝑡𝑟 =100, 𝐻 =65 for the numerical experiment (i), 𝑛𝑡𝑟 =500, 𝐻 =70 for the numerical experiment ii), 𝑛𝑡𝑟 =2500, 𝐻 =80 for the numerical experiment (iii) and 𝑛𝑡𝑟 =4500, 𝐻 =70 for the numerical experiment (iv). 19
Basis size POD-G POD-ANN Numerical experiment (i) 1 7.0e-4 4.9e-4 Numerical experiment (ii) 3 1.3e-2 4.8e-4 Numerical experiment (iii) 3 1.5e-2 4.9e-4 Numerical experiment (iv) 4 1.3e-2 5.1e-4 Table 5: Thermal model: online time (in s) for all the numerical experiments under investigation. Concerning POD-ANN, we use 𝑛𝑡𝑟 =100, 𝐻 =65 for the numerical experiment (i), 𝑛𝑡𝑟 =500, 𝐻 =70 for the numerical experiment ii), 𝑛𝑡𝑟 =2500, 𝐻 =80 for the numerical experiment (iii) and 𝑛𝑡𝑟 =4500, 𝐻 =70 for the numerical experiment (iv). (a) FOM solution (b) POD-G solution (c) POD-ANN solution Figure 11: Thermal model: comparison between the temperature field (in K) computed by the FOM and by the POD-G and POD-ANN methods related to the numerical experiment (iv) for Ξ = {2.365,0.6,0.6,0.5,3.2,14.10,8.50,9.2,9.9,10.6,10}. We consider 4 POD modes. For POD-ANN, we set 𝑛𝑡𝑟 =4500 and 𝐻=70. 6.3.2. Mechanical model We remark that for POD-ANN we refer to the (𝑊 𝑀)ℎmodel, whilst we consider (𝑊 𝑀1)ℎand (𝑊 𝑀2)ℎmodels for POD-G. As done for the thermal model, we consider four different numerical experiments having different kinds and numbers of parameters: •Numerical experiment (i): 4 physical parameters: Ξ = {𝑘, 𝜇, 𝜆, 𝛼}. •Numerical experiment (ii): 4 physical parameters and all 3 geometric parameters: Ξ = {𝑘, 𝜇, 𝜆, 𝛼, 𝑡0, 𝐷2, 𝐷4}. •Numerical experiment (iii): 4 physical parameters and 6 geometric parameters: Ξ = {𝑘, 𝜇, 𝜆, 𝛼, 𝑡0, 𝑡2, 𝑡4, 𝐷0, 𝐷2, 𝐷4}. •Numerical experiment (iv): 4 physical parameters and all (10) geometric parameters: Ξ = {𝑘, 𝜇, 𝜆, 𝛼, 𝑡0, 𝑡1, 𝑡2, 𝑡3, 𝑡4, 𝐷0, 𝐷1, 𝐷2, 𝐷3, 𝐷4}. Table 6 shows the total number of samples provided by the full order model, the number of samples used for training and the one used for testing of ANN. For all the numerical experiments, the POD space was computed by considering 1000 snapshots. The eigenvalue plot is shown in Figure 12. Like the thermal model, we observe that the numerical experiment (iv), characterized by the larger number of parameters, shows the lowest decay. On the other hand, as expected, among the different mechanical models we consider, the model (𝑊 𝑀)ℎexhibits the slowest eigenvalues decay including it both thermal and mechanical effects. 20
(a) Numerical experiment (i). (b) Numerical experiment (ii). (c) Numerical experiment (iii). (d) Numerical experiment (iv). Figure 12: Mechanical model: plot of the eigenvalues {𝜃𝑖 𝑀}50 𝑖=1sorted in descending order for all the numerical experiments considered. Figure 13 shows the relative error (69) both for POD-ANN and POD-G. The projection error (70) is also depicted. As observed for the thermal model, POD-G is able to provide more accurate results with respect to POD-ANN. Figure 14 shows the qualitative comparison between the computed FOM and MOR related to the numerical experiment (iv) for the parameters tuple Ξ = {2.365,0.6,0.6,0.5,3.2,14.10,8.50,9.2,9.9,10.6,10,2.08𝑒9,1.39𝑒9,1𝑒−6}. We use 7 POD basis. The POD-ANN solution was computed with 𝑛𝑡𝑟 =2500 and 𝐻=130. We could observe that both MOR approaches are able to provide a good reconstruction of the displacement field. In order to justify our choice to consider separately (𝑊 𝑀1)ℎand (𝑊 𝑀2)ℎin the Galerkin projection framework, we also highlight that, as shown in Figure 14, the scale difference between their displacements (derived from mechanical loads for the first and from the thermal ones for the second) is of one order of magnitude. Thus, the use of the model (𝑊 𝑀)ℎcould lead to a less accurate reconstruction of the displacement field. Finally, we briefly discuss the efficiency of our MOR approach. We report in Table 7 some estimations related to the offline time for all the numerical experiments carried out. We observe that, unlike the thermal model, the time taken by POD is comparable for all the numerical experiments. This is not surprising because in this case we consider the same number of snapshots for the computation of the reduced space for all the numerical experiments. Like the thermal model, the ANN training is faster for the numerical experiment (i), probably because of the minor complexity with respect to the other numerical experiments. Unlike the thermal model, the total offline cost of the POD-ANN method is comparable with the one related the POD-G method. This is because we train two models for the POD-G method that take a similar amount of time as training one model for the POD-ANN method. We report in Table 8 the online time related to the POD-G and POD-ANN methods for all the numerical experiments carried out. Like the thermal model, the online time of POD-G method increases significantly in presence of geometric parameters by moving from 8𝑒−4s (numerical experiment (i)) to 2.6/6.9𝑒−2s (numerical experiments (ii)-(iv)) for the model (𝑊 𝑀1)ℎand from 4.5𝑒−2s (numerical experiment (i)) to 1.9/2.6𝑒−1s (numerical experiments (ii)-(iv)) for the model (𝑊 𝑀2)ℎ. We could observe that the online time taken by the model (𝑊 𝑀2)ℎis significantly greater than that taken by the model (𝑊 𝑀1)ℎ. This is expected because for the model (𝑊 𝑀2)ℎa reduced basis approximation of temperature needs to be computed due to 21
the thermo-mechanical coupling. On the other hand, the POD-ANN, that does not need reduced basis approximation of temperature thanks to its non intrusive nature, is able to provide a higher computational efficiency. Moreover, like the thermal model, the POD-ANN online time remains relatively constant for all the numerical experiments under investigations, around 5𝑒−4, by showing a low sensitivity at varying of the kind and number of parameters considered. (a) Numerical experiment (i). For POD-ANN 𝑛𝑡𝑟 =500 and 𝐻=60. (b) Numerical experiment (ii). For POD-ANN 𝑛𝑡𝑟 =500 and 𝐻=80. (c) Numerical experiment (iii). For POD-ANN 𝑛𝑡𝑟 =1000 and 𝐻=170. (d) Numerical experiment (iv). For POD-ANN 𝑛𝑡𝑟 =2500 and 𝐻=130. Figure 13: Mechanical model: error analysis for POD-G and POD-ANN for all the numerical experiments considered. 𝑛𝑡𝑟 Training Testing Numerical experiment (i) 500 350 150 Numerical experiment (ii) 500 350 150 Numerical experiment (iii) 1000 700 300 Numerical experiment (iv) 2500 1750 750 Table 6: Mechanical model: number of total samples 𝑛𝑡𝑟 by FOM and number of samples used for training and testing of ANN. Basis size POD-G (𝑊 𝑀1)ℎPOD-G (𝑊 𝑀2)ℎPOD-ANN (𝑊 𝑀)ℎ Numerical experiment (i) 1 8e-4 4.5e-2 6.7e-4 Numerical experiment (ii) 3 2.6e-2 1.9e-1 5.3e-4 Numerical experiment (iii) 4 5.4e-2 2.1e-1 5.2e-4 Numerical experiment (iv) 7 6.9e-2 2.6e-1 4.9e-4 Table 8: Mechanical model: online time (in s) for all the numerical experiments under investigation. Concerning POD-ANN, we use 𝑛𝑡𝑟 =500, 𝐻 = 60 for the numerical experiment (i), 𝑛𝑡𝑟 =500, 𝐻 =80 for the numerical experiment (ii), 𝑛𝑡𝑟 =1000, 𝐻 =170 for the numerical experiment (iii) and 𝑛𝑡𝑟 =2500, 𝐻 =130 for the numerical experiment (iv). 22
𝑡𝑃𝑂𝐷−𝐺 𝑜 𝑓 𝑓 𝑡𝑃𝑂𝐷 𝑡𝐹𝑂𝑀 (𝑊 𝑀1)ℎ(𝑊 𝑀2)ℎ(𝑊 𝑀1)ℎ(𝑊 𝑀2)ℎ(𝑊 𝑀1)ℎ(𝑊 𝑀2)ℎ Numerical experiment (i) ≈4.2𝑒2≈6.2𝑒21.6e1 1.8e1 4e-1 6e-1 Numerical experiment (ii) ≈4.2𝑒2≈6.2𝑒21.6e1 1.7e1 Numerical experiment (iii) ≈4.2𝑒2≈6.2𝑒21.8e1 1.7e1 Numerical experiment (iv) ≈4.2𝑒2≈6.2𝑒21.7e1 1.7e1 (a) Problem (𝑊 𝑀1)ℎand Problem (𝑊 𝑀 2)ℎ 𝑡𝑃𝑂𝐷−𝐴𝑁 𝑁 𝑜 𝑓 𝑓 𝑡𝑃𝑂𝐷 𝑡𝑡𝑟 𝑡𝐹𝑂𝑀 𝑡𝑝𝑟𝑜 𝑗 (𝑊 𝑀)ℎ(𝑊 𝑀)ℎ(𝑊 𝑀)ℎ(𝑊 𝑀)ℎ(𝑊 𝑀)ℎ Numerical experiment (i) ≈1.1𝑒31.8e1 3.4 7e-1 8.1e-4 Numerical experiment (ii) ≈1.1𝑒31.8e1 3.7 9.2e-4 Numerical experiment (iii) ≈1.5𝑒31.6e1 2.9e1 1.0e-3 Numerical experiment (iv) ≈2.5𝑒31.8e1 7.6e1 9.1e-4 (b) Problem (𝑊 𝑀)ℎ Table 7: Mechanical model: time (in s) taken by (i) the entire offline stage (𝑡𝑜 𝑓 𝑓 ), (ii) the computation of the POD modes (𝑡𝑃𝑂𝐷), (iii) the training of ANN (𝑡𝑡𝑟 ), (iv) the computation of a FOM solution (𝑡𝐹𝑂𝑀 ) and (v) the projection of a FOM solution on the POD space (𝑡𝑝𝑟𝑜 𝑗 ). Concerning POD-ANN, we use 𝑛𝑡𝑟 =500, 𝐻 =60 for the numerical experiment (i), 𝑛𝑡𝑟 =500, 𝐻 =80 for the numerical experiment ii), 𝑛𝑡𝑟 =1000, 𝐻 =170 for the numerical experiment (iii) and 𝑛𝑡𝑟 =2500, 𝐻 =130 for the numerical experiment (iv). (a) FOM solution related to the problem (𝑊 𝑀1)ℎ (b) FOM solution related to the problem (𝑊 𝑀2)ℎ (c) FOM solution related to the problem (𝑊 𝑀)ℎ (d) POD-G solution related to the problem (𝑊 𝑀1)ℎ (e) POD-G solution related to the problem (𝑊 𝑀2)ℎ (f) POD-ANN solution related to the problem (𝑊 𝑀)ℎ Figure 14: Mechanical model: comparison between the displacement (in m) computed by FOM and by the POD-G and POD-ANN methods related to the numerical experiment (iv) for Ξ = {2.365,0.6,0.6,0.5,3.2,14.10,8.50,9.2,9.9,10.6,10,2.08𝑒9,1.39𝑒9,1𝑒−6}. We consider 7POD modes. For POD-ANN, we set 𝑛𝑡𝑟 =2500 and 𝐻=130. 23
7. Some concluding remarks In this work we propose a computational pipeline to obtain fast and reliable numerical simulations for one-way coupled steady state linear thermo-mechanical problems in a finite element environment. The test case is referred to a relevant industrial problem related to the investigation of the thermo-mechanical phenomena occurring in blast furnace heart walls. After introducing the main theoretical features of FOM, we detect customized benchmarks for the validation of its numerical implementation. Then we present our MOR framework: we apply POD for the computation of reduced basis space whilst for the evaluation of the modal coefficients we use two different methodologies, the one based on a classic Galerkin projection (POD-G) and the other one based on artificial neural networks (POD-ANN). We found that POD-G is generally more accurate than POD-ANN although POD-ANN exhibits a very higher efficiency, especially when geometric parameters are considered. The higher efficiency of POD-ANN in the case of mechanical model can also be attributed to the fact that the computation of reduced basis approximation of temperature field is not required to compute the reduced basis approximation of the displacement field. We believe that insights given in this work could help to develop advanced numerical tools in order to deal with complex industrial problems. As a follow-up of this work, we are going to enhance training capacity of the deep learning methods, i.e. to reduce the offline cost, as well as to move towards more complex thermo-mechanical problems involving heterogeneity, orthotropy and non-linearity. Some preliminary efforts, in the latter direction, have been already carried out both at full order [53] and reduced order level [52]. Concerning the POD-Galerkin method, it is expected that the computational efficiency furthermore decreases due to absence of affine expansion for assembling system of equations. Regarding the accuracy of the method, the interpolation of operators could introduce an additional source of error. On the other hand, the POD-ANN method can be properly set to take care of non-linearities by increasing the number of hidden layers and/or their depth. However, the accuracy and the computational efficiency are not likely to change to the significant extent [27, 41]. Acknowledgements We are grateful to Dr. Federico Pichi (SISSA mathLab) for insights and crucial support in the numerical implementation of artificial neural network. We would like to acknowledge the financial support of the European Union under the Marie Sklodowska-Curie Grant Agreement No. 765374. We also acknowledge the partial support by the European Union Funding for Research and Innovation - Horizon 2020 Program - in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” and INDAM-GNCS project “Advanced intrusive and non-intrusive model order reduction techniques and applications”, 2019. This work was also partially supported by FEDER and Xunta de Galicia [grant numbers ED431C 2017/60, ED431C 2021/15], and the Agencia Estatal de Investigación [PID2019-105615RBI00/AEI/10.13039/501100011033]. This work has focused exclusively on civil applications. It is not to be used for any illegal, deceptive, misleading or unethical purpose or in any military applications. This includes any application where the use of this work may result in death, personal injury or severe physical or environmental damage. Appendix A. Validation of the full order model In this section we verify the numerical implementation of the FOM introduced in Secs. 4 and 5. All the FOM computations have been performed by using the python finite element library FEniCS [1]. We use a mesh of 𝜔containing 121137 triangular elements and 61147 vertices. The minimum mesh size is 0.011 m and the maximum one is 0.045 m. Its minimum quality is 𝑞𝑒=0.25 (eq. 68). The pipeline that we follow for the design of reliable benchmark tests to be used for the FOM validation consists of three steps: •We set analytical expressions for temperature and displacement. •We calculate corresponding model data, including boundary conditions and source terms, in order to identify the FOM for which the analytical relationships are solutions. •Finally, we numerically solve the problem and compare the computational solutions with the analytical ones. 24
We consider the physical properties reported in Table A.9 for all numerical simulations shown in this section. Property Value Thermal conductivity 𝑘10 𝑊 𝑚𝐾 Convection coefficient ℎ𝑐,−2000 𝑊 𝑚2𝐾 Convection coefficient ℎ𝑐, 𝑓 200 𝑊 𝑚2𝐾 Convection coefficient ℎ𝑐,𝑜𝑢𝑡 2000 𝑊 𝑚2𝐾 Young’s modulus 𝐸5𝑒9𝑃𝑎 Poisson’s ratio 𝜈0.2 Thermal expansion coefficient 𝛼10−6/𝐾 Reference temperature 𝑇0298𝐾 Gravitational acceleration 𝑔9.81 𝑚 𝑠2 Table A.9: Physical properties values used for the FOM benchmark tests. Appendix A.1. Thermal model We consider the following analytical expression for the temperature, 𝑇𝑎(𝑟, 𝑦)=𝐶0𝑟2𝑦, with 𝐶0=1𝐾/𝑚3.(A.1) Then: •The corresponding source term 𝑄is obtained by using eq. (16), 𝑄(𝑟, 𝑦)=−𝑘𝜕2𝑇𝑎 𝜕𝑟2−𝑘𝜕2𝑇𝑎 𝜕𝑦2−𝑘 𝑟 𝜕𝑇𝑎 𝜕𝑟 =−4𝐶0𝑘𝑦 , (A.2) •The heat flux 𝑞+, as well as the temperatures 𝑇𝑓, 𝑇𝑜𝑢𝑡 and 𝑇−, are derived from eq. (17), on 𝛾+:𝑞+(𝑟, 𝑦)=−𝑘𝜕𝑇𝑎 𝜕𝑦 =−𝐶0𝑘𝑟2,(A.3a) on 𝛾𝑠 𝑓 :𝑇𝑓=𝑇𝑎+𝑘 ℎ𝑐, 𝑓 𝜕𝑇𝑎 𝜕𝑟 𝑛𝑟+𝜕𝑇𝑎 𝜕𝑦 𝑛𝑦 =𝐶0𝑟2𝑦+𝐶0𝑘 ℎ𝑐, 𝑓 (2𝑟𝑦𝑛𝑟+𝑟2𝑛𝑦),(A.3b) on 𝛾𝑜𝑢𝑡 :𝑇𝑜𝑢𝑡 =𝑇𝑎+𝑘 ℎ𝑐,𝑜𝑢𝑡 𝜕𝑇𝑎 𝜕𝑟 =𝐶0𝑟2𝑦+𝐶02𝑟𝑦𝑘 ℎ𝑐,𝑜𝑢𝑡 ,(A.3c) on 𝛾−:𝑇−=𝑇𝑎−𝑘 ℎ𝑐,− 𝜕𝑇𝑎 𝜕𝑦 =𝐶0𝑟2𝑦−𝐶0𝑟2𝑘 ℎ𝑐,− ,(A.3d) and it is verified that on 𝛾𝑠:𝜕𝑇𝑎 𝜕𝑟 =0.(A.4a) We solve the (𝑊𝑇)ℎproblem for the data 𝑄, 𝑞+, 𝑇𝑓, 𝑇𝑜𝑢𝑡 , 𝑇−given by equations (A.2)-(A.4a). We choose a discretized space of polynomial of degree 3. Analytical and numerical solutions are reported in Figure A.15 (left and center). As we can see, a very good agreement is obtained. For a more quantitative comparison, we also display the absolute error in Figure A.15 (right), and compute the relative error, ||𝑇𝑎−𝑇ℎ||𝐻1 𝑟(𝜔) ||𝑇𝑎||𝐻1 𝑟(𝜔) =7𝑒−13. 25
[26] Hernández-Becerro, P., Spescha, D., Wegener, K., 2021. Model order reduction of thermo-mechanical models with parametric convective boundary conditions: focus on machine tools. Computational Mechanics 67, 167–184. [27] Hesthaven, J., Ubbiali, S., 2018. Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics 363, 55 – 78. [28] Hesthaven, J.S., Rozza, G., Stamm, B., 2015. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. SpringerBriefs in Mathematics, Springer International Publishing. [29] Hlaváček, I., 1989. Korn’s inequality uniform with respect to a class of axisymmetric bodies. Aplikace Matematiky 34, 146–154. [30] Hlaváček, I., 1989. Shape optimization of elastic axisymmetric bodies. Aplikace matematiky 34, 225–245. [31] Hoang, K.C., Kim, T.Y., 2017. Fast and accurate two-field reduced basis approximation for parametrized thermoelasticity problems. Finite Elements in Analysis and Design 141, 96–118. [32] Huynh, D., Patera, A., 2007. Reduced basis approximation and a posteriori error estimation for stress intensity factors. International Journal for Numerical Methods in Engineering 72, 1219–1259. [33] Iman, R.L., 2008. Latin Hypercube Sampling. John Wiley & Sons, Ltd. [34] Kingma, D.P., Ba, J., 2015. Adam: A method for stochastic optimization, 3rd International Conference on Learning Representations (ICLR). arXiv:1412.6980. [35] Li, H., 2011. Finite element analysis for the axisymmetric laplace operator on polygonal domains. J. Computational Applied Mathematics 235, 5155–5176. [36] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., Anandkumar, A., 2020. Neural operator: Graph kernel network for partial differential equations. arXiv:2003.03485. [37] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., Anandkumar, A., 2021. Fourier neural operator for parametric partial differential equations. arXiv:2010.08895. [38] Meneghetti, L., Demo, N., Rozza, G., 2021. A dimensionality reduction approach for convolutional neural networks. arXiv:2110.09163. [39] Necas, J., 1967. Les Methodes Directes en Theorie Des Equations Elliptiques Jindrich Necas. Masson. [40] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S., 2019. Pytorch: An imperative style, high-performance deep learning library, in: Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., Garnett, R. (Eds.), Advances in Neural Information Processing Systems 32. Curran Associates, Inc., pp. 8024–8035. [41] Pichi, F., Ballarin, F., Rozza, G., Hesthaven, J.S., 2021. An artificial neural network approach to bifurcating phenomena in computational fluid dynamics. arXiv:2109.10765. [42] PyTorch, URL: www.pytorch.org. Accessed: 04-December-2021. [43] Quarteroni, A., Manzoni, A., Negri, F., 2016. Reduced Basis Methods for Partial Differential Equations. Number v. 1 in La Matematica per il 3+2, Springer International Publishing. [44] Raissi, M., Perdikaris, P., Karniadakis, G., 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686–707. [45] Raissi, M., Yazdani, A., Karniadakis, G.E., 2020. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science 367, 1026–1030. [46] RBniCS, URL: www.rbnicsproject.org. Accessed: 04-December-2021. 32
[47] Rojas, R., 1996. Neural Networks: A Systematic Introduction. Springer-Verlag, Berlin, Heidelberg. [48] Rozza, G., Huynh, D., Patera, A., 2007. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Archives of Computational Methods in Engineering 15. [49] San, O., Maulik, R., Ahmed, M., 2019. An artificial neural network framework for reduced order modeling of transient flows. Communications in Nonlinear Science and Numerical Simulation 77, 271–287. [50] Schilders, W., 2008. Introduction to Model Order Reduction. Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 3–32. [51] Schilders, W.H.A., Lutowska, A., 2014. A Novel Approach to Model Order Reduction for Coupled Multiphysics Problems. Springer International Publishing, Cham. pp. 1–49. [52] Shah, N.V., Girfoglio, M., Barral, P., Rozza, G., Quintela, P., Lengomin, A., 2022. Coupled parameterized reduced order modelling of thermomechanical phenomena arising in blast furnaces. Ph.D. thesis. Scuola Internazionale Superiore di Studi Avanzati. URL: hdl.handle.net/20.500.11767/127929. [53] Shah, N.V., Girfoglio, M., Rozza, G., 2021. Thermomechanical modelling for industrial applications. arXiv:2108.13366. [54] Sinha, T., Sikka, K., Lall, R., 2021. Artificial neural networks and bayesian techniques for flip-chip package thermo-mechanical analysis, 2021 IEEE 71st Electronic Components and Technology Conference (ECTC), pp. 1442–1449. [55] Swartling, M., Sundelin, B., Tilliander, A., Jönsson, P.G., 2010. Heat transfer modelling of a blast furnace hearth. steel research international 81, 186–196. [56] Vizzaccaro, A., Givois, A., Longobardi, P., Shen, Y., Deü, J.F., Salles, L., Touzé, C., Thomas, O., 2020. Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using threedimensional finite elements. Computational Mechanics 66, 1293–1319. [57] Vázquez-Fernández, S., García-Lengomín Pieiga, A., Lausín-Gónzalez, C., Quintela, P., 2019. Mathematical modelling and numerical simulation of the heat transfer in a trough of a blast furnace. International Journal of Thermal Sciences 137, 365 – 374. [58] Wang, Q., Hesthaven, J.S., Ray, D., 2019. Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem. Journal of Computational Physics 384, 289–307. [59] Zhang, G., Eddy Patuwo, B., Y. Hu, M., 1998. Forecasting with artificial neural networks:: The state of the art. International Journal of Forecasting 14, 35–62. [60] Zhang, S., Oskay, C., 2017. Reduced order variational multiscale enrichment method for thermo-mechanical problems. Computational Mechanics 59, 887–907. [61] Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z., 2005. The Finite Element Method: Its Basis and Fundamentals. 6 ed., Butterworth-Heinemann. 33