Mathematical modelling and computational reduction of molten glass fluid flow in a furnace melting basin

Computational and Applied Mathematics (2024) 43:407 https://doi.org/10.1007/s40314-024-02921-1 Mathematical modelling and computational reduction of molten glass fluid flow in a furnace melting basin Francesco Ballarin1·Enrique Delgado Ávila2·Andrea Mola3·Gianluigi Rozza4 Received: 11 September 2023 / Revised: 15 July 2024 / Accepted: 25 August 2024 / Published online: 13 September 2024 © The Author(s) under exclusive licence to Sociedade Brasileira de Matemática Aplicada e Computacional 2024 Abstract In this work, we present the modelling and numerical simulation of a molten glass fluid flow in a furnace melting basin. We first derive a model for a molten glass fluid flow and present numerical simulations based on the finite element method (FEM). We further discuss and validate the results obtained from the simulations by comparing them with experimental results. Finally, we also present a non-intrusive proper orthogonal decomposition (POD) based on artificial neural networks (ANN) to efficiently handle scenarios which require multiplesimulationsofthefluidflow upon changing parameters of relevant industrial interest. This approach lets us obtain solutions of a complex 3D model, with good accuracy with respect to the FEM solution, yet with negligible associated computational times. Keywords Molten glass flow ·Finite element method ·Artificial neural networks ·Proper orthogonal decomposition 1 Introduction Efficiently addressing computational fluid dynamics (CFD) problems arising from industrial problems with an appropriate numerical model in an affordable computational time is typically very challenging. While every application often comes with its own set of challenges, a BEnrique Delgado Ávila [email protected] Francesco Ballarin [email protected] Andrea Mola [email protected] Gianluigi Rozza [email protected] 1Department of Mathematics and Physics, Università Cattolica del Sacro Cuore, Brescia, Italy 2Departamento de Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, 41080 Seville, Spain 3MUSAM Lab, Scuola IMT Alti Studi Lucca, Piazza S. Ponziano, 6, 55100 Lucca, Italy 4mathLab, Mathematics area, SISSA, International School for Advances Studies, via Bonomea 265, 34136 Trieste, Italy 123 407 Page 2 of 21 F. Ballarin et al. common one is that the complexity of industrial mathematical models makes them impossible to solve in real time, especially in those cases where the model itself depends on several parameters. This work aims to exemplify the application of a non-intrusive reduced order modelling technique to a specific industrial application, namely the fluid flow of molten glass. In order to achieve this goal, first of all we will present a mathematical model for molten glass fluid flow in a furnace melting basin. We will discuss in particular which physics are the most relevant for the application at hand, and therefore must be incorporated by means of appropriate partial differential equations in the model; the remaining physical contributions will be included instead by means of forcing terms or boundary conditions. After this, we will introduce a post-processing stage needed to compute the performance parameters desired by design engineers and furnace operators. The molten glass is assumed to be a viscous fluid, which is heated in two ways: by methane combustion and by electric booster heating. The two heat sources will be modeled differently: on one hand, the heat flux due to the methane combustion flame is incorporated in the energy of the system by means of a boundary condition. On the other hand, a partial differential equation will be set up to model the booster heating, by considering an electrical potential equation. Thus, we are dealing with a (non-linear) multiphysics system, where the velocity, pressure and temperature of the fluid are coupled with the electric potential produced by the boosters. The resulting model is more complete than the one presented by previous authors, who discuss more basic models for molten glass, without the consideration of boosting heating (Abbassi and Khoshmanesh 2008;Maetal.2018; Pilon et al. 2006a,b;Shelby2007; Viskanta 1994; Xiqi and Viskanta 1986). Even though modelling of booster heating has already been presented in a few previous works (Chang et al. 2002; Choudhary 1985; Choudhary et al. 2010; Curran 1973; Giessler and Thess 2009), most of the simulations performed therein have been made by Finite Difference Methods (FDM). Here we consider instead the Finite Element Method (FEM), since it allows us to consider more complex geometries than those previously discussed with FDM, and enables us to perform local refinement close to specific regions of interest in the domain. The resulting FEM model will then be validated against experimental measurements. While the proposed FEM model allows us to obtain reliable simulations against experimental data, the CPU time required to run each simulation results in a computational model which is too expensive to be queried in any practical scenario, which ideally would require a real time response (e.g., to monitor the furnace conditions and take appropriate action in case of undesired situations). To this end, we identify five parameters that mimic a change in operating conditions in the furnace: one parameter corresponds to the amount of energy coming in the system by the methane combustion, and the other four correspond into the voltage of the boosters. Thus, we are parameterizing the energy that is incorporated to the system. In literature, many authors have proposed the use of Reduced Order Models (ROM) in order to reduce the computational time (Quarteroni 2014; Quarteroni et al. 2015,2011). Nowadays,twofamiliesofROMsaremainlystudied:theintrusiveROMandthenon-intrusive ROM. The classic intrusive ROMs consist of solving the reduced problem by performing a Galerkin projection onto the reduced spaces to compute the value of the modes weights, while the most recent non-intrusive ROMs approximate the modes weights, that only depend on the parameter, by interpolation or regression techniques. In order to achieve a real time visualization of the furnace state, we rely on a non-intrusive ROM based upon Artificial Neural Networks (ANNs) (Chen et al. 2021; Hesthaven and Ubbiali 2018) for the molten glass fluid flow model. During a training phase, we construct the reduced spaces by a Proper Orthogonal Decomposition (POD). During the evaluation 123 Mathematical modelling and computational reduction of molten… Page 3 of 21 407 phase, we seek a reduced order solution of the problem as linear combination of the POD modes, where an ANN is queried to produce the weights to be employed in the combination. On one hand, by foregoing the use of an intrusive ROM we do not achieve certification of the error committed by the ROM (Ballarin et al. 2020; Chacón Rebollo et al. 2017; Kunisch and Volkwein 2002; Manzoni 2014; Novo and Rubino 2021); on the other hand, since we do not have to consider a Galerkin projection of the model equations onto the reduced spaces, this non-intrusive ROM lets us compute on-line solutions for this model that has highly non-linear terms, avoiding the use of techniques such as the EIM or DEIM (Chaturantabut and Sorensen 2010;Barraultetal.2004). Non-intrusive ROMs (Chen et al. 2021; Hesthaven and Ubbiali 2018; Pichi et al. 2023) are being increasingly considered in recent years when handling complex problems, and they have proven to be able to solve almost in real-time realistic industrial problems, since we avoid the treatment of the non-linear terms of the model equations. Other non-intrusive techniques can be addressed in this context. In Ivagnes et al. (2023) and Tonicello et al. (2024) authors consider, in addition to ANN, the Radial Basis Functions (RBF) technique for the reconstruction of the POD coefficients, applied to CFD. In those works, authors compare the error between both techniques for the reconstruction of the POD coefficients, showing that the error of the POD-NN and POD-RBF are comparable. To summarize, this work presents an industrial application, which results in a nonlinear PDE problem coupling thermal, electrical and fluid dynamic fields in a complex domain. The spatial discretization has been successfully tackled with with Finite Elements Method, which resulted in a significantly lower numbers of the degrees of freedom of the numerical problem, with respect to the ones that would have been obtained with either Finite Differences or Finite Volumes methods typically used in these applications. Tomaximise the information provided todesignengineersand furnace operators, a post-processing step has beendevisedto compute the motion of air bubbles which are released by the melting glass and might reduce the final product quality if present at the throat upon extraction from the furnace. Given its elevated computational cost, such a complicated problem represents an ideal—albeit challenging— test bench for model order reduction techniques. In particular, the goal of providing real time predictions to furnace operators has been achieved by means of a non-intrusive ROM combining POD for the computation of modal basis and Artificial Neural Networks (ANNs) for the computation of the basis coefficients. Documenting the most critical aspects and the results obtained with the academic methods described when put to test on the present industrial problem represent critical information for the community. Hijazi et al. (2020) The structure of this paper is as follows. In Sect. 2, we present the problem modelling, with its mathematical formulation where the boundary conditions are also presented. In Sect. 3, we present the Finite Element discretization of the molten glass flow problem, and the Proper Orthogonal Decomposition Neural Network (POD-NN) reduced order method. In Sect. 4, is presented the validation of the Finite Element problem with experimental data, and a post-processing of the solution obtained and the numerical solution of the PODNN. 2 Fluid flow modelling of a glass furnace melting basin In the present section, we discuss the modeling aspects of the container glass furnace melting basin simulations carried out in this work, namely a molten fluid glass flow. 123 407 Page 4 of 21 F. Ballarin et al. 2.1 Geometry of an industrial melting basin The melting basin studied herein is depicted in Fig. 1, and will be denoted in the following as the spatial domain . The melting basin reported in the Figure is a slightly simplified geometry of the actual basin employed by our industrial partners Bormioli Pharma Srl. We will provide a description of each part of the basin geometry next, yet refrain from reporting precise quantitative information on the geometry itself since such information is privileged and covered by the intellectual property of the industrial partner. However, to allow for the qualitative reproduction of the results reported, we report that the overall length of the furnace considered is of the order of tenths of meters. The inflow of raw material occurs through the inflow section denoted as in in Fig. 1.The inflow section lies on top of a region, called the doghouse region in industrial practice, which in Fig. 1has the shape of a trapezoidal prism and which goal is to collect all the material that enters in the basin. The material then moves to the actual furnace, whose is modelled in Fig. 1as the parallelepiped delimited by w,top and bott .bott and ware the walls of the basin, placed respectively at its bottom and on its vertical sides. top is not an actual physical wall, but the free surface that the molten glass assumes during the melting process; for reasons that will be discussed in the following, it is reasonable in our model to assume such a surface to be flat. Furthermore, eight pairs of boosters are installed on the bottom surface; the surface of each booster is denoted by bi,fori=1,...,8. Finally, a barrier wris placed towards the end of the furnace to create a sudden contraction in the basin, with the goal of enhancing the recirculation of the molten glass before extraction of the finished product. The throat, represented by the elongated parallelepiped delimited by duct and out in Fig. 1, constitutes the channel downstream of the furnace by which the finished product is extracted; in particular, the surface of the throat is denoted by duct and acts as a physical wall, while the outflow section of the finished product by out . The different boundaries of the domain are indicated in Fig. 1. Fig. 1 Industrial melting basin employed in this work. The boundaries indicated in the figure will be referred in this work as in:Inflow of material, top:Top of the furnace, bi:Surface of booster i, w:The walls of the furnace, bott :The bottom of thefurnace, wr:The wall for the recirculation, duct:The walls of the throat, out :The exit of the throat 123 Mathematical modelling and computational reduction of molten… Page 5 of 21 407 2.2 Rationale and desiderata for the fluid flow model The main goal of the simulation tool developed herein is to model physical processes occurring in the melting basin in order to assess how furnace design parameters affect the quality of the molten glass extracted from the furnace at the end of the melting process. The main aspect of glass quality typically monitored in industrial practice is the presence of air and gas bubbles within the molten glass. In fact, such bubbles could be visible in the final product if they were still present in the fluid glass at the moment when the material is eventually molded into the desired containers, and thus they would become an undesired flaw. The presence of bubbles within the fluid glass is partly due to the progressive melting of the air-filled granular raw material lying on the basin surface in the doghouse region, and partly due to chemical reactions occurring in the basing and releasing gas products. Thus, in principle a suitable model for the application at hand should be able to reproduce, along with fluid dynamic and thermal fields, the melting mechanism of granular material, as well as the chemical reactions leading to gas release in the fluid. However, given the fact that gas and air bubbles released in the molten glass are quite difficult to prevent, many of the design solutions and control parameters available to engineers and technicians focus instead on obtaining glass residence times in the furnace long enough for the air and gas to be released through the free surface of the basin. The presence of walls and electrodes (or boosters) which introduce heat through electric current are among the solutions put in place so as to create and control convective and recirculating flow patterns that increase the glass residence time. Therefore, in order to test the effects of new designs and control parameters, in this work we concentrate our efforts on devising a model that results in the accurate reproduction of the molten glass motion in the basin through the mechanisms affecting the electric potential, thermal and fluid dynamic fields, but not the production of air bubbles. However, as it will be detailed in the following sections, once the fluid velocity field has been obtained, we will evaluate the motion of small air bubbles by means of a post-processing step, under the assumption that their presence does not significantly affect the fluid dynamic field. 2.3 Mathematical model Fortheaforementionedreasons,andconsideringthelowReynoldsnumberoftheflowathand, thegoverning equation adopted in this workis thetransient three-dimensional incompressible Navier–Stokes equations, which describe the motion of a viscous Newtonian fluid. To include thermal convection effects even in the presence of an incompressible fluid model, we resort to Boussinesq approximation for the buoyancy terms. Finally, the heat source generated by boosters is considered in the model as a one-way coupling term between the electric potential field and the thermal field. Let ⊂R3be the bounded polyhedral computational domain introduced in Sect. 2.1, and [0,T]⊂Rthe time interval. The proposed mathematical model is governed by the following partial differential equations ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ρ∂tu+ρ(u·∇u)−∇·(μ∇u)+∇p−ρgβ(θ −θ0)=0in ×[0,T], ∇·u=0in×[0,T], ρcp∂tθ+ρcp(u·∇θ) −∇·(kef f ∇θ) =σ|∇φ|2in ×[0,T], −∇·(σ∇φ) =0in×[0,T]. (1a) (1b) (1c) (1d) 123 407 Page 6 of 21 F. Ballarin et al. Table 1 Physical properties of the model βThermal expansion 7.5·10−5[1/K] θ0Reference temperature 1617 [K] ρDensity 2250 [kg/m3] cpSpecific heat 1381 [J/(kgK)] μViscosity exp10425 θ−500 −6.0917[kg/(ms)] kef f Thermal conductivity 1.73 +2.5·10−8θ3[W/(mK)] σElectrical conductivity exp(7.605 −7200/θ) [1/(m)] The unknowns of the model are the fluid velocity u[m/s], the fluid pressure p[Pa], the fluid temperature θ[K], and the electric potential φ[V]. Here, gdenotes the gravity. The expressions for the physical properties appearing in (1) are summarized in Table 1,whereρ, cp,μ,β,kef f ,andσare, respectively, the density, the specific heat, the viscosity, the thermal expansion coefficient, the thermal conductivity, and the electrical conductivity of the fluid. Weremark thatviscosity,thermal andelectricalconductivityshowdependenceonthe local temperature. Their expressions as reported in Table 1are based upon the pyshical properties considered in Choudhary (1985). We point out that the expression for the viscosity, μ(θ) =exp 10425 θ−500 −6.0917kg ms(2) presents a vertical asymptote at θa=500K, which would prevent to reach low temperatures for the material, as the viscosity would become infinite as the temperature approaches θa. Considering that the viscosity model reported in Table 1is no longer valid for low temperatures as it would be nonphysical for the viscosity to reach an arbitrarily large value, we set a threshold at θc=973K, and modify the viscosity expression as follows: if θ≥θc, than the viscosity law is modelled as in Table 1, otherwise we assume a constant viscosity equal to the one obtained at θc. Moreover, the power at the electrodes is modeled by a source term in the energy equation, as σ(θ)|∇φ|2. This corresponds with the heat generated by Joule effect by the boosters. Taking a look at Eq. 1, it would be in principle possible to consider a one-way coupling between electric potential and temperature field, σshowed only dependence on θ. In this work, variations of electrical conductivity as a function of temperature (cf. Choudhary (1985)) are expressed by σ(θ) =exp(7.605 −7200/θ) [1/( m)].(3) In the temperature range T∈[973◦K;1900◦K]observed in the furnace, the value of σthen varies between 1.2−1m−1and 44.7−1m−1. Given this extreme temperature sensitivity of the electrical conductivity, Eq. 1d is strongly affected by temperature changes, which can only be tackled considering a coupling between the fields. 2.4 Boundary conditions We present next the boundary conditions accompanying the partial differential equations (1). We discuss separately boundary conditions for velocity/pressure, temperature and electric equations. 123 Mathematical modelling and computational reduction of molten… Page 7 of 21 407 2.4.1 Velocity and pressure boundary conditions For the boundary conditions associated with the Navier–Stokes momentum equation, we consider a prescribed inlet velocity in the normal direction to in. Assuming the normal to in in Fig. 1to be aligned with the vertical direction ey, the velocity at the inlet is given by: uin =uyeyon in.(4) Since our simulations start from an empty furnace, we consider a time-modulated function uy=uy(t), in which during the first 72h we gradually increase the quantity of material introduced in the furnace, after which the desired inflow velocity is reached and kept constant. Thus, the function uy(t)in (4)isdefinedas uy(t)=⎧ ⎨ ⎩ 1 2sin πt 72 ·3600 −π 2+1 2¯uyif t≤72 ·3600 s ¯uyif t>72 ·3600 s, (5) where, in order to respect no slip boundary conditions at the surrounding walls, ¯uyis a parabolic profile scaled in such a way to obtain a mass flow of 13.7tons per day. With respect to the side walls, bottom, the wall for the recirculation, throat, and boosters walls, the velocity is set to zero (no-slip condition), thus, u=0on w∪bott ∪wr∪duct ∪bi. For the top of the furnace, since there are no walls on the top boundary, we have considered slip boundary conditions for the velocity, that is u·n=0ontop. Finally, at the outlet, we consider the outflow conditions, given by ∂nu−pn=0on out , where we are denoting by nthe outward normal vector, and by ∂nthe normal derivative. 2.4.2 Thermal boundary conditions Concerning the temperature field boundary conditions, we have considered a constant temperature on the inlet boundary. Thus, we consider θ=θin on in. In particular, we consider that the temperature on the inlet boundary is fixed at θin =1073K. On the top boundary, we consider a heat flux boundary condition, given by kef f ∂nθ=qtop on top. where qtop is a spatial function for the thermal boundary condition on the top boundary. This function is qtop =max{Q(z−zmin)(z−zmax)(x−xmax)(x−xmin), 0},(6) where xmin,xmax,zmin,and zmax are the coordinates of the boundaries of the basin top surface hit by the flame, and Qis computed such that the total heat flux is 269,662kW/(m2s). 123 407 Page 8 of 21 F. Ballarin et al. This flux intends to model the heat flux of methane combustion flame as a parabolic function along both xand ycoordinates. We point out that making use of different values of xmin,xmax,zmin,and zmax is is also possible to deform the shape of the flame. On the side walls and on the bottom of the furnace, we consider convection boundary conditions corresponding to Robin boundary conditions. By supposing the external temperature to be θext =300K, these conditions are given by kef f ∂nθ=hw(θ −θext)on w, kef f ∂nθ=hbott (θ −θext)on bott , and kef f ∂nθ=hduct(θ −θext)on duct, with hw=6.123W/m2K, hbott =0.531W/m2K, and hduct =0.05W/m2K respectively. These heat transfer coefficients depend on the insulation material of each wall considered. Finally, both on the outflow and the boosters we have considered adiabatic boundary conditions, that is ∂nθ=0onout ∪bi. 2.4.3 Electric boundary conditions Concerning the electric boundary conditions, we consider a fixed potential en each electrode, that is φ=˜ Vion bii=1,...,8.(7) These ˜ Viare imposed so as we have the same Vat each pair of electrodes. In particular, we consider that each pair of booster has ˜ V2j−1=Vj/2and ˜ V2j=−Vj/2, j=1,2,3,4, with V1=150.6V, V2=132.4V, V3=112V, V4=124.6V, corresponding from the first pair of boosters to the fourth, respectively. Finally,weconsider homogeneous Neumann boundary condition,on the rest of theboundary, i.e., ∂nφ=0onin ∪out ∪w∪top ∪wr∪bott ∪duct. 2.5 Air bubbles motion as a post-processing Once the glass velocity and thermal field are computed, it is possible to make use of the flow field to predict the motion of air bubbles in the molten glass. Assuming that the air bubbles are spheres with a radius rthat is smaller than the scales of the velocity field computed, the relative motion of the sphere in the flow field results equivalent to the simpler problem of a sphere moving in a fluid at rest. In addition, given the small relative velocity between the bubble and the fluid, we will assume that air bubbles’ motion can be modelled under Stokes flow and null inertial contribution. In such a case, the air bubble is subjected to the hydrostatic buoyancy force and to the fluid dynamic resistance, which in the case of Stokes flow has closed form. The force balance along the vertical direction for an air bubble of radius rrising with relative speed Vwith respect to the surrounding glass then reads 4 3πr3g(ρg(θ) −ρa(θ)) =6πrμ(θ)V(8) 123 Mathematical modelling and computational reduction of molten… Page 9 of 21 407 where ρa(θ) is the air density. ρga(θ) is the glass density and gthe gravity acceleration. The equation results in an expression for the vertical velocity of the rising bubble V=2 9πr2g(ρg−ρa) μ,(9) where the temperature dependency of the densities and of the glass viscosity is henceforth omitted to lighten the notation. In this work, however, we used a slight modification of this formula, namely Vb=gρr2 3μ,(10) based on experimental evidence reported in Shelby (2007), where ρ =ρg−ρa.Eq.(10) applies to bubbles with radius smaller than 1mm. Making use of the glass density ρgand viscosity μvalues reported in Table 1, the Reynolds number ReVb=2 9ρ2 μ2gr3obtained for 1mm bubbles is significantly smaller than 1, which justifies Stokes flow assumptions. Thus, the experimental correction in Eq. (10) is more likely justified by a non perfectly spherical shape of the rising bubbles measured in the experiments. As said, the bubble rising displacement is the only relative motion between the air sphere and the surrounding fluid. For such reason, the local air bubble velocity field is obtained using the local horizontal components of the fluid velocity, and obtaining the vertical component as the sum of Vband the local fluid velocity vertical component. 3 Numerical model for full order and reduced order solver In the present section, we will first describe the numerical solver based on Finite Element Method (FEM) used for the resolution of the mathematical model equations described in the previous section. After this, we will also provide details of the fully data driven POD-NN reduced order model implemented to accelerate the molten glass flow simulations. 3.1 Finite element modelling To discretize the mathematical problem (1) by means of the Finite Element Method, we start defining the variational formulation of problem (1). We consider the space of square integrable functions L2() and the Sobolev space H1(), equipped with the standard norm. For simplicity in the notation, let us denote Dv=w∪bott ∪wr∪duct ∪bi. Thus, we consider the following spaces for velocity, pressure, temperature and electric potential, respectively W={v∈[H1()]3:v·n=0ontop,v=0onDv,v=uin on in}, M={q∈L2()},={z∈H1() :z=θin on in}, and, ={ω∈H1() :ω=˜ Vion bi}. 123 407 Page 16 of 21 F. Ballarin et al. in the furnace. Heat is in fact released in the molten glass basin both through methane combustion in the upper chamber, and through Joule effect associated with the electric currents generated by the boosters. Thus, in the numerical test described in the present section, we consider that the total heat flux of combustion, qtop (see Eq. (6)) ranges in the interval [240,000, 300,000]kW/(m2s), while the voltage (see Eq. (7)), of each of the four pairs of boosters ranges in the interval [100,160]V. With this configuration, we are assuming that our parametrized PDE depend on five parameters, which vary in the parameter space D=[240,000,300,000]×[100,160]4.A total of 200 FOM simulations have been carried out so as to assemble and validate the PODNN model described. The results of the first 150 snapshots were obtained with parameters in the training set, namely μ∈Dtrain. These solutions were first used to produce the snapshot matrix from which the POD modes were extracted, and then train the feed-forward Neural Network used to compute the modal coefficients at the reduced order level. The residual 50 solutions computed, associated with the testing parameters μ∈Dtest have then been used for validation purposes. We have computed the FOM solutions in a cluster with IBM nodes: Xeon E5-2680 v2 of 40GB of RAM. We have computed the FOM solution with 40 processors, and each FOM solution takes about 3h to be computed each one. Thus, to perform the offline phase (FOM computation plus training of the Neural Network), it took about one month to be completed. Then, the online phase, as can be viewed at the website https://youtu.be/yRNQaGVmXF4, it is performed on real time. Figure 7displays a comparison between the ROM and FOM solutions obtained for a parameter setting corresponding to 283,961.719kW/(m2/s) of combustion heat flux, and 146.681V, 116.542V, 150.134V, and 106.326V for the voltage of each pair of booster, respectively. The plot on the top depicts contours of temperature field predicted with the FOM solver, on three vertical planes within melting the basin. Using an identical layout, the middle plot illustrates the the contours of the temperature field computed by means of the ROM algorithm. Finally, the bottom plot depicts the contours of the relative error between ROM and FOM temperature on the three vertical planes used in the temperature plots. We can observe the similarity of both FOM and ROM solution. Actually, in this case, the maximum of the local relative error of both solutions is about 1%. In addition, we compute the overall POD-NN relative error, eu,ep,eθ,eφ, of each field for a given parameter value μ∈D,as eu(μ)=uh(μ)−uNN(μ) uh(μ),ep(μ)=ph(μ)−pNN(μ) ph(μ), eθ(μ)=θh(μ)−θNN(μ) θh(μ),eφ(μ)=φh(μ)−φNN(μ) φh(μ). We consider a set of parameter values, Dtest, for which we compute the POD-NN solution relative error with respect to the full order solutions. In Fig. 9, we show the mean of the PODNN relative errors for velocity, pressure, temperature and electric potential, depending on the number of POD modes we have considered to construct the reduced spaces. The plot shows that, among the variables considered, velocity is the one comparing the order of magnitude of them. We can observe that the POD-NN relative error, with respect to the high fidelity solution, ranges between 6 and 1%, noticing that the error in temperature is the lowest one. We can observe also a high variability on the mean of the relative error depending on the number of POD modes. That means that we do not need a high number of POD modes in order to have 123 Mathematical modelling and computational reduction of molten… Page 17 of 21 407 Fig. 7 From top to bottom: Temperature of FOM solution, temperature of ROM solution, and the relative error between FOM and ROM of the temperature field. An animation of the solution variations resulting from modification of the value of different process parameters can be viewed at the website https://youtu.be/ yRNQaGVmXF4 a good approximation, with respect to the high fidelity solution, with the POD-NN approach considered. In addition, the plot clearly shows that the POD-NN approximation error with respect to the FOM solutions is growing as the number of modes considered in the ROM increases. Of course, this might seem a counter intuitive occurrence, as the error between ROM and FOM is ideally expected to be reduced by considering more modes in the ROM. However, we must remark that in this work we only had 200 FOM solutions available, 150 of which were used for training. Thus, as we increased the number of modes considered, and as the dimension of 123 407 Page 18 of 21 F. Ballarin et al. 051001050 10-20 10-15 10-10 10-5 100Nomalized eigenvalues Velocity eigenvalues Pressure eigenvalues Temperature eigenvalues Electric potential eigenvalues Fig. 8 Normalized eigenvalues for velocity, pressure, temperature and electric potential 150100500 POD modes 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Velocity error Pressure error Temperature error Electric potential error Fig. 9 Comparative of relative errors for velocity, pressure, temperature and electric potential the NN output space was growing, we could not increase accordingly the amount of samples usedfortraining.In lightof this consideration, the observeddegradationoftheROMaccuracy has to be ascribed to a loss of performance of the NN training. In Fig. 8we show the decay of the normalized eigenvalues for each field. We can observe how at the beginning, the decay is quite fast, while for the last values this decay get slower. Moreover, the decay rate of each field is quite similar. Since the main decay of the eigenvalues is produces for the first nodes, it could be that there is when the performance of the NN training is lost. Clearly, we could have increased the offline computational effort and provide the Neural Network with additional training samples and snapshots to produce decreasing error curves, as the ones observed for the first 10 modes. However, we feel that having a cap on the number of full order simulations is a more realistic scenario in real world applications. In fact, in most situations the additional com- 123 Mathematical modelling and computational reduction of molten… Page 19 of 21 407 putational time—and money—required for each FOM simulation, and even the amount of storage required to save snapshots for each field considered in the PDE solution, put a hard limit to the amount of FOM simulations that can be used for training. A plot like that presented in Fig. 9is then representing a relatively common situation in the realm of industrial applications. 5 Conclusions In this work we have presented the modelling of a 3D flow of a molten glass fluid inside a furnace. We have presented the numerical results of a FEM simulation and compared them with experimental data. We have also modelled the behaviour of air bubbles inside the fluid depending on their size. Moreover, we have presented a POD-NN model in order to reduce the computational time of the numerical solution. We have performed this POD-NN model with different number of POD modes and we have compared the solution obtained from the non-intrusive POD-NN model with respect to the FE solution. For all the fields, we have obtained an error ranging between 6% and 1%, the error in temperature being the lowest one. Also, we have seen that that the consideration of non-intrusive ROM for industrial problems can provide acurrate approximations of the solution for a parametrized setting. As future work, we intend to consider alternative non-intrusive methods, such as a POD-RBF, in order to recover the POD coefficients. In order to improve the Neural Network training, we should need to increase the offline computational effort and provide the Neural Network with additional training samples and snapshots to produce decreasing error curves in the POD-NN model. Acknowledgements We acknowledge the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA). Francesco Ballarin acknowledges the INdAM-GNCS project “Metodi numerici per lo studio di strutture geometriche parametriche complesse” (CUP E53C22001930001), the project “Reduced order modelling for numerical simulation of partial differential equations” funded by Universitá Cattolica del Sacro Cuore and the PRIN 2022 PNRR project “ROMEU: Reduced Order Models for Environmental and Urban flow” (CUP J53D23015960001) funded by MUR. Andrea Mola also acknowledges the funding by project ROMEU (CUP D53D2301888001). Data availability Data subject to third party restrictions. Declarations Conflict of interest Not applicable. References Abbassi A, Khoshmanesh K (2008) Numerical simulation and experimental analysis of an industrial glass melting furnace. Appl Therm Eng 28(5–6):450–459. https://doi.org/10.1016/j.applthermaleng.2007.05. 011 Arnold DN, Brezzi F, Fortin M (1984) A stable finite element for the Stokes equations. Calcolo 21(4):337–344. https://doi.org/10.1007/bf02576171 Ballarin F, Chacón Rebollo T, Delgado Ávila E, Gómez Mármol M, Rozza G (2020) Certified reduced basis VMS-Smagorinsky model for natural convection flow in a cavity with variable height. Comput Math Appl 80(5):973–989. https://doi.org/10.1016/j.camwa.2020.05.013 Barrault M, Maday Y, Nguyen NC, Patera AT (2004) An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. CR Math 339(9):667–672. https:// doi.org/10.1016/j.crma.2004.08.006 123 407 Page 20 of 21 F. Ballarin et al. Chacón Rebollo T, Delgado Ávila E, Gómez Mármol M, Ballarin F, Rozza G (2017) On a certified Smagorinsky reduced basis turbulence model. SIAM J Numer Anal 55(6):3047–3067. https://doi.org/10.1137/ 17m1118233 Chang S-L, Zhou C, Golchert B (2002) A numerical investigation of electric heating impacts on solid/liquid glass flow patterns. In: 8th AIAA/ASME joint thermophysics and heat transfer conference. American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.2002-2879 Chaturantabut S, Sorensen DC (2010) Nonlinear model reduction via discrete empirical interpolation. SIAM J Sci Comput 32(5):2737–2764. https://doi.org/10.1137/090766498 Chen W, Wang Q, Hesthaven JS, Zhang C (2021) Physics-informed machine learning for reduced-order modeling of nonlinear problems. J Comput Phys 446:110666. https://doi.org/10.1016/j.jcp.2021.110666 Choudhary M (1985) Three-dimensional mathematical model for flow and heat transfer in electric glass furnaces. Heat Transfer Eng 6(4):55–65. https://doi.org/10.1080/01457638508939639 Choudhary MK, Venuturumilli R, Hyre MR (2010) Mathematical modeling of flow and heat transfer phenomena in glass melting, delivery, and forming processes. Int J Appl Glass Sci 1(2):188–214. https://doi.org/ 10.1111/j.2041-1294.2010.00018.x Curran RL (1973) Mathematical model of an electric glass furnace: effects of glass color and resistivity. IEEE Trans Ind Appl IA 9(3):348–357. https://doi.org/10.1109/tia.1973.349916 Giessler C, Thess A (2009) Numerical simulation of electromagnetically controlled thermal convection of glass melt in a crucible. Int J Heat Mass Transf 52(13–14):3373–3389. https://doi.org/10.1016/j. ijheatmasstransfer.2009.01.004 Hesthaven JS, Ubbiali S (2018) Non-intrusive reduced order modeling of nonlinear problems using neural networks. J Comput Phys 363:55–78. https://doi.org/10.1016/j.jcp.2018.02.037 Hijazi S, Stabile G, Mola A, Rozza G (2020). In: D’Elia M, Gunzburger M, Rozza G (eds) Non-intrusive polynomial chaos method applied to full-order and reduced problems in computational fluid dynamics: a comparison and perspectives. Springer, Cham, pp 217–240. https://doi.org/10.1007/978-3-030-48721- 8_10 Ivagnes A, Demo N, Rozza G (2023) Towards a machine learning pipeline in reduced order modelling for inverse problems: neural networks for boundary parametrization, dimensionality reduction and solution manifold approximation. J Sci Comput. https://doi.org/10.1007/s10915-023-02142-4 Kunisch K, Volkwein S (2002) Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal 40(2):492–515. https://doi.org/10.1137/s0036142900382612 Ma Q, Fang H, Shang C, Liu Z, Wang J (2018)Analysis of melt flows in an electric heating furnace for quartz glass synthesis. In: Volume 7: fluids engineering. American Society of Mechanical Engineers. https:// doi.org/10.1115/imece2018-86112 Manzoni A (2014) An efficient computational framework for reduced basis approximation and a posteriori error estimation of parametrized Navier–Stokes flows. ESAIM Math Model Numer Anal 48:1199–1226 Novo J, Rubino S (2021) Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows. SIAM J Numer Anal 59(1):334–369. https://doi.org/10.1137/20m1341866 Pichi F, Ballarin F, Rozza G, Hesthaven JS (2023) An artificial neural network approach to bifurcating phenomena in computational fluid dynamics. Comput Fluids 254:105813 Pilon L, Zhao G, Viskanta R (2006a) Three-dimensional flow and thermal structures in glass melting furnaces part I. Effects of the heat flux distribution. Glass Sci Technol 75(3):115–124 Pilon L, Zhao G, Viskanta R (2006b) Three-dimensional flow and thermal structures in glass melting furnaces part. II. Effects of the heat flux distribution. Glass Sci Technol 75(2):55–68 Pinnau R (2008) Model reduction via proper orthogonal decomposition. In: Schilders WHA, Vorst HA, Rommes J (eds) Mathematics in industry. Springer, Berlin, pp 95–109. https://doi.org/10.1007/978- 3-540-78841-6_5 Quarteroni A (2014) Numerical models for differential problems. Springer, Berlin. https://doi.org/10.1007/ 978-88-470-5522-3 Quarteroni A, Manzoni A, Rozza G (2011) Certified reduced basis approximation for parametrized partial differential equations and applications. J Math Ind 1(3):1–49 Quarteroni A, Manzoni A, Negri F (2015) Reduced basis methods for partial differential equations: an introduction. Springer, Berlin Shelby JE (2007) Introduction to glass science and technology. R Soc Chem. https://doi.org/10.1039/ 9781847551160 Tonicello N, Lario A, Rozza G, Mengaldo G (2024) Non-intrusive reduced order models for the accurate prediction of bifurcating phenomena in compressible fluid dynamics. Comput Fluids 278:106307. https:// doi.org/10.1016/j.compfluid.2024.106307 Viskanta R (1994) Review of three-dimensional mathematical modeling of glass melting. J Non-Cryst Solids 177:347–362. https://doi.org/10.1016/0022-3093(94)90549-5 123 Mathematical modelling and computational reduction of molten… Page 21 of 21 407 Xiqi W, Viskanta R (1986) Modeling of heat transfer in the melting of a glass batch. J Non-Cryst Solids 80(1–3):613–622. https://doi.org/10.1016/0022-3093(86)90454-0 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 123