Thermodynamics-informed neural network for recovering supercritical fluid thermophysical information from turbulent velocity data

International Journal of Thermofluids 20 (2023) 100448 Available online 23 August 2023 2666-2027/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect International Journal of Thermofluids journal homepage: www.elsevier.com/locate/ijtf Thermodynamics-informed neural network for recovering supercritical fluid thermophysical information from turbulent velocity data Núria Masclans a, Fernando Vázquez-Novoa b, Marc Bernades a, Rosa M. Badiab, Lluís Jofre a,∗ aDepartment of Fluid Mechanics, Universitat Politècnica de Catalunya . BarcelonaTech (UPC), Barcelona 08034, Spain bDepartment of Computer Science, Barcelona Supercomputing Center, Barcelona 08034, Spain ARTICLE INFO Keywords: Deep learning Thermodynamics-informed neural network Supercritical fluid Turbulent flow ABSTRACT Recent research has highlighted the potential of supercritical fluids under high-pressure transcritical conditions to achieve microconfined turbulence as a result of the thermophysical properties they exhibit in the vicinity of the pseudo-boiling region. This has led to increased interest in understanding their hybrid thermophysical properties when operating near the pseudo-boiling transitioning region. However, despite the potential benefits of microfluidic systems working under transcritical conditions, limited experimental data is available due to the inherent challenges of performing experiments at high-pressure conditions. In addition, traditional experimental methods, such as particle image velocimetry and particle tracking velocimetry, are inadequate for measuring thermophysical properties under such conditions, since they are primarily designed for velocity-related data acquisition. In this regard, this work introduces an efficient thermodynamicsinformed neural network framework for reconstructing thermophysical information from velocity data in high-pressure turbulent transcritical regimes. The proposed model incorporates thermophysical constraints through a thermodynamics-informed loss function consisting of the residual of the real-gas equation of state and integrates boundary conditions into the network’s architecture to ensure their satisfaction. The performance of the proposed framework is evaluated through the analysis of two test cases and compared against non-physically informed models. The results demonstrate the superior accuracy, robustness, and satisfaction of physical constraints achieved by the proposed model, as well as its ability to reconstruct averaged thermophysical profiles and preserve bulk quantities with a relative error reduction of approximately 2×. In addition, the physically-consistent predictions provided by the model enable a more accurate reconstruction of dependent thermophysical properties. 1. Introduction High-pressure supercritical fluids have received significant attention for a wide range of engineering applications over the past years, including gas turbines, supercritical water-cooled reactors, and liquid rocket engines [1,2]. In this regard, Bernades et al. [3,4] have recently proposed a distinctive approach to achieve microconfined turbulence by leveraging the hybrid thermophysical properties of supercritical fluids operating in the vicinity of the pseudo-boiling region. Upon crossing the pseudo-boiling line, supercritical fluids experience rapid variations of thermophysical properties such as density, enthalpy, isothermal compressibility and heat capacity. Consequently, Banuti et al. [5] has proposed a three-regime-model for supercritical pseudoboiling that distinguishes between supercritical liquid-like (LL), twophase-like (TPL) and vapor-like (VL) fluids. According to this model, a (i) supercritical liquid-like fluid presents large densities combined with ∗Corresponding author. E-mail address: [email protected] (L. Jofre). small isothermal compressibilities and is characterized by transport coefficients similar to a liquid, a (ii) supercritical gas-like fluid has smaller density, but larger isothermal compressibility and transport coefficients similar to gases, and a (iii) supercritical two-phase-like fluid includes the smooth transition from liquid-like to gas-like across the so-called pseudo-boiling or Widom line. Therefore, it is of high relevance to accurately characterize the working regime of supercritical fluids to understand their physical behavior and enable the design and optimization of engineering applications & solutions. The utilization of microfluidic systems under transcritical conditions offers significant potential benefits [6,7], yet acquiring experimental data in this regime remains challenging due to the inherent difficulties associated with (i) physical/visual access to microscale devices [8–10] and (ii) the high-pressure and high-temperature conditions considered [11]. To investigate fluid flows under these demanding conditions, https://doi.org/10.1016/j.ijft.2023.100448 Received 22 June 2023; Received in revised form 4 August 2023; Accepted 17 August 2023 International Journal of Thermofluids 20 (2023) 100448 2 N. Masclans et al. microfluidics technology has developed in situ characterization methods that aim to provide non-intrusive space- and time-resolved information of hydrodynamic properties. These include optical characterization approaches, such as optical microscopy and high-speed cameras [12], micro-Particle Image Velocimetry (𝜇-PIV) [13], ultraviolet–visible (UV– Vis) spectroscopy [14,15], infrared spectroscopy [16], Particle Tracking Velocimetry (PTV) [17,18], Laser Doppler Velocimetry (LDV) [19], and Schlieren imaging [20,21]. These techniques are not restricted only to low-pressure subcritical fluids; viz. Valori et al. [22] achieved PIV measurements of instantaneous and time-averaged velocity distributions in a thermally convective supercritical fluid. However, these optical methods are not suitable for directly measuring thermophysical properties such as fluid temperature and density, as they are primarily designed for velocity-related data acquisition. To that end, several strategies have been developed to directly obtain thermophysical information at microconfined scales. These include the use of temperature sensors based on fluorescence microscopes, thermocouples, thermistors and custom-build temperature sensors [23], viscosity measurements by the co-flow method [24], and density determination by measuring the pressure drop between a capillary inlet–outlet flow [25] or by Micro-Electro-Mechanical Systems (MEMS) [26]. Nevertheless, these techniques for the direct measurement of temperature and pressure often provide segmented or averaged measurements and require specific channel geometries and mixing strategies. In particular, direct contact measurements can achieve good accuracy, but inherently affect the targeted thermophysical properties of the microfluidic system. Alternatively, no temperature and/or pressure exchange is caused by the non-contact methods, which are based on capturing the changes in the chemical and physical properties of the fluid, and therefore lead to less accurate measurements. Additionally, compatibility issues may be encountered when integrating these sensors into the fabrication process of the microfluidic device [27]. In contrast, this work proposes a novel approach to recover instantaneous and/or averaged local thermophysical properties by fitting them from experimentally accessible velocity-based data. In the field of thermofluid mechanics, accurate assessment of thermophysical fields from velocity data remains a challenge. While exploratory analyses between relevant thermophysical properties have been carried out for high-pressure transcritical fluids based on physicsbased models [3], an accurate fitting method is still needed to obtain thermophysical-related fields from measurable velocity data. In this regard, traditional fitting methods, such as nonlinear regression approaches, assume a predefined mapping function between input and output variables, limiting their applicability when the function type is unknown. In recent years, however, the field of fluid mechanics has witnessed the increased popularity of Machine Learning (ML) methods for regression, optimization and modeling. These methods have been successfully applied to various flow problems, including data-driven inference of turbulence closure models [28–30], improvement of computational fluid dynamics solvers [31,32], super-resolution and/or fluid flow reconstruction [33–35], modal decomposition and dimensionality reduction [36–38], active flow control [39,40], and extraction of sophisticated flow physics insight [41,42]. These problems have been tackled using various data-driven ML methods, such as support vector machines [28], random forests [28,30], and Bayesian optimization & Gaussian processes [29]. Among these methods, Artificial Neural Networks (ANNs) have raised significant interest due to their ability to handle nonlinear and high-dimensional problems without assumptions about the hidden patterns within the data. Particularly, ANNs excel at handling large datasets, exhibit computational efficiency once trained, and have a strong capacity to capture nonlinear mapping relationships. However, despite their widespread use, ANNs often lack physics interpretability and have limited applicability beyond interpolation tasks, particularly when dealing with new data with different probability distributions than the training data [43,44]. This limitation, which is not a major concern for many modern deep-learning applications such as image classification, benefiting from extensive training datasets, proves challenging for the fluid mechanics community. In particular, given that this work targets high-pressure transcritical experimental applications, ANNs may not be well-suited due to the limited (and noisy) available datasets. Physics-driven mechanistic models guided by physical laws/models can overcome these limitations, but rely on accurate physical models that may not be always available. Hence, bridging the gap between data-driven ML methods and physics-driven mechanistic models is crucial to address this limitations. This has led to the development of Physics-Informed ML (PIML) approaches that integrate domain knowledge and physical constraints into ML algorithms, enabling a more efficient data utilization while enhancing prediction accuracy & stability. In the context of fluid mechanics, PIML offers potential to address the challenges of interpretability and extrapolation associated with traditional (data-driven) ML methods. One promising approach is the use of ANNs to solve Partial Differential Equations (PDEs) [45,46]. These physics-embedded networks incorporate physical information and tensorial properties into the learning process, ensuring adherence to the underlying physics of the problem. Successful examples include Long- Short-Time Memory (LSTM) networks [47,48], equivariant or tensor field networks that exhibit rotation and translation equivariance [49– 51], and Tensor Basis Neural Networks (TBNN) [52] that incorporate selected physical invariances and symmetries into the network architecture. In particular, the Physics-Informed Neural Network (PINN) [53] approach introduced by Raissi et al. has raised considerable attention for solving supervised problems while satisfying physical laws and constraints, such as initial and boundary conditions. Moreover, the capabilities of PINNs have been successfully explored in numerous studies in the field of fluid mechanics. For example, PINNs have been effectively used for solving incompressible Navier–Stokes equations to simulate laminar and turbulent flows from computational Direct Numerical Simulation (DNS) data [54–57], and to densely reconstruct experimental scattered PIV/PTV measurements [58]. PINNs have been also extended to solve inverse design problems by adding hard constraints (hPINN) [59,60], achieving simpler and smoother results than conventional PDE-constrained optimization methods for problems with multiple solutions. Other extensions of PINNs include the fractional PINN (fPINN) to solve space–time fractional advection–diffusion equations [61], and the physics-informed generative adversarial networks (PI-GAN) to solve stochastic differential equations [62]. In this work, therefore, the reconstruction of thermophysical information from velocity data is approached as a thermodynamicsinformed regression task, where the aim is not to solve PDEs as typically performed in PINNs implementations. Instead, the solution strategy is based on incorporating physical knowledge into the model to enhance applicability, robustness, and stability. This can be achieved by introducing hypothesis that respect certain physical constraints or by incorporating constraints into the optimization process, either through the network architecture or as penalizations in the loss function. For example, Fiore et al. [63] have recently utilized both approaches to develop a physics-constrained neural network (PCNN) to predict turbulent heat flux for low-Prandtl-number fluids, while enforcing Galilean and rotational invariances. Zhang et al. [64] proposed a thermodynamicsinformed neural network (TINN) to calculate the phase equilibrium properties of a multi-component multi-phase flow to enhance fluid mixtures simulations, and Hernandez et al. [65] extended the TINN framework for computing the dynamic response of deformable objects due to real-time user interactions. However, although numerous TINN implementations exist in the literature, their utilization has been little exploited for regression tasks in flow problems. In this regard, the objective of this work is to reconstruct three-dimensional (3-D) temperature and density fields from velocity data that accurately satisfy the thermodynamics of wall-bounded high-pressure transcritical fluid turbulence. Special emphasis is placed on the physical consistency of the International Journal of Thermofluids 20 (2023) 100448 3 N. Masclans et al. predicted fields, enabling the calculation of other thermophysical properties and the accurate reconstruction of bulk quantities and profiles. To that end, the paper is organized as follows. In Section 2, a general physics-informed regression approach is presented. This is followed by a detailed description of the problem setup and computational dataset used to train and validate the model in Section 3, introducing the flow physics modeling utilized to generate the high-pressure transcritical fluid turbulence data. Next, in Section 4, the TINN framework for the thermophysical reconstruction problem is thoroughly presented. Results regarding the performance of the proposed method are discussed and compared to non-physically informed approaches in Section 5. The goal is to characterize how the incorporation of physical constraints affects the proposed method accuracy and robustness through evaluating metric results, Joint Probability Density Functions (JPDFs), and reconstructed thermophysical profiles. Finally, the work is concluded and future directions are proposed in Section 6. 2. Physics-informed regression Embedding physical laws into a ML algorithm ensures that the predicted solution lies within the space of admissible solutions described by the equations of fluid motion and initial & boundary conditions of the problem. In this regard, the subsections below present a general physics-informed regression framework as a foundation of the proposed TINN framework presented in Section 4. In particular, this work considers a feed-forward Fully-Connected Neural Network (FCNN) architecture, with physics constraints included both in the architecture and in the loss function. 2.1. Fully-connected neural network Given a general regression problem of 𝑛𝑥features and 𝑛𝑦targets, a regression algorithm aims to find a mapping 𝐲=𝑓(𝐱)that accurately fits the training dataset {𝐱∗,𝐲∗}and extrapolates to the validation dataset {𝐱∗∗,𝐲∗∗}, where 𝐱∈R𝑛𝑥and 𝐲∈R𝑛𝑦. In the present work, a feed-forward FCNN has been considered consisting of a large number of neurons (or units) organized in layers [66]. This distributed architecture is characterized by (i) the dense connection between neurons of adjacent layers, where each neuron is connected to all neurons of the following layer (fully connected), and (ii) the output of each layer is feed-forwarded as the input of the next layer. Specifically, the FCNN mapping can be defined by ℒ𝑙(ℒ𝑙−1) = 𝝈𝑙(𝑾𝑙ℒ𝑙−1 +𝒃𝑙),for 𝑙∈ {1,…, 𝐿},(1) which represents the mapping from the input layer ℒ0∈R𝑛𝑥to the output layer ℒ𝐿∈R𝑛𝑦of a FCNN architecture of 𝐿+ 1 layers. Using multiple layers of transformations, FCNN can capture complex and hierarchical interactions between features. Here, ℒ𝑙∈R𝑛𝑙is the nonlinear map for the l-th hidden layer containing 𝑛𝑙neurons, with 𝝈𝑙the non-linear element-wise activation function, and 𝑾𝑙∈R𝑛𝑙×𝑛𝑙−1 and 𝒃𝑙∈ R𝑛𝑙the weights and biases of the linear transformation, respectively. Hence, the FCNN output prediction for a given input vector can be obtained recursively by  𝒚=ℒ𝑙(𝒙;𝜽),(2) where 𝜽=⋃𝐿 𝑖=0 (𝑾𝑖,𝒃𝑖)collects the weights and biases across the network. Generally, the method’s ability to undergo a certain task, i.e., a regression task, is assessed by the so-called loss function , and the learning process aims to minimizing such loss by formulating the optimization loss problem as  𝜽= arg min 𝜃 (𝒙,𝜽),(3) where  𝜽is the optimal set of weights and biases. Typically, the optimization problem back-propagates the loss gradient along the network layers as their weights 𝑾𝑙and biases 𝒃𝑙are updated, and this iterative process is repeated for each training epoch. After training, the optimized network parameters are set constant, and the predicted vector obtained in the output layer ℒ𝐿can be directly calculated for a given input vector ℒ0from Eq. (1). 2.2. Data-driven loss function Traditionally, one can build a regression model of physical quantities without directly enforcing the satisfaction of any governing physical laws or other domain constraints. This purely data-driven approach aims to fit the regression solution 𝒚(𝒙,𝑺), which in FCNN framework is formulated as 𝒚(𝒙,𝑺) ≈  𝒚(𝒙,𝑺)≜ℒ𝐿(𝒙,𝑺;𝜽),(4) where 𝑺is the intrinsic problem parameters that can include fluid and geometry parameters, and  𝒚(𝒙,𝑺)is the network prediction provided by the last output layer ℒ𝐿. The training process of a purely data-driven regression approach aims to find the set of (sub)optimal FCNN parameters  𝜽such that the discrepancy between training data 𝒚𝑑(𝒙,𝑺)and predictions  𝒚(𝒙,𝑺)is locally minimized by the optimization problem as  𝜽= arg min 𝜃 𝐷(𝒙,𝑺;𝜽),(5) 𝐷(𝒙,𝑺;𝜽) = ‖‖‖𝒚𝑑(𝒙,𝑺) − ℒ𝐿(𝒙,𝑺;𝜽)‖‖‖𝛺,(6) where 𝐷(𝒙,𝑺;𝜽)is the data-driven loss function, and ‖⋅‖𝛺is the 𝐿2 norm or prediction mean-squared error over the data domain 𝛺. 2.3. Physics-informed loss function The PINNs framework describes a multi-object loss function that leverages the known PDEs, constitutive relationships, initial conditions (ICs) and boundary conditions (BCs) for a physics-informed optimization. As a result, the associated physics-based loss function 𝑃(𝒙,𝑺;𝜽,𝜶) can be formulated as 𝑃(𝒙,𝑺;𝜽,𝜶) = 𝛼𝑝𝑝(𝒙,𝑺;𝜽) ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ PDEs loss +𝛼𝑐𝑐(𝒙,𝑺;𝜽) ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ Constitutive loss +𝛼𝑏‖‖‖(𝒙,𝒚𝑝)−(𝒙, 𝒚;𝜽)‖‖‖𝛤,𝑡 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ BCs loss +𝛼𝑖‖‖‖(𝒙,𝒚𝑝)−(𝒙, 𝒚;𝜽)‖‖‖𝛺,𝑡0 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ICs loss , (7) where the coefficients 𝛼𝑝, 𝛼𝑝, 𝛼𝑏, and 𝛼𝑖are the weights associated with the PDEs, constitutive models, and BCs & ICs losses, respectively. In detail, 𝜶= (𝛼𝑝, 𝛼𝑐, 𝛼𝑏, 𝛼𝑖)represents the coefficients vector,  𝒚is the network prediction, 𝒚𝑝is the prescribed problem solution in the boundaries of the domain 𝛤for BCs formulated as , and at initial conditions 𝑡0in the domain 𝛺for ICs formulated as . The idea behind physicsinformed methods is to regularize the data-driven NNs solutions to satisfy known physical relations and constraints. Therefore, physics information is often introduced as soft constraints by modifying the original data-driven loss (Eq. (6)) to include the physics-based loss contributions (Eq. (7)). The resulting physics-informed loss function, or hybrid loss function (), is formulated as a weighted sum of data-driven and physics-based loss functions, and the optimization problem is posed as  𝜽= arg min 𝜃 (𝒙,𝑺;𝜽,𝝀,𝜶), (𝒙,𝑺;𝜽,𝝀,𝜶) = 𝜆𝐷(𝒙,𝑺;𝜽) + 𝜆𝑃𝑃(𝒙,𝑺;𝜽,𝜶), (8) where 𝝀= (𝜆𝐷, 𝜆𝑃)is the weights’ vector with 𝜆𝐷and 𝜆𝑃the loss weights corresponding to data-driven and physics-based loss functions, respectively. Note that for a modeling problem to be well-posed, proper boundary and initial conditions need to be defined and imposed as constraints. Nevertheless, the thermodynamics-informed regression problem aimed at this work does not require the enforcing of PDEs and ICs International Journal of Thermofluids 20 (2023) 100448 4 N. Masclans et al. as commonly done for physics-informed neural networks (PINNs), as the problem is not solved in time but the data is fitted for a certain time instant; see Section 4.1 for details on the thermodynamics-informed loss function of the proposed TINN method. 2.4. Physics-constrained architecture The introduction of IC/BCs within the loss function through a soft enforcement approach (Eq. (7)) may present limitations, including (i) the absence of a quantitative guarantee regarding the accuracy of the imposed IC/BCs, and (ii) difficulties in setting appropriate weights to each loss term, being the optimization problem (Eq. (8)) dependent on their relative importance [67]. In order to overcome these difficulties, the IC/BCs can be alternatively imposed in a hard manner through their incorporation into the network architecture, i.e., physicsconstrained architecture. Contrarily to the loss functions approach, this technique ensures the exact satisfaction of constraints, and thus reduces training computational cost. In this regard, the proposed TINN approach imposes BCs through a thermodynamics-constrained architecture as hard constraints, while the real-gas EoS is introduced into the thermodynamics-informed loss function as a soft constraint; see Section 4.3 for details on the thermodynamics-constrained architecture of the proposed TINN method. 3. High-pressure transcritical turbulence database The present work considers a high-pressure transcritical channel flow problem using direct numerical simulation (DNS) data obtained from the in-house compressible flow solver RHEA [68]. The computational dataset in terms of (i) problem setup, (ii) flow physics modeling, (iii) supercritical flow states, and (iv) thermophysical data is described below. 3.1. Problem setup As illustrated in Fig. 1, a microfluidic high-pressure transcritical turbulent channel flow is chosen to evaluate the performance of the regression methods. The fluid selected is N2, with critical pressure 𝑃𝑐= 3.4MPa, critical temperature 𝑇𝑐= 126.2K, and molecular weight 𝑊= 2.801⋅10−2 kg/mol. The fluid system at supercritical bulk pressure 𝑃𝑏∕𝑃𝑐= 2 is confined between cold (𝑐𝑤) and hot (𝑡𝑤) isothermal walls at 𝑇𝑐𝑤∕𝑇𝑐= 0.75 and 𝑇ℎ𝑤∕𝑇𝑐= 1.5, respectively. These walls are separated at a distance 𝐻= 2𝛿, with 𝛿= 100 μm the channel half-height. The friction Reynolds number at the cold wall is selected as 𝑅𝑒𝜏,𝑐𝑤 =𝜌𝑐𝑤 𝑢𝜏,𝑐𝑤 𝛿∕𝜇𝑐𝑤 = 100, with 𝑢𝜏,𝑐𝑤 =√(𝜇∕𝜌)(𝑑⟨𝑢⟩∕𝑑𝑦) the friction velocity, and 𝜌𝑐𝑤 and 𝜇𝑐𝑤 values obtained from the highpressure thermophysical model described in Section 3.2. The mass flow rate in the streamwise direction is imposed through a body force 𝒇, controlled by a proportional feedback loop with gain 𝑘𝑝= 0.1to reduce the difference between the desired (𝑅𝑒𝜏,𝑐𝑤 = 100) and numerical 𝑅𝑒𝜏,𝑐𝑤 values. This setup consisting on (i) high-pressure conditions over the critical pressure, and (ii) a temperature difference between cold and hot walls force the fluid to undergo a transcritical trajectory by operating within a thermodynamic region across the pseudo-boiling line [3,6]. The imposition of a friction Reynolds number of 𝑅𝑒𝜏,𝑐𝑤 = 100 under such high-pressure and high-temperature conditions results in fullydeveloped turbulent flow at the microscale. In detail, based on DNS data, Bernades et al. [4] characterized the turbulent flow physics of such microconfined channel flow as a result of (i) the thermophysical properties of supercritical fluids upon crossing the pseudo-boiling line, and (ii) the introduction of vorticity by a localized baroclinic torque close to the hot wall. The results obtained suggest that microconfined turbulence can be achieved, with mixing and heat transfer increments of 20× and 100× for the cold and hot walls, respectively, compared to that of equivalent low-pressure systems. The computational domain is 4𝜋𝛿×2𝛿×4∕3𝜋𝛿 in the streamwise (𝑥), wall-normal (𝑦), and spanwise (𝑧) directions, respectively. The streamwise and spanwise boundaries are set periodic, and no-slip conditions are imposed on the horizontal boundaries (𝑥−𝑧planes). The computational grid is uniform in the streamwise and spanwise directions with resolutions in wall units (based on 𝑐𝑤 values) of 𝛥𝑥+≈ 9.8and 𝛥𝑧+≈ 3.3. The grid is stretched toward the walls in the vertical direction with the first grid point at 𝑦+=𝑦𝑢𝜏,𝑐𝑤∕𝜈𝑐𝑤 ≈ 0.1and with sizes in the range 0.2≲ 𝛥𝑦+≲2.3. Thus, based on preliminary studies, the proposed grid arrangement corresponds to a DNS of size 128 ×128 ×128 grid points. The simulation strategy starts from a linear velocity profile with random fluctuations [69], which is advanced in time to reach turbulent steady-state conditions after approximately 5flow-through-time (FTT) units; a FTT is defined as 𝑡FTT =𝐿𝑥∕𝑢𝑏∼𝛿∕𝑢𝜏with 𝑢𝑏the bulk velocity and 𝐿𝑥= 4𝜋𝛿 the channel length. 3.2. Flow physics modeling The framework utilized for studying microconfined supercritical fluids turbulence is described below in terms of (i) equations of fluid motion, (ii) real-gas thermophysical properties, and (iii) numerical method. 3.2.1. Equations of fluid motion The turbulent flow motion of supercritical fluids is described by the following set of the conservation equations of mass, momentum, and total energy 𝜕𝜌 𝜕𝑡 + ∇ ⋅(𝜌𝒖)= 0,(9) 𝜕(𝜌𝒖) 𝜕𝑡 + ∇ ⋅(𝜌𝒖𝒖)= −∇𝑃+ ∇ ⋅𝝉+𝒇,(10) 𝜕(𝜌𝐸) 𝜕𝑡 + ∇ ⋅(𝜌𝒖𝐸)= −∇ ⋅𝒒− ∇ ⋅(𝑃𝒖)+∇⋅(𝝉⋅𝒖)+𝒇⋅𝒖,(11) where 𝜌is the density, 𝒖is the velocity vector, 𝑃is the pressure, 𝝉=𝜇(∇𝒖+ ∇𝒖𝑇)− (2𝜇∕3)(∇ ⋅𝒖)𝑰is the viscous stress tensor with 𝜇the dynamic viscosity and 𝑰the identity matrix, 𝒇represents a body force utilized to drive the flow in the streamwise direction, 𝐸=𝑒+|𝒖|2∕2 is the specific total energy with 𝑒=𝑒(𝑃 , 𝑇 )the specific internal energy as detailed in Jofre & Urzay [2], and 𝒒= −𝜅∇𝑇is the Fourier heat conduction flux with 𝜅the thermal conductivity. 3.2.2. Real-gas thermophysical properties Equations of state (EoS) describe the thermodynamic space of solutions for the state variables pressure 𝑃, temperature 𝑇, and density 𝜌of a fluid. In this regard, the Peng–Robinson real-gas EoS [70] is commonly utilized for high-pressure flow systems, which writes as 𝑃=𝑅𝑢𝑇 𝑣 −𝑏−𝑎 𝑣2+ 2𝑏 𝑣 −𝑏2,(12) where 𝑅𝑢is the universal gas constant, 𝑣 =𝑊∕𝜌is the molar volume, and 𝑎and 𝑏are coefficients accounting for real-gas effects related to attractive forces and finite packing volume, respectively. Their analytical expressions at high pressures are derived from the Peng–Robinson equation of state as detailed, for example, in Jofre & Urzay [2], and are related to system parameters such as the fluid critical temperature 𝑇𝑐, critical pressure 𝑃𝑐, and acentric factor 𝜔. The Peng–Robinson realgas EoS is supplemented with high-pressure departure functions [71] for the specific internal energy and heat capacity at constant pressure 𝑐𝑝. These functions account for the difference between the ideal-gas conditions (low pressure - only temperature dependent) and the supercritical conditions (high pressure). Thereby, the NASA 7-coefficient polynomial [72] is employed to calculate the ideal-gas contributions, whereas high-pressure analytical departure functions are derived from the Peng–Robinson EoS, as described in Jofre & Urzay [2]. Finally, in terms of transport coefficients, the dynamic viscosity 𝜇 and thermal conductivity 𝜅are computed using correlation expressions International Journal of Thermofluids 20 (2023) 100448 5 N. Masclans et al. Fig. 1. Schematic illustration of the microconfined high-pressure transcritical turbulent flow setup, including an instantaneous snapshot of the streamwise velocity for a 𝑥−𝑦slice at 𝑡∕𝑡FTT ≈ 13.5. proposed by Chung et al. [73,74], which are applicable to Newtonian fluids at high pressures. Further comprehensive information is available in dedicated works, like for example Jofre & Urzay [2] and Poling et al. [75]. 3.2.3. Numerical method The equations of fluid motion presented in Section 3.2.1 are numerically solved using a standard semi-discretized procedure. In particular, the spatial discretization uses a second-order central-differencing scheme, while the time-advancement is achieved through a third-order strong-stability preserving (SSP) Runge–Kutta explicit approach [76]. The convective terms are calculated by means of the Kennedy–Gruber– Pirozzoli (KGP) splitting [77,78], which has been previously validated for high-pressure supercritical fluid turbulence [79,80]. This numerical framework (i) ensures conservation of kinetic energy by convection, (ii) guarantees local conservation of mass, momentum, and total energy, and (iii) delivers stable computations without requiring artificial dissipation or stabilization procedures. 3.3. Supercritical fluid states Supercritical fluids operate within high-pressure thermodynamic spaces in which intermolecular forces and finite packing volume effects become important. In this regard, it is important to distinguish between supercritical gas-like and liquid-like fluids: (i) a supercritical liquid-like fluid is one whose density is large, and whose transport coefficients behave similarly to a liquid; whereas (ii) the density of supercritical gas-like fluids is smaller, and their transport coefficients vary similar to gases. These regions are divided by a coexistence line, commonly referred to as Widom or pseudoboiling (pb) line [81,82]. When crossing this line, the fluid undergoes a supercritical pseudoboiling (continuous) transition, which is the supercritical analog to subcritical phase change. Contrarily to subcritical boiling, the supercritical state transition does not occur in phase equilibrium but over a finite temperature interval. Consequently, this study utilizes the three-regime-model proposed by Wang et al. [83], which characterizes the supercritical fluid as in liquid-like (LL), two-phase-like (TPL) and gas-like (GL) regimes. These three regimes interface at an onset of pseudo-boiling temperature 𝑇−and at a termination of pseudo-boiling temperature 𝑇+, i.e., supercritical fluids at 𝑇 < 𝑇 −are considered LL, 𝑇 > 𝑇 +are GL, and 𝑇−< 𝑇 < 𝑇 +are TPL. Note that TPL includes the pseudo-boiling line, therefore it is satisfied that 𝑇−< 𝑇𝑝𝑏 < 𝑇 +with 𝑇𝑝𝑏 the pseudo-boiling temperature. As proposed by several authors [5,84], the maximum of the specific isobaric heat capacity 𝑐𝑝,𝑝𝑏 is used as a marker of the pseudo-boiling line. In this regard, the pseudo-boiling temperature 𝑇𝑝𝑏 and enthalpy ℎ𝑝𝑏 are defined at the values where the heat capacity reaches its maximum value 𝑐𝑝,𝑝𝑏. For the determination of 𝑇−and 𝑇+parameters, Banuti et al. [5] employ the approximation of enthalpy evolution during supercritical isobaric transition as ℎ𝑝𝑏 =𝑐𝑝,𝑝𝑏(𝑇−𝑇𝑝𝑏) + ℎ0,𝑝𝑏. Thereby, 𝑇−and 𝑇+are obtained at the intersection between ℎ𝑝𝑏(𝑇)and the enthalpy asymptotes for the ideal liquid ℎ𝐿(𝑇)and gas ℎ𝑖𝐺(𝑇)reference states, respectively. For the case at hand, the three-regimes-model is used to locally classify the supercritical fluid state with respect to the predicted temperature, with 𝑇−= 123.33 K and 𝑇+= 176.80 K for N2at 𝑃∕𝑃𝑐= 2. As an example, at time instant 𝑡∕𝑡FTT ≈ 13.5all three supercritical regimes are represented in the dataset: (i) LL (15.1%) is found in the cold-wall viscous sublayer and partially in the buffer layer, (ii) GL (0.7%) is located at the hot-wall viscous sublayer, and (iii) TPL (84.2%) is found elsewhere, using volume fraction units. Nevertheless, it is worth noting that, at such time instant, these different supercritical regimes represent the 22.2%, 3.2%, and 74.6% of the dataset points, respectively, due to the grid stretching toward the walls, where the presence of LL and GL supercritical fluid is more significant. Previous works have also explored the distribution of supercritical fluid regimes along channels and the location of the pseudo-boiling region at different bulk pressures. In this regard, Kim et al. [85] performed DNS of transcritical turbulent channel flow with differentially heated walls with R-134a at (i) slightly supercritical pressure (𝑃∕𝑃𝑐= 1.1) and (ii) significantly higher supercritical pressure (𝑃∕𝑃𝑐= 2), as utilized in this work. Their results indicate a direct relationship between the increase in walls heat difference and the displacement of the average pseudoboiling transitioning region towards the hot/top wall. Additionally, higher supercritical pressures result in broader pseudo-boiling transitioning regions, where the supercritical fluid exhibits smoother density and temperature gradients along the channel height upon undergoing pseudo-boiling. In contrast, significantly larger thermodynamic gradients are observed at the pseudo-boiling region near critical pressures [3,85]. Therefore, and as a consequence of the significantly higher supercritical pressure and large wall heat differences employed in this work, as depicted in Fig. 5, a broad pseudo-boiling transitioning region (TPL represents the 84.2% of the channel volume) and gradual change of thermophysical properties along the channel is observed. 3.4. Thermophysical properties data The computational dataset includes multiple temporal snapshots of the main thermophysical fields with spatial resolution of 128 ×128 × 128. These snapshots are distributed along a total simulation time of 13.5 FTTs, with the fluid achieving fully turbulent flow characteristics at 𝑡∕𝑡FTT ≈ 5. As an example, Fig. 2 illustrates the turbulent contours of normalized stream-wise velocity (𝑢+), reduced density (𝜌∕𝜌𝑐) and reduced temperature (𝑇∕𝑇𝑐) at 𝑡∕𝑡FTT ≈ 13.5. The dataset snapshots consist of instantaneous and time-averaged velocity-related fields. From these fields, additional velocity-based quantities can be derived, such as normalized (∗) turbulent kinetic energy (TKE∗), velocity-gradient invariants (𝑃𝐴∗, 𝑄𝐴∗, 𝑅𝐴∗, 𝑃𝑆∗, 𝑄𝑆∗, 𝑅𝑆∗, 𝑃𝑊∗, 𝑄𝑊∗, 𝑅𝑊∗), vorticity (𝜔∗), and enstrophy (∗). Note that these quantities are normalized as 𝑇 𝐾𝐸∗=𝑇 𝐾𝐸∕𝑢2 𝑏,𝑃∗ 𝑖=𝑃𝑖(𝛿∕𝑢𝑏),𝑄∗ 𝑖= International Journal of Thermofluids 20 (2023) 100448 6 N. Masclans et al. Fig. 2. Snapshots of instantaneous streamwise velocity in wall units 𝑢+(top), reduced density (center), and reduced temperature (bottom) for a 𝑥−𝑦slice at 𝑡∕𝑡FTT ≈ 13.5. 𝑄𝑖(𝛿∕𝑢𝑏)2,𝑅∗ 𝑖=𝑅𝑖(𝛿∕𝑢𝑏)3,𝜔∗=𝜔(𝛿∕𝑢𝑏), and ∗=(𝛿∕𝑢𝑏)2, where subscript 𝑖refers to the different tensors 𝐴,𝑆and 𝑊, and 𝑢𝑏is the bulk streamwise velocity. To reduce the high-dimensionality of the dataset, while retaining its essential information, a feature selection process has been performed using Mutual Information (MI) analysis [86]. This approach selects the most informative features for building a reduced dimensionality dataset, leading to improved computational accuracy, redundancy reduction, improved model performance and deeper insights into underlying patterns and relationships within the data. The MI value quantifies the degree of dependency between random variables, with high MI values indicating higher dependency, and 0for independent variables. The MI analysis depicted in Fig. 3 revealed that the thermophysical quantities targeted exhibit the highest dependency on certain features, ranked in order of importance as follows: the non-dimensional channel height (𝑦∕ℎ), the streamwise instantaneous velocity (𝑢), the normalized turbulent kinetic energy (TKE∗), and the normalized second invariants of the symmetric and skew-symmetric velocity-gradient tensor (𝑄𝑆∗and 𝑄𝑊∗). Notably, the rate-of-rotation 𝑊invariants 𝑃𝑊∗and 𝑅𝑊∗are null by definition as 𝑊 is skew-symmetric, and vorticity and enstrophy are linearly dependent to 𝑄𝑊∗by construction. Nevertheless, the computation of invariants involves finite-difference calculations, making it more computationally expensive compared to other readily available velocity-based quantities and possibly introducing computationally errors. Therefore, based on MI results, only three features, namely 𝑢∕𝑢𝑏,𝑦∕𝛿and 𝑇 𝐾𝐸∗, are initially selected. However, to ensure consistency and generalization across various problem configurations and symmetries, the other velocity components 𝑣and 𝑤are subsequently added. Consequently, the reduced feature set considered comprises the normalized distance to the bottom wall (𝑦∕𝛿), the three normalized velocity components (𝑢∕𝑢𝑏, 𝑣∕𝑢𝑏,𝑤∕𝑢𝑏), and the normalized turbulent kinetic energy (TKE∗). 4. Thermodynamics-informed neural network framework As introduced in Section 1, the goal of this work is to develop and implement a TINN regression model for reconstructing thermophysical properties from velocity-based data of microconfined high-pressure transcritical fluid turbulence. The TINN takes spatial and velocityrelated variables as inputs 𝒙= (𝑦, 𝑢, 𝑣, 𝑤, TKE), and predicts the hidden thermophysical properties of interest as outputs 𝒛= (𝜌, 𝑇 ), with all inputs and outputs normalized between 0and 1. Additionally, the predicted density and temperature values are used to calculate other relevant thermophysical properties 𝝓= (𝑐𝑝, 𝑐, 𝜇, 𝜅)to assess the model’s ability to capture the underlying physics of the problem. More details on the calculation of these thermophysical quantities can be found in Section 3.2.2. In particular, as illustrated in Fig. 4, it is proposed a straightforward design with two nodes in the output layer corresponding to the two targeted thermophysical fields. Physical knowledge is introduced into the ML algorithm as soft and hard constraints. On the one hand, the satisfaction of the high-pressure real-gas equation of state (Eq. (12)) is introduced as a soft constraint into the thermodynamicsinformed loss function (Section 4.1), which requires the use of a custom penalty method (Section 4.2). On the other hand, boundary constraints are incorporated as hard constraints through a custom thermodynamics-constrained architecture (Section 4.3). In addition, a localized-physics weighting is added to the physical soft constraint to enhance the performance of the model (Section 4.4). 4.1. Thermodynamics-informed loss function To leverage both the supervised data of thermophysical properties and the physical knowledge of the real-gas equation of state that describes the flow, the thermodynamics-informed optimization problem is formulated as  𝜽= arg min 𝜃 (𝒙,𝑺;𝜽,𝝀), (𝒙,𝑺;𝜽, 𝝀) = 𝜆𝐷𝐷(𝒙;𝜽) + 𝜆𝑃𝑃(𝒙,𝑺;𝜽), (13) where is the hybrid thermodynamics-informed loss function, 𝜆𝐷and 𝜆𝑃are the fixed weights, or penalty coefficients, relative to data-driven and thermodynamics-based loss/penalty terms, respectively, and 𝑺= {𝑊 , 𝜔, 𝑇𝑐, 𝑃𝑐, 𝑃𝑏}is the set of parameters prescribed by the problem setup, namely the molecular weight 𝑊, acentric factors 𝜔, critical temperature 𝑇𝑐and critical pressure 𝑃𝑐of the fluid, and uniform bulk International Journal of Thermofluids 20 (2023) 100448 7 N. Masclans et al. Fig. 3. MI analysis of the dataset between main reduced velocity-based quantities (features) and normalized thermophysical properties (targets). Fig. 4. Schematic depicting the TINN framework for reconstructing thermophysical information from velocity data of microconfined high-pressure transcritical fluid turbulence. pressure 𝑃𝑏. The data-based loss function 𝐷incorporates the mean squared error of the normalized thermophysical targets as 𝐷(𝒙;𝜽) = 𝜌+𝑇=1 𝑁 𝑁 ∑ 𝑖=1 |||𝜌𝑖−𝜌𝑖||| 2+1 𝑁 𝑁 ∑ 𝑖=1 |||𝑇𝑖− 𝑇𝑖||| 2,(14) with 𝜌𝑖=𝜌 (𝑥𝑖;𝜽)and  𝑇𝑖= 𝑇(𝑥𝑖;𝜽)the normalized predicted values, and 𝑁the number of domain samples. The thermodynamics-based loss function 𝑃guides the training to physics-consistent predictions that satisfy the real-gas equation of state (Eq. (12)), which can be expressed as the normalized mean squared error in the real-gas EoS as 𝑃(𝒙,𝑺;𝜽) = 1 𝑁 𝑁 ∑ 𝑖=1 ||||| 𝑃𝑏− 𝑅𝑢 𝑇𝑖 𝑝ℎ𝑦  𝑣𝑖 𝑝ℎ𝑦 −𝑏 +𝑎𝑖 ( 𝑣𝑖 𝑝ℎ𝑦)2+ 2 𝑏 𝑣𝑖 𝑝ℎ𝑦 −𝑏2||||| 2/𝑃2 𝑏,(15) where 𝜌𝑖 𝑝ℎ𝑦 and  𝑇𝑖 𝑝ℎ𝑦 are the predicted values transformed to the original physical scaling,  𝑣𝑖 𝑝ℎ𝑦 = 𝑣𝑝ℎ𝑦(𝑥𝑖,𝑺;𝜽) = 𝑊∕𝜌𝑖 𝑝ℎ𝑦 is the calculated molar volume, and 𝑎𝑖=𝑎(𝑥𝑖, 𝑆;𝜽) = 𝑓( 𝑇𝑖, 𝑊 , 𝜔, 𝑇𝑐, 𝑃𝑐)and 𝑏(𝑆) = 𝑓(𝑊 , 𝑇𝑐, 𝑃𝑐)are the coefficients of the Peng–Robinson real-gas EoS, being  𝑣 and 𝑎 dependent on predicted outputs 𝜌 and  𝑇, respectively. Unlike typical PINNs modeling implementations (Eq. (7)), the proposed TINN framework deviates from including governing PDEs and ICs in the thermodynamics-informed loss function (Eq. (15)). Instead, TINN is utilized as a thermodynamics-informed regression method, where the hybrid loss function incorporates the thermodynamic equation term to enforce compliance with the real-gas EoS. Note that BCs are not explicitly included in the loss function either as they are already incorporated into the network architecture as hard constraints to alleviate the optimization problem and guarantee its satisfaction; further details regarding the incorporation of BCs can be found in Section 4.3. 4.2. Penalty method Optimizing neural networks with hybrid loss functions can pose difficulties when dealing with large penalty coefficients or unbalanced penalty terms. The weighting of the thermodynamics-based and the data-driven loss contributions must be carefully considered, which can be a time-consuming task. Large penalty coefficients can result in an illconditioned optimization problem with slow convergence, while small penalty terms may lead to poor solutions, with the gradient descent optimization converging to suboptimal local minima. In particular, for the current problem, the proposed optimization problem (Eq. (13)) yields constant thermophysical fields throughout data domain when the soft International Journal of Thermofluids 20 (2023) 100448 8 N. Masclans et al. constrain penalty term is excessively large, as shown in the results presented in Section 5.7. To overcome these challenges, the introduction of a penalty method (PM) is required, as commonly performed in PINNs applications for solving forward & inverse problems [59,60,87]. The penalty method considers a sequence of 𝑁𝑃𝑀 unconstrained optimization problems with increasing penalty coefficients, in contrast to the soft constraints optimization approach (8) that defines a constrained optimization problem with fixed penalty coefficients 𝜆𝐷and 𝜆𝑃. The 𝑘th unconstrained optimization problem corresponds to  𝜽𝑘= arg min 𝜃 𝑘(𝒙,𝑺;𝜽,𝝀), 𝑘(𝒙,𝑺;𝜽,𝝀) = 𝜆𝐷𝐷(𝒙,𝑺;𝜽) + 𝜆𝑘 𝑃𝑃(𝒙,𝑺;𝜽), (16) where 𝜆𝐷is fixed and 𝜆𝑘 𝑃is the soft constraint penalty coefficient at 𝑘th iteration of the penalty method. The network parameters are optimized until the loss function converges with tolerance 𝜀, with patience of 𝑁𝑝𝑎𝑡𝑖𝑒𝑛𝑐𝑒 epochs, after which the physics penalty coefficient is increased by a constant factor 𝛽𝑃for the next penalty method iteration as 𝜆𝑘+1 𝑃=𝛽𝑃𝜆𝑘 𝑃,for 𝑘= 0,…, 𝑁𝑃𝑀 − 1.(17) 4.3. Thermodynamics-constrained architecture The problem considered in this work prescribes temperature at coldwall (𝑇𝑐𝑤) and hot-wall (𝑇ℎ𝑤), as described in Section 3.1. In this regard, let Dirichlet BCs for temperature 𝑇be written as 𝑇(𝒙) = 𝑔𝑇(𝑦) = {0.75𝑇𝑐, 𝑦 =𝑦𝑐𝑤 1.5𝑇𝑐, 𝑦 =𝑦ℎ𝑤 , 𝑦 ∈𝛿𝛺𝑦,(18) where the 𝑦-coordinate of the domain is 𝛺𝑦∈ [𝑦𝑐𝑤, 𝑦ℎ𝑤], and with domain boundaries 𝛿𝛺𝑦= {𝑦𝑐𝑤, 𝑦ℎ𝑤}. In addition, an extended linear formulation is proposed for the BCs function 𝑔𝑇as 𝑔𝑇(𝑦) = 𝑇ℎ𝑤 −𝑇𝑐𝑤 𝑦ℎ𝑤 −𝑦𝑐𝑤 (𝑦−𝑦𝑐𝑚) + 𝑇𝑐𝑤, 𝑦 ∈𝛺𝑦.(19) This formulation facilitates the direct and computationally efficient integration of the BCs into the network architecture, as the function is defined and continuous along the entire domain. Then, using the uniform pressure value 𝑃𝑏∕𝑃𝑐= 2 and the real-gas EoS (Eq. (12)), the prescribed density at the wall boundaries (𝜌𝑐𝑤,𝜌ℎ𝑤) is calculated by means of an equivalent equation written in terms of density. The introduction of Dirichlet BCs is achieved by modifying the network architecture, specifically the network output, as depicted in Fig. 4. This approach guarantees the exact satisfaction of BCs by constructing the predicted solutions as  𝒛(𝒙,𝑺;𝜽) = (1 − 𝑙(𝑦)) 𝒈(𝑦) + 𝑙(𝑦)ℒ𝐿(𝒙,𝑺;𝜽),𝒙∈𝛺, (20) where  𝒛= (𝜌,  𝑇)is the predictions vector, 𝒈= (𝑔𝜌, 𝑔𝑇)is the vector of extended boundary conditions, ℒ𝐿(𝒙,𝑺;𝜽)is the network output, 𝛺is the input features vector domain, and 𝑙is a function satisfying {𝑙(𝑦)=0, 𝑦 ∈𝛿𝛺𝑦, 𝑙(𝑦)>0, 𝑦 ∈𝛺𝑦−𝛿𝛺𝑦.(21) The regular geometry of the problem allows to analytically choose the 𝑙-function as [60] 𝑙(𝑦) = [4 𝐿2 𝑦 (𝑦−𝑦𝑐𝑤)(𝑦ℎ𝑤 −𝑦)]𝛼 , 𝑦 ∈𝛺𝑦,(22) which for the normalized 𝑦in the domain 𝛺𝑦∈ [0,1], and with domain boundary 𝛿𝛺𝑦= {0,1} and domain length 𝐿𝑦= 1, takes the coefficient values 𝑦𝑐𝑤 = 0 and 𝑦ℎ𝑤 = 1. Note that the scaling factor 4∕𝐿𝑦is added to scale the function to be of order one inside 𝛺𝑦. Finally, the exponent 𝛼 < 1is added to reduce the BCs influence on the data points close to the walls; the specific 𝛼value chosen is provided in Table 1. 4.4. Localized-physics weighting Convergence to inaccurate solutions has been observed as the penalty coefficient for the real-gas EoS penalty coefficient increased, particularly for data points near the boundaries of the domain. To better understand this phenomenon, it is crucial to consider two factors: (i) as depicted in Fig. 2, the fluid experiences turbulent flow behavior in the central part of the channel, and (ii) temperature (and density) boundary conditions are enforced at the channel walls (Eq. (18)). Accordingly, the central part of the channel presents higher complexity and time-dependency compared to the regions in the vicinity of wall boundaries. Consequently, a physics-informed regression is particularly necessary to be applied from the buffer layers of the walls inwards, while a data-driven approach already provides satisfactory performance for the wall viscous sublayers. To address this issue, a localized-physics weighting approach has been proposed in the thermodynamics-based penalty term which utilizes non-uniform weights for the real-gas EoS penalty loss term as  𝑃(𝒙,𝑺;𝜽) = 𝑝(𝑦)𝑃(𝒙,𝑺;𝜽),(23) where the localized-physics weighting function 𝑝(𝑦)satisfies {𝑝(𝑦)=0, 𝑦 ∈𝛿𝛺𝑦, 𝑝(𝑦)>0, 𝑦 ∈𝛺𝑦−𝛿𝛺𝑦,(24) and the smooth 𝑝-function is defined in terms of sigmoid functions as 𝑝(𝑦) = 1 1 + 𝑒−𝑘(𝑦𝑏𝑤−𝛿𝑦)−1 1 + 𝑒−𝑘(𝑦𝑡𝑤−𝛿𝑦),(25) which for the normalized 𝑦-coordinate data takes values 𝑦𝑐𝑤 = 0 and 𝑦ℎ𝑤 = 1, and the chosen 𝑘and 𝛿𝑦hyper-parameters are specified in Table 1. 5. Results & discussion The performance of the proposed TINN is thoroughly compared against two traditional methods, namely a Random Forest (RF) [88,89] regressor and a Fully-Connected Neural Network (FCNN) [66,90], to demonstrate the enhanced performance of including physical information and constraints into the ML algorithm. The discussion of the results is presented below. 5.1. Datasets For fair unbiased results, all ML methods are trained and tested using the same dataset. To replicate 𝜇-PIV data, each time snapshot is represented by a single plane, specifically the central 𝑥−𝑦slice of the microchannel, consisting of 128 ×128 data points. Additionally, all inputs and outputs are normalized between 0and 1to enhance the methods performance [66]. The evaluation of the proposed FCNN method, including its extrapolation capabilities, uses two study cases of increasing complexity: Case I and Case II. These cases allow for a comprehensive analysis of the method’s performance under different training conditions. Specifically, Case I involves training and test data with a small time difference of less than 1 FTT, while Case II presents a more challenging scenario with a time difference of 8 FTTs between training and test data. This systematic approach enables a thorough assessment of the method’s ability to handle complex flow phenomena and accurately predict beyond the training data. Fig. 5 illustrates the 𝑥-averaged profiles along the wall-normal direction in viscous units 𝑦+of the true values of reduced density (a) and temperature (b). These plots highlight the similarities in the profiles between training and test datasets of Case I, whereas the profiles for Case II exhibit clear distinctions within the buffer and log-wall regions for both cold (cw) and hot (hw) walls. International Journal of Thermofluids 20 (2023) 100448 9 N. Masclans et al. Fig. 5. Averaged profiles of reduced density 𝜌∕𝜌𝑐(a) and temperature 𝑇∕𝑇𝑐(b) along the wall-normal direction in viscous units 𝑦+for the cold (cw) and hot (hw) walls calculated from ground-truth data of train & test datasets for Case I & II. 5.2. Implementation details The choice of hyper-parameters in ML methods is problemdependent [91]. Consequently, the selected values for this study are specified in Table 1. It is important to note that the hyper-parameters chosen fall within the general recommended range and have not undergone extensive optimization specifically tailored to this problem. The goal of this work is to analyze the impact of incorporating physics in a NN-based framework for a general understanding, rather than achieving the most optimal performance on this specific problem. The RF regressor method employs 48 estimators and evaluates the best split quality based on two features and the mean squared error (MSE). For a fair comparison, both the FCNN and the TINN employ a 24-layer 16-neuron architecture, with the latter including an additional layer to incorporate Dirichlet BC (Section 4.3). The NN-based models are trained using back-propagation with the Adam optimization method [92], and initial tests suggest using a learning rate (𝜂) of 0.0001 and batch size of 64 samples for fast convergence and optimal accuracy. This work uses the popular rectified linear unit (ReLU) activation function [66], for which both the inputs and outputs are scaled within a range of [0,1] and the networks weights are initialized by Kaiming- or He-normal initialization [93] as recommended. Additionally, the proposed TINN uses prescribed data-driven penalty coefficient (𝜆𝐷) set to 1, while the penalty method (PM) is applied to the thermodynamics-based penalty coefficient (𝜆𝑘 𝑃). The PM shows good results with an initial penalty coefficient (𝜆0 𝑃)of 0.001 and an amplification factor (𝛽𝑃) of 1.5. The best balance between data-driven and thermodynamics-based contributions is achieved after 35 − 40 PM iterations with tolerance (𝜀) of 0.005 and patience (𝑁𝑝𝑎𝑡𝑖𝑒𝑛𝑐𝑒) of 10 epochs. The thermodynamics-constrained architecture is implemented with exponent 𝛼= 0.1to ensure a smooth transition between the hard-constrained values and the network predictions. The localizedphysics weighting parameters are chosen as 𝑘= 200 and 𝛿𝑦= 0.05. The FCNN is trained with epoch shuffling for 500 epochs (𝑁𝑇), while the TINN undergoes 40 PM iterations (𝑁𝑃 𝑀 ), which corresponds to approximately 400 − 600 training epochs. Each training is repeated 5 times in order to account for the effect of initialization and shuffling randomness. That is to say, due to the random processes involved in the methods initialization and training, a model trained twice using the same training dataset will generate different predictions when evaluated on the test dataset. In particular, the results in Tables 2 and 3represent the average ±standard deviation of the evaluated metrics from 5 training repetitions. 5.3. Evaluation metrics A quantitative comparison of the synthetic dataset is performed to demonstrate the feasibility and applicability of the proposed method Table 1 Main hyper-parameter values for the ML regression models evaluated. ML model Hyper-parameter RF 𝑛𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑜𝑟𝑠 = 48; criterion = MSE; 𝑛𝑓𝑒𝑎𝑡𝑢𝑟𝑒𝑠,𝑚𝑎𝑥 = 2 FCNN optimizer = Adam; 𝜂= 1 ⋅10−4;𝑁𝑇= 800 epochs TINN optimizer = Adam; 𝜂= 1 ⋅10−4; 𝑁𝑃𝑀 = 40 iterations; 𝑁𝑝𝑎𝑡𝑖𝑒𝑛𝑐𝑒 = 10 epochs; 𝜀= 0.05; 𝜆𝐷= 1;𝜆0 𝑃= 0.001;𝛽𝑃= 1.5;𝛼= 0.1;𝑘= 200;𝛿𝑦= 0.05 compared to other non-physically informed methods. This analysis primarily focuses on the reconstruction of the targeted thermophysical properties, namely density 𝜌and temperature 𝑇, which will be referred to as the predicted fields hereafter. To assess the physical consistency of the predictions, additional thermophysical properties are calculated by applying the predicted 𝜌and 𝑇to the real-gas thermophysical relations, including pressure (𝑃), heat capacity (𝑐𝑝), speed of sound (𝑐), dynamic viscosity (𝜇) and thermal conductivity (𝜅); see Section 3.2.2 for details on the thermophysical model. As a result, accurate reconstruction of these calculated properties provide a strong indication of the physical consistency in the predicted 𝜌and 𝑇. Three well-known evaluation metrics are employed to assess and compare the performance of the ML regression methods to accurately reconstruct the predicted and calculated thermophysical fields: (i) the coefficient of determination (R2) (Eq. (26)), (ii) the Normalized Root Mean Squared Error (NRMSE) (Eq. (27)), and (iii) the Mean Absolute Percentage Error (MAPE) (Eq. (28)). The physical coherence of the solutions is further evaluated by means of the bulk values of reconstructed temperature, density and pressure fields by the Relative Percentage Error (RPE) (Eq. (29)). Additionally, the supercritical fluid is locally classified as LL, TPL or GL state by the three-regime-model (Section 3.3) using the predicted temperature  𝑇. This classification task is evaluated by commonly used accuracy and mean-dice metrics. In detail, defining a thermophysical field 𝑧with predicted field 𝑧 and bulk quantity 𝑧𝑏, the regression evaluation metrics are calculated as R2(𝑧) = 1 − ∑𝑛 𝑖=1(𝑧𝑖−𝑧𝑖)2 ∑𝑛 𝑖=1(𝑧𝑖−𝑧)2, 𝑧 =1 𝑛 𝑛 ∑ 𝑖=1 𝑧𝑖,(26) NRMSE(𝑧) = √1 𝑛∑𝑛 𝑖=1(𝑧𝑖−𝑧𝑖)2 𝑧 ,(27) MAPE(𝑧) = 1 𝑛 𝑛 ∑ 𝑖=1 |||| 𝑧𝑖−𝑧𝑖 𝑧𝑖||||× 100%,(28) RPE(𝑧𝑏) = |||𝑧𝑏−𝑧𝑏||| 𝑧𝑏 × 100%.(29) International Journal of Thermofluids 20 (2023) 100448 16 N. Masclans et al. [37] I. Scherl, B. Strom, J.K. Shang, O. Williams, B.L. Polagye, S.L. Brunton, Robust principal component analysis for modal decomposition of corrupt fluid flows, Phys. Rev. Fluids 5 (2020) 054401, http://dx.doi.org/10.1103/PhysRevFluids.5. 054401. [38] K. Taira, M.S. Hemati, S.L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S.T.M. Dawson, C.-A. Yeh, Modal analysis of fluid flows: applications and outlook, AIAA J. 58 (2020) 998–1022, http://dx.doi.org/10.2514/1.J058462. [39] J. Park, H. Choi, Machine-learning-based feedback control for drag reduction in a turbulent channel flow, J. Fluid Mech. 904 (2020) A24, http://dx.doi.org/10. 1017/jfm.2020.690. [40] F. Pino, L. Schena, J. Rabault, M.A. Mendez, Comparative analysis of machine learning methods for active flow control, J. Fluid Mech. 958 (2023) A39, http://dx.doi.org/10.1017/jfm.2023.76. [41] P. Ma, Y. Tian, Z. Pan, B. Ren, D. Manocha, Fluid directed rigid body control using deep reinforcement learning, ACM Trans. Graph. 37 (2018) http://dx.doi. org/10.1145/3197517.3201334. [42] S. Verma, G. Novati, P. Koumoutsakos, Efficient collective swimming by harnessing vortices through deep reinforcement learning, Proc. Natl. Acad. Sci. USA 115 (2018) 5849–5854, http://dx.doi.org/10.1073/pnas.1800923115. [43] S.L. Brunton, B.R. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annu. Rev. Fluid Mech. 52 (2020) 477–508, http://dx.doi.org/10.1146/ annurev-fluid-010719-060214. [44] A. Zhang, Z.C. Lipton, M. Li, A.J. Smola, Dive into deep learning, 2021, arXiv: 2106.11342. [45] M. Dissanayake, N. Phan-Thien, Neural-network-based approximations for solving partial differential equations, Commun. Numer. Methods. Eng. 10 (1994) 195–201, http://dx.doi.org/10.1002/cnm.1640100303. [46] R. González-García, R. Rico-Martínez, I.G. Kevrekidis, Identification of distributed parameter systems: A neural net based approach, Comput. Chem. Eng. 22 (1998) 965–968, http://dx.doi.org/10.1016/S0098-1354(98)00191-4. [47] P.R. Vlachas, W. Byeon, Z.Y. Wan, T.P. Sapsis, P. Koumoutsakos, Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proc. Math. Phys. Eng. Sci. 474 (2018) 20170844, http://dx.doi.org/ 10.1098/rspa.2017.0844. [48] K. Li, J. Kou, W. Zhang, Deep neural network for unsteady aerodynamic and aeroelastic modeling across multiple mach numbers, Nonlinear Dynam. 96 (2019) 1–21, http://dx.doi.org/10.1007/s11071-019-04915-9. [49] N.C. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, P. Riley, Tensor field networks: rotation- and translation-equivariant neural networks for 3D point clouds, 2018, arXiv:1802.08219. [50] H. Frezat, G. Balarac, J.J. Le Sommer, R. Fablet, R. Lguensat, Physical invariance in neural networks for subgrid-scale scalar flux modeling, Phys. Rev. Fluids 6 (2021) 024607, http://dx.doi.org/10.1103/PhysRevFluids.6.024607. [51] R. Wang, R. Walters, R. Yu, Incorporating symmetry into deep dynamics models for improved generalization, 2021, arXiv:2002.03061. [52] J. Ling, R. Jones, J. Templeton, Machine learning strategies for systems with invariance properties, J. Comput. Phys. 318 (2016) 12567, http://dx.doi.org/10. 1016/j.jcp.2016.05.003. [53] M. Raissi, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707, http://dx.doi.org/10.1016/j.jcp.2018.10.045. [54] S. Cai, Z. Wang, S. Wang, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks for heat transfer problems, J. Heat Transfer 143 (2021) 060801, http://dx.doi.org/10.1115/1.4050542. [55] X. Jin, S. Cai, H. Li, G.E. Karniadakis, Nsfnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations, J. Comput. Phys. 426 (2021) 109951, http://dx.doi.org/10.1016/j.jcp. 2020.109951. [56] G. Karniadakis, Y. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physicsinformed machine learning, Nat. Rev. Phys. 3 (2021) 422–440, http://dx.doi. org/10.1038/s42254-021-00314-5. [57] H. Eivazi, T. Mojtaba, P. Schlatter, R. Vinuesa, Physics-informed neural networks for solving Reynolds-averaged Navier–Stokes equations, Phys. Fluids 34 (2022) 075117, http://dx.doi.org/10.1063/5.0095270. [58] H. Wang, Y. Liu, S. Wang, Dense velocity reconstruction from particle image velocimetry/particle tracking velocimetry using a physics-informed neural network, Phys. Fluids 34 (2022) 017116, http://dx.doi.org/10.1063/5.0078143. [59] Z. Mao, A.D. Jagtap, G.E. Karniadakis, Physics-informed neural networks for high-speed flows, Comput. Methods Appl. Mech. Engrg. 360 (2020) 112789, http://dx.doi.org/10.1016/j.cma.2019.112789. [60] L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, S. Johnson, Physics-informed neural networks with hard constraints for inverse design, SIAM J. Sci. Comput. 43 (2021) 1105–1132, http://dx.doi.org/10.1137/21M1397908. [61] G. Pang, L. Lu, G.E. Karniadakis, fPINNs: Fractional physics-informed neural networks, SIAM J. Sci. Comput. 41 (2019) A2603–A2626, http://dx.doi.org/10. 1137/18m1229845. [62] L. Yang, D. Zhang, G.E. Karniadakis, Physics-informed generative adversarial networks for stochastic differential equations, SIAM J. Sci. Comput. 42 (2020) A292–A317, http://dx.doi.org/10.1137/18M1225409. [63] M. Fiore, L. Koloszar, C. Fare, M.A. Mendez, M. Duponcheel, Y. Bartosiewicz, Physics-constrained machine learning for thermal turbulence modelling at low Prandtl numbers, Int. J. Heat Mass Transfer 194 (2022) 122998, http://dx.doi. org/10.1016/j.ijheatmasstransfer.2022.122998. [64] T. Zhang, S. Sun, Thermodynamics-informed neural network (TINN) for phase equilibrium calculations considering capillary pressure, Energies 14 (2021) 7724, http://dx.doi.org/10.3390/en14227724. [65] Q. Hernádez, A. Badías, F. Chinesta, E. Cueto, Thermodynamics-informed neural networks for physically realistic mixed reality, Comput. Methods Appl. Mech. Engrg. 407 (2023) 115912, http://dx.doi.org/10.1016/j.cma.2023.115912. [66] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016. [67] L. Sun, H. Gao, S. Pan, J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112732, http://dx.doi.org/10.1016/j.cma.2019. 112732. [68] L. Jofre, A. Abdellatif, G. Oyarzun, RHEA - an open-source reproducible hybridarchitecture flow solver engineered for academia, J. Open Source Softw. 8 (2023) 4637, http://dx.doi.org/10.21105/joss.04637. [69] K.S. Nelson, O.B. Fringer, Reducing spin-up time for simulations of turbulent channel flow, Phys. Fluids 29 (2017) 105101, http://dx.doi.org/10.1063/1. 4993489. [70] D.Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59–64, http://dx.doi.org/10.1021/i160057a011. [71] W.C. Reynolds, P. Colonna, Thermodynamics: Fundamentals and Engineering Applications, first ed., Cambridge University Press, Cambridge (UK), 2019, http: //dx.doi.org/10.1017/9781139050616. [72] A. Burcat, B. Ruscic, Third Millennium Ideal Gas and Condensed Phase Thermochemical Database for Combustion with Updates from Active Thermochemical Tables, Technical Report, Argonne National Laboratory, 2005, URL: publications. anl.gov/anlpubs/2005/07/53802.pdf. [73] T.H. Chung, L.L. Lee, K.E. Starling, Applications of kinetic gas theories and multiparameter correlation for prediction of dilute gas viscosity and thermal conductivity, Ind. Eng. Chem. Fund. 23 (1984) 8–13, http://dx.doi.org/10.1021/ i100013a002. [74] T.H. Chung, M. Ajlan, L.L. Lee, K.E. Starling, Generalized multiparameter correlation for nonpolar and polar fluid transport properties, Ind. Eng. Chem. Fund. 27 (1988) 671–679, http://dx.doi.org/10.1021/ie00076a024. [75] B.E. Poling, J.M. Prausnitz, J.P. O’Connell, Properties of Gases and Liquids, fifth ed., McGraw Hill, New York (USA), 2001, http://dx.doi.org/10.1021/ja0048634. [76] S. Gottlieb, C.W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev. 43 (2001) 89–112, http://dx.doi.org/10.1137/ S003614450036757X. [77] G. Coppola, F. Capuano, S. Pirozzoli, L. de Luca, Numerically stable formulations of convective terms for turbulent compressible flows, J. Comput. Phys. 382 (2019) 86–104, http://dx.doi.org/10.1016/j.jcp.2019.01.007. [78] G. Coppola, F. Capuano, L. de Luca, Discrete energy-conservation properties in the numerical simulation of the Navier-Stokes equations, Appl. Mech. Rev. 71 (2019) 010803, http://dx.doi.org/10.1115/1.4042820. [79] M. Bernades, F. Capuano, F.X. Trias, L. Jofre, Energy-preserving stable computations of high-pressure supercritical fluids turbulence, in: 8th European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS, 2022, pp. 1–12, http://dx.doi.org/10.23967/eccomas.2022.110. [80] M. Bernades, L. Jofre, F. Capuano, Investigation of a novel numerical scheme for high-pressure supercritical fluids turbulence, in: Proceedings of the Summer Program 2022, Center for Turbulence Research, Stanford University, 2022, pp. 225–234, https://web.stanford.edu/group/ctr/ctrsp22/iii11_Bernades.pdf. [81] K. Okamoto, J. Ota, K. Sakurai, H. Madarame, Transient velocity distributions for the supercritical carbon dioxide forced convection heat transfer, J. Nucl. Sci. Technol. 40 (2003) 763–767, http://dx.doi.org/10.3327/jnst.40.763. [82] M. Oschwald, J.J. Smith, R. Branam, J. Hussong, Injection of fluids into supercritical environments, Combust. Sci. Technol. 178 (2006) 49–100, http: //dx.doi.org/10.1080/00102200500292464. [83] Q. Wang, X. Ma, J. Xu, M. Li, Y. Wang, The three-regime-model for pseudoboiling in supercritical pressure, Int. J. Heat Mass Transfer 181 (2021) 121875, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121875. [84] L. Xu, P. Kumar, S.V. Buldyrev, S.H. Chen, P.H. Poole, F. Sciortino, H.E. Stanley, Relation between the Widom line and the dynamic crossover in systems with a liquid–liquid phase transition, Proc. Natl. Acad. Sci. USA 102 (2005) 16558–16562, http://dx.doi.org/10.1073/pnas.0507870102. [85] K. Kim, J.P. Hickey, C. Scalo, Pseudophase change effects in turbulent channel flow under transcritical temperature conditions, J. Fluid Mech. 871 (2019) 52–91, http://dx.doi.org/10.1017/jfm.2019.292. [86] A. Kraskov, H. Stögbauer, P. Grassberger, Estimating mutual information, Phys. Rev. E 69 (2004) 066138, http://dx.doi.org/10.1103/PhysRevE.69.066138. [87] A.D. Jagtap, E. Kharazmi, G.E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Comput. Methods Appl. Mech. Engrg. 365 (2020) 113028, http://dx.doi.org/10.1016/j.cma.2020.113028. [88] T.K. Ho, Random decision forests, in: Proceedings of 3rd International Conference on Document Analysis and Recognition, Vol. 1, 1995, pp. 278–282, http://dx. doi.org/10.1109/ICDAR.1995.598994. International Journal of Thermofluids 20 (2023) 100448 17 N. Masclans et al. [89] F. Vázquez-Novoa, J. Conejero, C. Tatu, R.M. Badia, Scalable random forest with data-parallel computing, in: 29th International European Conference on Parallel and Distributed Computing, 2023. [90] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (2015) 436–444, http://dx.doi.org/10.1038/nature14539. [91] D.P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods, Academic Press, 2014. [92] D.P. Kingma, J. Ba, Adam: A method for stochastic optimization, 2017, arXiv: 1412.6980. [93] K. He, X. Zhang, S. Ren, J. Sun, Delving deep into rectifiers: surpassing human-level performance on imagenet classification, 2015, arXiv:1502.01852.