DATA-DRIVEN AND NONLOCAL APPROACHES IN MODELING, ANALYSIS AND SIMULATION OF TURBULENT MIXING PHENOMENA By Ali Akhavan Safaei A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy Computational Mathematics, Science and Engineering – Dual Major 2022 ABSTRACT The overreaching goal of this study is utilizing data-driven methods and sophisticated math- ematical tools for modeling and simulation of turbulent transport of passive scalars. We focus on embedding the intrinsic nonlocal nature of the turbulence into our models. We study the nonlocal dynamics in the context of (i) subgrid-scale (SGS) modeling for large- eddy simulation (LES), and (ii) the turbulent cascade under large-scale anisotropic sources. Moreover, we implement stochastic modeling methodologies to systematically investigate the contributing mechanisms leading a high-speed hydrodynamic transport system into insta- bility and chaos, as well as discovering the anomalies in the featured characteristics of the transport. First, we present a computational-statistical framework to obtain high-fidelity data for homogeneous isotropic turbulent (HIT) flow and passive scalar transport. A parallel im- plementation of the well-known pseudo-spectral method in addition to the comprehensive record of the statistical and small-scale quantities of the turbulent transport are offered for executing on distributed memory CPU-based supercomputers. Afterwards, we investigate the inherent nonlocal behavior of the SGS passive scalar flux through studying its two-point statistics obtained from the filtered direct numerical simula- tion (DNS) data for passive scalar transport in HIT flow. We propose a statistical model for microscopic SGS motions by considering the filtered Boltzmann transport equation (FBTE) for passive scalar. In FBTE, we approximate the filtered equilibrium distribution with an α-stable Lévy distribution that incorporates a power-law behavior to resemble the observed nonlocal statistics of SGS scalar flux. Through generic ensemble-averaging of FBTE, we formulate a continuum-level closure model for the SGS scalar flux appearing in terms of a fractional-order Laplacian that is a nonlocal operator. Moreover, we revisit the spectral transfer model for the turbulent intensity in the passive scalar transport (under large-scale anisotropic forcing), and a subsequent modification to the scaling of scalar variance cascade is presented. Accordingly, we obtain a revised scalar transport model using fractional-order Laplacian operator that facilitates the robust inclusion of the nonlocal effects originated from large-scale anisotropy transferred across the multitude of scales in the turbulent cascade. We provide an a priori estimate for the nonlocal model, and examine the model through a new DNS. We conduct a detailed analysis on the evolution of the scalar variance, high-order statistics of scalar gradient, and two-point statistical metrics of the turbulent transport to compare the developed nonlocal model and its standard version. In another study, a deep learning surrogate model in the form of fully connected feed- forward neural networks is developed to predict the SGS scalar flux in the context of large- eddy simulation of turbulent transport. The deep neural network (DNN) model is trained and validated using filtered DNS dataset at P eλ = 240, Sc = 1 that includes the filtered scalar and velocity gradients as input features. Using the transfer learning concept, we generalize the performance of this trained model to turbulent scalar transport regimes with higher P eλ and Sc numbers with a relatively low amount of data and computations. Finally, in stochastic modeling of hydrodynamic transport, we study the flow dynamics inside a high-speed rotating cylinder after introducing strong symmetry-breaking disturbance factors at cylinder wall motion. We perform a statistical analysis on the fluctuating fields characterizing the fingerprints and measures of intense and rapidly evolving non-Gaussian behavior through space and time. Such non-Gaussian statistics essentially emerge and evolve due to an intensified presence of coherent vortical motions initially triggered by the flow instability due to symmetry-breaking rotation of the cylinder. We show that this mechanism causes significant memory effects in the flow so that noticeable anomaly in the time-scaling of enstrophy record is observed in the long run apart from the onset of instability. Copyright by ALI AKHAVAN SAFAEI 2022 To my beloved wife, Afsoun v ACKNOWLEDGEMENTS First of all, I would like to deeply thank my academic advisor, Dr. Mohsen Zayernouri for providing me this great opportunity in FMATH group. I am very grateful for his wonderful support, mentoring, as well as his great academic and research guidance. These key elements have been invaluable resources that immensely helped me to grow as a researcher, always learn more, and tackle challenging obstacles throughout the course of my PhD studies. My sincere gratitude goes to the rest of my committee members, Dr. Michael Murillo, Dr. Junlin Yuan, and Dr. Tong Gao for their very nice and helpful research discussions, comments, and feedback. Moreover, I am grateful to the support, collaborative works, and friendship of all my former and current peers in the FMATH group: Dr. Ehsan Kharazmi, Dr. Mehdi Samiee, Dr. Jorge Suzuki, Dr. Eduardo de Moraes, Ms. Pegah Varghaei, and Dr. Hadi Seyedi. I always cherish our nice memories and friendships in graduate school. I am grateful to the faculty and staff of the Department of Mechanical Engineering (ME) and the Department of Computational Mathematics, Science, and Engineering (CMSE), and especially thank Ms. Jessica Pung from ME department and Ms. Heather Williamson from CMSE department for always being nice and helpful. I am beyond grateful to my family and close friends, especially to my parents Maryam and Morteza for their unconditional love and support. Most importantly, I thank my beloved wife Afsoun, for her true love, emotional support, and exemplary patience that have always been the most valuable assets in my life. She has been my driving force during my years in the PhD program especially to overcome the unimaginable hardships and challenging times. Therefore, I dedicate this dissertation to her. Finally, I would like to acknowledge the support for high-performance computing re- sources and services that were provided by the Institute for Cyber-Enabled Research (ICER) at Michigan State University. vi TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Turbulence and Anomalous Transport . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Data-driven Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Nonlocal Modeling with Fractional-order Operators . . . . . . . . . . . . . . 7 1.4 Outline of the Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 CHAPTER 2 A PARALLEL COMPUTATIONAL-STATISTICAL FRAMEWORK FOR SIMULATION OF TURBULENCE . . . . . . . . . . . . . . . 12 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Solver Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 An Illustrative Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Impact and Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 CHAPTER 3 FRACTIONAL-ORDER SUBGRID-SCALE MODELING FOR LES OF SCALAR TURBULENCE . . . . . . . . . . . . . . . . . . . . . 25 3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 Why SGS Dynamics are Statistically Nonlocal? . . . . . . . . . . . . . . . . 27 3.4 Boltzmann Transport Framework . . . . . . . . . . . . . . . . . . . . . . . . 32 3.5 Data-driven Nonlocal SGS Modeling . . . . . . . . . . . . . . . . . . . . . . 39 3.6 A Priori Testing via SGS Dissipation of the Resolved Scalar Variance . . . . 44 CHAPTER 4 A NONLOCAL SPECTRAL TRANSFER MODEL AND SCALING LAW FOR SCALAR TURBULENCE . . . . . . . . . . . . . . . . . 45 4.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Turbulent transport of passive scalars . . . . . . . . . . . . . . . . . . . . . . 46 4.3 Statistical analysis of the nonlocal scalar turbulence model . . . . . . . . . . 53 4.4 Reconciliation with the fractional-order SGS modeling for LES . . . . . . . . 61 CHAPTER 5 DEEP LEARNING MODELING FOR SGS FLUX IN THE LES OF SCALAR TURBULENCE AND GENERALIZATION TO OTHER TRANSPORT REGIMES . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Modeling Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.3 Numerical and Optimization Methods . . . . . . . . . . . . . . . . . . . . . . 67 5.4 Model Assessments in the Inference Mode . . . . . . . . . . . . . . . . . . . 74 CHAPTER 6 ANOMALOUS TRANSPORT IN INTERNAL CYLINDER FLOW INSTABILITIES SUBJECT TO UNCERTAINTY . . . . . . . . . . . 86 6.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.3 Stochastic Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . 88 6.4 Stochastic Computational Fluid Dynamics Framework . . . . . . . . . . . . 90 6.5 Variance-Based Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . 94 vii 6.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 CHAPTER 7 SUMMARY AND FUTURE WORKS . . . . . . . . . . . . . . . . . 112 7.1 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.2 Future Directions and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . 118 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 APPENDIX A FRACTIONAL-ORDER DIFFERENTIAL OPERATORS . . . . 140 APPENDIX B DERIVATION OF PASSIVE SCALAR FLUX . . . . . . . . . . 141 APPENDIX C VALIDATION OF NUMERICAL SETUP . . . . . . . . . . . . . 143 APPENDIX D COMPUTATIONAL WORKFLOW . . . . . . . . . . . . . . . . 145 viii CHAPTER 1 INTRODUCTION 1.1 Turbulence and Anomalous Transport Understanding, quantifying, and exploiting anomalous transport open up a rich field, which can transform our perspective towards the extraordinary processes in thermo-fluid problems. This emerging class of physical phenomena refers to fascinating and realistic pro- cesses that exhibit non-Markovian (long-range memory) effects, non-Fickian (nonlocal) in- teractions, non-ergodic statistics, and non-equilibrium dynamics [1]. It is observed in a wide variety of complex, multi-scale, and multi-physics systems such as: sub-/super-diffusion in brain, kinetic plasma turbulence, aging of polymers, glassy materials, amorphous semicon- ductors, biological cells, heterogeneous tissues, and disordered media. Of particular interest, the structure of chaotic and turbulent flows is in a way that nonlocal and memory effects cannot be ruled out [2, 3]. In fact, anomalous transport can essentially manifest in heavy-tailed and asymmetric distributions, sharp peaks, jumps, and self-similarities in the time-series data of fluctuating velocity/vorticity fields. On the other hand, large-scale natural flows such as atmospheric ones, as well as a wide variety of engineering applications, are among many systems that are substantially influenced by turbulence. Nonlinearity and stochastisity are two inherent elements of fluid dynamics that when significantly triggered lead the flow into turbulence [4, 5]. Turbulence is characterized by persistent fluctuating field variables that are immensely non-Gaussian and have a multi-scale and ubiquitous influence on the fluid dynamics with a great impact on the quality of transport and mixing (see Figure 1.1. Moreover, notable emergence of the extreme and anomalous events reflected in the statistical measurements of turbulent fields and quantities intensify the level of complexity in turbulent flows [6, 7]. Therefore, taking into account the effects of turbulence cannot be compromised in predictions and design procedure for a fluid system affected by the turbulent regime. Filtering turbulent fields and reduced-order modeling: Although considerable ad- 1 Figure 1.1: Complex nature of gradient of turbulent temperature fluctuations obtained from high-resolution simulation [8]. vancements in the modern computational architectures and high-performance computing (HPC) over the past decade have greatly facilitated the high-fidelity predictions of turbulent transport through direct numerical simulation (DNS), those efforts have mostly remained in the area of canonical and fundamental turbulent transport. Nevertheless, large-eddy simulation (LES) of turbulence has shown a promising path towards robust, accurate, and computationally affordable predictions of the turbulent flow behavior in large-scale and real- world applications [9]. In fact, LES is considered as a reliable trade-off between the DNS and the low-fidelity simulations with Reynolds-Averaged Navier-Stokes (RANS) models. The main idea in the LES is that for sufficiently high-Reynolds flows that the statistics of tur- bulent fluctuations associated with small-scale motions are isotropic and hence we expect a universal behavior, one can numerically resolve the large-scale motions while dealing with the subgrid-scale (SGS) effects through proper closure modeling means that utilize resolved- 2 Figure 1.2: Development of turbulent fluctuations under large-scale anisotropy [19]. scale variables. In practice, there is a spatial filtering acting on the conservation equations of transport that represents the LES equations [10, 11]. A vital point to credibly certify an SGS model for the LES is its capability to accurately encode the statistics of turbulent transport and SGS dynamics [12, 13]. Therefore, developing a statistically consistent LES closure model that is structurally capable of capturing the nonlocal interactions in the turbulent dissipation [14, 15] is a vital, that are intensified in the SGS effects during LES [16, 17, 18]. Large-scale anisotropy effects in turbulence: Understanding the mechanisms respon- sible for transport of a passive scalar, e.g. temperature field, in high-speed turbulent flow medium is of fundamental importance for scientific and engineering applications. For exam- ple, turbulent regime is known to enhance the mixing in passive scalars by the molecular diffusion where is a result of growth of small-scale fluctuations, distortion of scalar inter- faces, and the occurrence of highly intermittent scalar gradients at small scales of transport. 3 Advancement in understanding of such phenomena is closely dependent on unraveling the complexity that is enforced by the strong nonlinear couplings over a vast range of scales that are also accompanied by the stochastic nature of turbulence. Therefore, considerable efforts has been devoted to study the structure of turbulent transport in passive scalars in scale-space description at high Reynolds regime. These efforts were originated from the Kolmogorov’s scale-space description of turbulence [20, 21] that related the statistics of ve- locity increments to the average dissipation rate of the turbulent kinetic energy (TKE). Kolomogorov’s theory was fundamentally constituted based on a local model for turbulent energy cascade as demonstrated in Onsager’s cascade model for turbulent spectra [22, 23]. This theory was later extended to the turbulent transport of the passive scalars by [24], [25], and [26]. Afterwards, the analogy for the different regimes of passive scalar transport given the diffusivity range was developed by [27, 28]. The initial Kolmogorov’s theory developed in 1941 was later refined in order to take into account the strong intermittency in local energy dissipation rate [29, 30]. Similarly, highly intermittent fluctuations in the local energy dissipation and also scalar dissipation rates led to the development of the refined similarity hypotheses for passive scalars as presented in [31]. One of the main pillars of Kolmogorov’s theory and its extension to the passive scalars is the local isotropy at small-scales. In the case of passive scalar turbulence with large-scale anisotropy (e.g. non-zero mean gradients), it has been shown that the statistics of small-scale scalar fluctuations remain anisotropic [32] as [33] showed that for the case of homogeneous shear turbulent flows. On the other hand, passive scalars are known to exhibit anomalies such as the large-scale behavior that cannot be ruled out with the AD equation directly [34]. Nonlocal modeling of the effects of large-scale anisotropy (such as an imposed mean gradient, see Figure 1.2) throughout the cascade process would be a promising approach to construct a predictive computational framework for passive scalars. Symmetry-breaking effects and modeling stochasticity: Flow within and around cylinders is a rich physical problem that involves complex geometry and nonlinear flow 4 instabilities, with unsolved questions on flow/vortex structures and anomalous turbulent mixing [35]. Numerous researchers have studied the flow and heat transfer characteristics when a fluid flow encounters a cylinder. These studies include fixed, cross-flow oscillations, inline oscillations, and rotation of the cylinder cases. Studies related to the interactions of the flow and moving bodies were first conducted by Strouhal in 1878. Gerrard [36] proposed a model for the vortex shedding mechanism and the resulted von Kárámn vortex street. Effects of cross-flow and inline oscillations of a cylinder on vortex shedding frequency were first determined by Koopman [37] and Griffin and Ramberg [38], respectively. These studies are categorized as external flows around cylinders and some significant contributions in this regard may be found in [39, 40, 41, 42, 43]. However, flow inside systems with fast rotation including cylinders, squares and annulus geometries are also of great importance. Turbo- machinery, mixing process, gravity-based separators, geophysical flows, and journal bearing lubrication are all clear examples for these types of internal flows [44, 45]. Moreover, in rotational cylinder flows, the flow may face a concave wall and centrifugal instabilities may be developed when the thickness of boundary layer is comparable to the radius of the curvature. Consequently, centrifugal instabilities lead to formation of stream- wise oriented vortices that commonly called Taylor-Görtler vortices. These vortices can change the flow regime through a transition process to turbulence [46, 47, 48]. In particular, the Taylor problem in Couette flow between two concentric rotating cylinders is another well-known example of centrifugal instabilities in rotating systems, which have been studied experimentally [49, 50, 51] and numerically [52, 53, 54, 55]. In such problems, emergence of the adverse angular momentum is an important mechanisms, which initiates flow instability. More specifically, Lopez et al. [56] studied flow in a fully-filled rotating cylinder, which is driven by the counter-rotation of the endwall and found out that in the presence of considerably large counter-rotation, the separation of the Ekman layer from the endwall generates an unstable free shear layer that separates flow regions against the azimuth velocity. In fact, this shear layer is highly sensitive to the sources of disturbance appearing in the 5 azimuth velocity, which essentially breaks the symmetry in the flow. Other symmetry- breaking effects were further investigated when they are originated from other sources such as inertial waves [57], oscillating sidewalls [58] and, precessional forcing [59]. 1.2 Data-driven Modeling Abundance of data, and public access to modern machine learning (ML) libraries have shaped a data-driven era virtually for every modeling disciplines such as those in physics- based and engineering [60]. Over the past decade, the advent of scalable and computation- ally efficient implementations of deep learning libraries has sparked a growing interest in developing and using a variety of deep neural network (DNN) architectures that constitute predictive models from latent features of datasets even with nonlinear complexities. For instance, learning parametric spaces [61, 62], reduced-order modeling [63], inverse problems [64, 65, 66, 67], solving forward differential equations [68, 69, 65, 67], discovering governing equations [70, 71], and data generation [72, 73, 74] are a few examples among important applications in the physics-based modeling via deep learning. Deep Learning for Fluid Mechanics and Turbulence Modeling Fluid mechanics is among the challenging branches of physics-based modeling that has notably benefited from diverse deep learning architectures, recently [75, 76]. In turbulence modeling, promising advances has been made mainly due to prevalence of high-fidelity simu- lations and experimental data [77, 78, 79, 80, 81, 82]. Many of these contributions were made in the area of Reynolds-Averaged Navier-Stokes (RANS) modeling [see e.g., 83, 84, 85, 86, 87]. Particularly, in the Large-Eddy Simulation of turbulent transport, researchers have tra- ditionally proposed models for the closure terms that appear in the filtered governing equa- tions, and tried to model the effects of closure (unresolved) terms using resolved-scale flow quantities. These modeling strategies are predominantly categorized into “functional” or “structural” approaches [88]. In a functional model, the closure term takes the form of a mathematical operator acting on the resolved-scale field. Therefore, functional model are only reproducing the net transfer of turbulence intensity from the resolved to unresolved 6 scales. However, a structural modeling approach would approximate the LES closure in terms of the resolved-scale field, where the subgrid-scale (SGS) structures and statistical characteristics are recovered from the resolved scale information. However, novel DNN surro- gate modeling approaches have recently opened a thriving “data-driven” direction to closure modeling for LES. The study by Sarghini et al. is among the earliest reported works focus- ing on utilizing a neural network architecture for SGS modeling of turbulence [89]. More recently, a study by Beck et al. [90] managed to test the SGS modeling of closure term appearing in the filtered momentum equation with deep neural networks in the context of decaying homogeneous isotropic turbulent (HIT) flows. For more recent advancement of the LES modeling in the context of HIT flows, the reader is referred to [91, 92, 93, 94]. Moreover, studies such as Gamahara and Hattori’s [95] on modeling SGS stress tensor and LES in tur- bulent wall flows, have employed more advanced implementations of a deep neural network [for more recent works see: 96, 97, 98, 99], while Yang et al. [100] presented a physics-based neural network for wall-modeled LES. Interested readers about constraining and embedding the physical laws to the deep learning LES models are referred to [101, 102, 103]. 1.3 Nonlocal Modeling with Fractional-order Operators Fractional calculus is a powerful mathematical tool to generalize the integer-order calculus to its fractional-order counterpart. This essentially results in a broader class of mathemati- cal models, namely fractional ordinary differential equations (FODEs) and fractional partial differential equations (FPDEs) [104, 105]. Due to the remarkable and diverse applications of the FPDEs, development of high-order numerical methods [106, 107, 108, 109, 110, 111, 112] and data-driven numerical schemes [113], as well as numerical studies on the stochastic fractional PDEs [114, 115] have been an active area of research. Despite many important progresses in developing nonlocal mathematical models through the standard methods, frac- tional calculus appears to elevate the capabilities of standard mathematical tools so that the model can take in account the nonlocal effects such as power-law or logarithm singular kernels. Therefore, due to these intrinsic potentials of the fractional calculus in representing 7 the nonlocal interactions, sharp peaks, and memory effects, the fractional-order derivatives seem to be an effective solution to model physical phenomena such as anomalous diffusion processes and hydrodynamic transport (see e.g., [116]). Applications in Engineering In engineering applications, fractional-order differential operators provide a promising and predictive direction in mathematical modeling of the nonlocal behavior such as me- chanics of materials [117, 118, 119, 120], nonlinear vibration analysis [121], biomechanics [122, 123, 124], subgrig-scale modeling for large-eddy simulation of turbulence [16, 125, 126, 17, 127, 128], modeling the near-wall turbulence [129], and Reynolds-averaged Navier-Stokes modeling for wall-bounded turbulent flows [130, 131]. 1.4 Outline of the Study The overall goal of this research is to address the emergence of nonlocal effects in hydro- dynamic transport, and propose modeling, computational, and statistical approaches that are inherently designed to predict nonlocality and its signatures in extreme flow regimes such as turbulence and chaos. In the following, we briefly introduce these contributions: Chapter 2: Our goal is to offer an extensible computational framework that carries out the high-fidelity simulations of homogeneous isotropic turbulence (HIT) for an incompressible flow and also obtains the transport of a passive scalar (temperature or concentration of species) in such flow [132, 133] while it keeps track of statistical quantities of turbulent transport. Here, we numerically solve the incompressible Navier-Stokes (NS) equations in addition to the advection-diffusion (AD) to sufficiently resolve the fluid velocity and passive scalar fields, respectively. Spatial homogeneity of fluctuating fields makes this problem well- suited for pseudo-spectral implementation of the NS and AD equations based on the Fourier collocation discretization as employed in our work. This computational platform is based upon programming in PYTHON and leveraging Message Passing Interface (MPI) library for parallel implementation. Chapter 3: In the context of LES for turbulent flows, a recent study by Samiee et al. 8 [16] introduced a nonlocal model for the divergence of SGS stress tensor in terms of a frac- tional Laplacian acting on the resolved-scale velocity field. In order to derive such a model, filtered Boltzmann transport equation was considered, where the filtered equilibrium distri- bution is approximated with an α-stable Lévy distribution. Moreover, Di Leoni et al. [125] proposed a nonlocal eddy-viscosity SGS model that employs a fractional gradient operator. Their modeling strategy is based upon the high-fidelity observations of nonlocal two-point correlation between the SGS stress and strain-rate tensors (inspired by the derivation of fil- tered Kárman-Howarth equation), and proposing a proper nonlocal convolution kernel that yields the fractional gradient operator. They sufficiently captured the nonlocal SGS effects through proper fractional orders for different turbulent flows including the anisotropy and inhomogeneity effects. These studies demonstrated that the fractional-order operators are sophisticated candidates for modeling of the SGS stresses in the LES of turbulent flows. Of particular interest, we aim to study the nonlocal SGS modeling for the conserved passive scalars in turbulent flows [132, 133, 34]; thus, we seek to model the SGS scalar flux arising as the closure term in the filtered scalar transport equation. Due to promising po- tential of Boltzmann transport framework to investigate the sources of spatial nonlocality appearing in the SGS dynamics [16], we manage to study the filtered version of the Boltz- mann transport equation for the passive scalars in turbulent flow. Using proper statistical assumptions at the kinetic level, we try to derive a continuum level closure model in terms of fractional-order Laplacian of the resolved scalar concentration. Through a statistical data- driven procedure our model is being calibrated to its optimal form so that it is capable of capturing the nonlocal statistics embedded in the ground-truth data. Chapter 4: We show that using a well-resolved standard DNS data for the transport of a passive scalar with uniform mean-gradient in a moderately high-Reynolds fully turbulent flow, the 3-D scalar spectrum does not precisely obey the k −5/3 scaling, and follows a scaling that is enforced by the large-scale anisotropy. Utilizing the Corrsin’s generalized cascade model [134], we propose that the modification to the local time-scale associated with the 9 eddies of size ℓ in a way to account for the nonlocal interactions. This modification returns a scaling relation that matches the scalar spectrum after parameterization. Subsequently, the total scalar dissipation is revised and an additional term in the form of a fractional Laplacian of the scalar concentration is obtained. The performance of the AD equation that is equipped with this nonlocal term is assessed in a seamless DNS setting. The resulting statistical analysis on the fully-developed turbulent scalar field shows that considering the effects of nonlocal interactions in the mathematical model (AD equation) provides a better a more pronounced prediction of small-scale scalar intermittecny along the direction of large- scale anisotropy, and it provides a consistent scaling for the third-order mixed longitudinal structure function over a wide range of scales. Chapter 5: We develop a DNN surrogate model for the SGS scalar flux, and utilize the gradients of the resolved-scale scalar and velocity fields as input features. We propose a spatio-temporal sampling approach for constructing a diverse training/validation dataset from filtered DNS (FDNS) results, and seamlessly optimize the DNN model for a specified turbulent transport regime. Using the transfer learning concept, we generalize our pre- trained DNN model to accurately predict on different scalar transport regimes with higher Schmidt and higher Péclet numbers. We demonstrate that compared to the pre-trained model, the amount of required data and the computational cost of model optimization are significantly reduced in the re-training steps. We conduct multiple tests to evaluate the performance of optimized model in the inference mode on full-size 3-D FDNS grids at multiple time instances. In a priori analysis, we show that the base DNN model decreases the time-averaged spatial error in the predicted SGS flux and the SGS dissipation of scalar variance compared to two traditional models. Similar comparisons are presented between the performance of the pre-trained and transfer learned models for higher Schmidt and higher Péclet number scalar transport regimes. Finally, we test the performance of our base DNN model when utilized in LES setting, and we show that unlike two other conventional SGS models, the DNN is successful in maintaining the accuracy for predicted time-record resolved- 10 scale scalar variance, and time-averaged two-point structure functions of LES-resolved scalar concentration field. Chapter 6: We seek to fill a gap in the rich literature of investigating flow instabilities in- side rotating flow systems by emphasizing on the stochastic modeling of the fluid dynamics and later focusing on the anomalous transport features of such system through statistical and scaling analysis of the response. This goal is achieved through a comprehensive com- putational framework that employs high-fidelity flow simulator as “forward solver” in our stochastic model. As a result, the main contributions of our study are highlighted in the following items: • We formulate stochastic Navier-Stokes equations subject to random symmetry-breaking inputs, affecting the incompressible flow within a high-speed rotating cylinder. We employ spectral element method (SEM) along with the probabilistic collocation method (PCM) to formulate a stochastic computational framework. • We perform a global sensitivity analysis and reduce the dimension of random space to the dominant stochastic directions. We compute the expected velocity field enabling us to obtain the fluctuating part of the velocity at the onset of flow instabilities induced by the modeled symmetry-breaking effects. Computing the velocity fluctuations lets us study the temporal evolution of their probability distribution function, which sheds light on the instability dynamics and anomalous transport features. • Obtaining the fluctuating vorticity field, we identify a well-pronounced and evolving non- Gaussian statistical behavior at the onset of flow instability essentially implying that the disturbances (influencing the cylinder rotation) cause generation of “coherent vortical struc- tures”. These vortices increase the memory effects in the hydrodynamics and we charac- terize their impact as long-time “anomalous” time-scaling of enstrophy leading to effective enhancement in the mixing capacity of the system. Chapter 7: The dissertation is concluded with a discussion over findings and possible future directions to extend the current investigations. 11 CHAPTER 2 A PARALLEL COMPUTATIONAL-STATISTICAL FRAMEWORK FOR SIMULATION OF TURBULENCE 2.1 Background Understanding the complex and random nature of turbulent flows, mixing, and trans- port is a vital step in predictions and the design of systems interacting with such hetero- geneous medium. Turbulence is inherently consisted of multi-scale processes that requires high-accurate measurements at the smallest scales of transport [135, 4]. Direct numer- ical simulation (DNS) of turbulent transport as a rigorous scientific tool is supposed to fully resolve the smallest scales of the motion resulted from the fluctuating fields in spa- tial domain while maintaining a high-order temporal accuracy as turbulence evolves in time [136]. Therefore, developing an open-source, sustainable, portably parallel, and integrated computational-statistical framework with high-order spatial and temporal accuracy provides a useful academic ground for better understanding of complex standard to anomalous tur- bulent transport across a multitude of scales. Moreover, from the educational point of view, developing such a user-friendly scientific software will essentially fill the existing training gap in the subjective trinity, i.e., fluid mechanics, computational fluid dynamics, and turbulent transport; hence, leading to a more cohesive ramp up in training the future generation of researchers in a variety of academic-to-industrial disciplines. Among the current open-source computational frameworks, Nektar++ [137, 138] (the spectral/hp element method flow software), HERCULIS [139] and Xcompact3D [140] (the high- order finite difference flow solvers), GRINS [141] (the adaptive mesh refinement finite element method software), spectralDNS [142] (the spectral method computational package for DNS), and OpenFOAM [143] are the notable contributions to the DNS of turbulent transport. On the other hand, the random nature of turbulence requires a thorough statistical analysis on the fluctuating fields and their gradients so that one can identify when the realistic and fully-developed turbulent state is obtained in DNS during an ongoing simulation. This ne- 12 (a) (b) Figure 2.1: Schematic of the architecture of the software. Fig. (a) Illustrates the pseudo- spectral NS solver to archive fully-developed turbulent state, statistical records, velocity output and restarting the simulation from output file. Fig. (b) Shows the pseudo-spectral NS and AD solvers to reach fully-developed turbulent scalar state. cessitates development of a comprehensive computational platform that includes computing and recording of such statistical quantities of turbulent transport as time-series format. The rest of this chapter is organized as follows: in section 2.2, we describe the details and capabilities of the developed platform as a scientific software and we point out the theoretical backgrounds briefly. Furthermore, in section 2.3, we go over a comprehensive example illustrating the results of a fully-developed turbulent flow and passive scalar field with proper statistical testing and verification. In section 2.4, we outline broader impacts of the current work onto the research in turbulent transport. 2.2 Solver Description Governing Equations The incompressible HIT considered in the present software is governed by the NS equa- tions ∂U + U · ∇U = −∇p + ν ∆U + A U , (2.1) ∂t 13 subject to the continuity ∇ · U = 0. In (6.1), U = (U1 , U2 , U3 ) and p are the instantaneous velocity and modified pressure (pressure divided by the constant density of fluid) fields in the Cartesian coordinate system x = (x1 , x2 , x3 ), respectively. Moreover, ν denotes the dynamic viscosity of the Newtonian fluid, and A is a dynamically evaluated coefficient corresponding to the artificial forcing scheme we employ in order to obtain statistically stationary and fully-turbulent state. From the Reynolds decomposition of instantaneous velocity field, U (x, t) = ⟨U (x, t)⟩ + u(x, t), where ⟨·⟩ represents the ensemble-averaging operator, and u(x, t) denotes the fluctuating part of the velocity field. In HIT, ⟨U (x, t)⟩ = 0; therefore, the instantaneous velocity field equals the fluctuating part that is governed by (6.1). Introducing a passive scalar Φ(x, t) transported in the considered fully-developed HIT flow, the AD equation governing the passive scalar concentration may be formulated as ∂Φ + u · ∇Φ = D ∆Φ, (2.2) ∂t where D denotes the diffusivity of passive scalar. Applying the Reynolds decomposition on the total passive scalar, Φ = ⟨Φ⟩+ϕ, and ϕ is the fluctuating part of the scalar concentration. Considering a uniform mean gradient for the passive scalar as ∇⟨Φ⟩ = (0, β, 0), where β is a constant, the AD equation in (2.2) is rewritten as ∂ϕ + u · ∇ϕ = −β u2 + D ∆ϕ. (2.3) ∂t Fourier Pseudo-Spectral Method Here, we consider spatial homogeneity for the fluctuating velocity and scalar concentra- tion, which allows periodic boundary conditions for these fluctuating fields as u(x + L ei , t) = u(x, t), ϕ(x + L ei , t) = ϕ(x, t), (2.4) where ei , i=1,2,3 , is the unit vector for the i-th direction of the Cartesian coordinate, and L is the periodicity length that defines the spatial domain as Ω = [0, L]3 . Discretizing Ω using a uniform three-dimensional grid returns N 3 grid points with grid spacing along each direction as ∆x = L/N . Transforming this discretization into spectral domain let us have a standard 14 pseudo-spectral representation of the governing equations (6.1) and (2.3). Subsequently, k = (k1 , k2 , k3 ) represents the coordinate system in the spectral space and using Fourier collocation method the discretized representation of k would be ki = (−N/2 + 1, . . . , N/2), i = 1, 2, 3. Accordingly, the discrete Fourier transform of any field variable such as ϕ(x, t) is written as 1 X ϕ(x, t) = 3 ϕ̂k (t) eik·x , (2.5) N k √ where i = −1, and eik·x are the Fourier basis functions. Subsequently, the Fourier coef- P ficients associated with k are represented as ϕ̂k (t) = x ϕ(x, t) e−ik·x . Standard pseudo- spectral formulation of the NS equations based upon the Fourier collocation method is ob- tained after taking the Fourier transform of (6.1), dûk + (u · ∇u)k = −ik p̂k − ν|k|2 ûk + Aûk ; (2.6) dt ik · ûk = 0, By taking the divergence of momentum equation in (2.6) and applying the continuity, modified pressure is explicitly represented in terms of the velocity field. Considering that k · k = |k|2 , one can derive p̂k = ik · (u · ∇u)k /|k|2 ; hence, (2.6) may be reformulated as dûk k · (u · ∇u)k + (u · ∇u)k = k 2 − ν|k|2 ûk + Aûk . (2.7) dt |k| Similarly, the pseudo-spectral representation of the AD equation for passive scalar (2.3), is written as dϕ̂k + (u · ∇ϕ)k = −β ûk2 − D|k|2 ϕ̂k . (2.8) dt Employing the fourth-order Runge-Kutta (RK4) scheme, the time-stepping for both NS and AD equations is explicitly done since the nonlinear (advective) terms are evaluated in the physical space and then transformed into the spectral space. 15 Table 2.1: Statistical characteristics of turbulent flow to be recorded from the NS solver as time-series within user-defied time-intervals. TKE K = 12 ⟨u · u⟩ Dissipation ε = 2 ν ⟨S : S⟩ 1/4 Kolmogorov length-scale η = ν 3 /ε Taylor-scale Reynolds Reλ = λ urms /ν Large-eddy turnover time Te = lo /urms   p p ∗ S = 12 ∇u + ∇uT , urms = 2K/3, λ = urms 15ν/ε, lo = u3rms /ε. (a) (b) 300 0.2 Reλ 250 0.0 Taylor-scale Re VGT Skewness −0.2 200 −0.4 150 −0.6 100 −0.8 ∂u1/∂x1 ∂u2/∂x2 ∂u3/∂x3 50 −1.0 0 5 10 15 20 25 30 35 40 t/Te (c) Figure 2.2: (a) Snapshot of fully-developed turbulent velocity field, u1 component. (b) Time-averaged TKE spectrum. (c) Time-series of Reλ (red dashed line), and VGT skewness factors, Su1,1 , Su2,2 , Su3,3 . Programming Architecture The structure of the prepared software is schematically illustrated in Figure 2.1. Ac- cording to Figure 2.1a, a user starts from a pre-processing step, where the isotropic and 16 random velocity initial condition (IC) is constructed based on a prescribed spectrum for turbulent kinetic energy (TKE). The procedure is the straightforward implementation of the well-known work by Rogallo to generate divergence-free isotropic velocity state [144]. According to Lamorgese et al. [145], the initial TKE spectrum is chosen to be    2 u2rms (κ/kF ) , if κ ≤ kF , E(κ, 0) = × (2.9) kF   (κ/kF )−5/3 , if κ > kF . where κ represents the wavenumber associated with spherical shells, kF denotes the maxi- mum wavenumber of TKE shell we apply artificial forcing to, and urms specifies the initial root-mean-square (rms) intensity of velocity fluctuations. In construction of velocity IC, urms is set to be unity while kF and the number of Fourier collocation points, N are taken as input parameters. In a UNIX/LINUX environment, these inputs are taken as arguments in the execution commandline that are imported through sys library in PYTHON. Once the velocity IC is obtained, it is partitioned into Np slabs according to the slab decomposition method. We completely adopted the implementation of Mortensen and Langtangen [142] for domain decomposition in addition to the parallel implementation of forward and inverse three-dimensional fast Fourier transform (FFT) in PYTHON programming language. Here, the MPI communications depend on the mpi4py library [146, 147, 148]. Having the partitioned velocity IC prepared in the pre-processing step, it is fed to the main body of the software where the initial velocity field might be magnified by a user-defined input argument so that a target TKE is considered for the simulation. The viscosity of fluid, ν, is also is taken as another user-defined input argument. Next, the magnified velocity field is passed into the solver where ûk and (∇u)k = ik ûk are separately transformed back into the physical space so that u · ∇u is simply computed and then transformed into Fourier space. The aliasing error that appears due to this procedure is removed by phase-shifting √ and truncation according to 2 2N/3 as the maximum wavenumber [149]. Afterwards, all of the terms in equation (2.7) are directly evaluated in every stage of RK4 time-integration; however, the last term in the right-hand side of (2.7) is only evaluated after the last stage 17 103 30 102 −2 hφ u2i/hχi 20 Kφ,2 101 10 100 0 10−1 0 5 10 15 20 25 t/Te (a) (b) Figure 2.3: (a) Time-records of production over dissipation of scalar variance (blue dashed line), and the flatness factor for the scalar gradient vector component along the direction of mean scalar gradient (green solid line). (b) Snapshot of fully-developed turbulent passive scalar field. during the artificial forcing by keeping the energy of the low wavenumbers constant, which is associated with the sphere of 0 < |k| ≤ kF . In this procedure, A is computed in a way that the the dissipated energy of turbulent motion is injected to the large-scales. This scheme prevents the flow to undergo decay process before the realistic and fully-developed turbulent state is achieved, nevertheless, the artificial forcing scheme could be turned off through a user-defined input argument if one seek to obtain decaying HIT data. The forcing coefficient could be determined either deterministically [150] or stochastically [151, 152] and both of these methods are supported in the software as would be specified as an input option. Moreover, regarding the stable time-integration, the Courant-Friedrichs-Lewy (CFL) number is dynamically checked through a user-defined time-frequency. According to Eswaran and Pope’s work [151], CFL for this problem is demonstrated as ∆t CFL := max (|u1 | + |u2 | + |u3 |) . (2.10) ∆x In (2.10), ∆x is the uniform grid spacing in each direction and ∆t is the user-defined constant time-interval used in the RK4 time-stepping. In practice, CFL is required to be less than unity to ensure a stable time-integration. Since the fully-developed turbulent state is characterized by a meticulous tracking of 18 statistical quantities of the flow, the present software provides a comprehensive framework for computing and recording the statistical quantities of turbulent flow. Given the homo- geneity of the fluctuating fields, spatial averaging is employed for computing these records at user-defined time intervals. These statistical quantities are categorized into turbulent characteristics of small-scale motion reported in Table 2.1, and high-order central moments of diagonal components of velocity gradient tensor (VGT), ∇u. For instance,  ,  ∂u1 3 ∂u1 2 3/2 Su1,1 = , (2.11) ∂x1 ∂x1  ,  ∂u1 4 ∂u1 2 2 Ku1,1 = , (2.12) ∂x1 ∂x1 where Su1,1 and Ku1,1 indicate the skewness factor and flatness factor (or kurtosis) associated with the first diagonal component of VGT, u1,1 = ∂u1 /∂x1 , respectively. Fully-turbulent flow state would be identified when the time-series of these records reach to a statistically stationary state after long enough time-integration, i.e. approximately 10 to 15 large-eddy turnover times (see Table 2.1). The parallel implementation for computing and collecting these statistical quantities and later recording them in time as different time series which were performed by point-to-point and collective MPI directives. Furthermore, the velocity and pressure fields might be written as output files stored in directories named Out ∗ based on a user-defined time-interval that might be useful for any post-processing after the flow reaches to fully-developed turbulent state. A “restart from file” capability is also designated so that once the statistical record is written out on file, the latest state of velocity field and its related time-integration information are also output on files, which are stored in a directory named Restart. Starting a simulation from either a prescribed IC or restarting it to continue an ongoing simulation that was stopped is specified by a user-defined input argument. All the parallel I/O to store the velocity and pressure fields is done by employing scipy.io library and using loadmat() and savemat() routines for the partition of data that is resolved inside each MPI process in the compressed format and with the machine precision accuracy. This is a fast and efficient I/O while it maintains 19 the simplicity for user for any post-processing step such as time-averaging on the statistically stationary turbulent data. According to Figure 2.1b, once the fully-turbulent velocity state is achieved, the user might be able to use a restart or output instance as the velocity IC to introduce a passive scalar transport with a directional constant mean-gradient as described in (2.3) while the fluctuating concentration is assumed to be zero, ϕ0 (x) = 0. Here, the goal is to resolve the fluctuating scalar concentration field, equation (2.8), transported on the fully-turbulent incompressible flow for long enough time-span so that the fully-developed and realistic tur- bulent state for the passive scalar is obtained. Subsequently, similar procedure as what de- scribed based on the schematic Figure 2.1a is followed while the pseudo-spectral AD solver is fed by the resolved velocity field from the NS solver. The diffusivity of passive scalar, D, is specified by a user-defined input argument for Schmidt number, Sc = ν/D. Accordingly, similar to the NS solver, the advective scalar flux, (u · ∇ϕ)k , is computed in the physical space by the inverse FFT of ûk and (∇ϕ)k and forward FFT computation of u·∇ϕ. Similar dealiasing procedure as described for NS solver is employed in pseudo-spectral AD solver and the RK4 time-integration scheme is utilized to numerically perform explicit time-stepping. In homogeneous scalar turbulence, time evolution of the scalar variance ⟨ϕ2 ⟩ is governed by d 2 ⟨ϕ ⟩ = − 2 ⟨q⟩ · ∇⟨Φ⟩ − ⟨χ⟩, (2.13) dt where q = ϕ u denotes the scalar flux vector, and the scalar dissipation is defined as χ = 2 D ∇ϕ · ∇ϕ [153]. According to (2.3), the first term in the right-hand side of (2.13) is simplified to −2 β ⟨ϕ u2 ⟩ that denotes the scalar variance production (by uniform mean scalar gradient, β). The present software is capable of computing and recording of rate of scalar variance in addition to the production and dissipation terms. This is useful in terms of the checking if the balance of both sides of equation (2.13) holds throughout a simulation so that one ensures that the implementation of the solver works seamlessly. On the other hand, as a measure to evaluate that the statistically stationary state for the passive scalar is 20 achieved is to check if −2 β ⟨ϕ u2 ⟩/⟨χ⟩ ∼ 1 throughout the simulation. Moreover, recording the skewness and flatness factors for the components of fluctuating scalar gradient vector (e.g., Sϕ,2 and Kϕ,2 for ϕ,2 = ∂ϕ/∂x2 similar to 2.11 and 2.12 for VGT) is another statistical indicator measure for fully-developed turbulent passive scalar state. Therefore, in the current computational platform, the user would be able to recognize the statistically stationary state through monitoring the explained time-series data that is written out according to the user- defined time-interval as a software input. The field data output and restart capability for the AD solver is designated similar to the described strategy for the NS solver so that the user would be able to resume an interrupted/stopped simulation and use the output field data for desired applications or post-processing. In the following section, we present a comprehensive example that step-by-step walks through using the present software. 2.3 An Illustrative Example This comprehensive example is mainly consisted of construction of isotropic velocity IC, obtaining well-resolved fully-turbulent velocity field, and simulating well-resolved passive scalar turbulence with imposed mean scalar gradient. IC construction and DNS of the HIT According to the descriptions in section 2.2, the isotropic and divergence-free velocity IC is constructed based upon a prescribed energy spectrum given in (2.9). Considering the periodicity length L = 2π, this pre-processing step is done through serial execution of Gen IC.py script that takes the following input arguments, respectively: N (spatial resolu- tion along each direction), kF (forcing wavenumber), and Np (number of slab partitions). We need to emphasize that Np must be chosen in a way that N be a multiple of Np . The resulting velocity field is located in a directory named IC, where Np number of .mat velocity files are stored. In this example, we take N = 520, kF = 2, and Np = 40. All the compo- nents of velocity IC in addition to the VGT components have Gaussian distribution. This 21 Table 2.2: Input arguments for the PScHIT.py and the specified values for the example case. The order of the arguments in the execution command-line are as listed here. Input Argument Value t end 40 Output frequency 1000 Stats frequency 100 TKE magnification 6.0 ν 0.0008 kF 2 forcing type deterministic N 520 ∆t 0.0005 If Restart 0 or 1 velocity IC is being passed into the NS solver written in PScHIT.py script that takes the fol- lowing input arguments given in Table 2.2. Here, Output frequency and Stats frequency are multiplied by the specified ∆t. Moreover, If Restart argument could be either 0 or 1, where 0 indicates it is a simulation starting from the constructed IC while 1 specifies resum- ing a simulation from restart files. In this example, we perform the simulation for t/Te ∼ 15 to ensure the fully-turbulent flow state is achieved and Figure 2.2a portrays the first com- ponent of the velocity field. Figure 2.2b shows the radial TKE spectrum averaged over 5 large-eddy turnover times. Moreover, Figure 2.2c includes the time-records of the Taylor- scale Reynolds number and VGT skewness factor for diagonal components computed and recorded over 40 large-eddy turnover times. This shows that the statistically stationary state is achieved through the long-time DNS where Reλ ∼ 240, Su1,1 = Su2,2 = Su3,3 ∼ −0.55, and Ku1,1 = Ku2,2 = Ku3,3 ∼ 6.8 at fully-turbulent state. The statistical records of VGT clearly show that the resolved velocity field is isotropic. Finally, kmax η > 1.45 ensures that √ the small scale turbulent motions are well-resolved (kmax = 2N/3) [154]. DNS of passive scalar transport Similar to starting the NS solver from a prescribed velocity IC, we take a fully-turbulent velocity output (velocity state at t/Te = 15 in section 2.3) and continue the simulation under 22 the artificial forcing while we introduce a passive scalar field where its fluctuating part is initialized at zero. The Schmidt number, Sc, is specified by user through an input argument. According to the problem setting for the mean scalar gradient, we let β = 1 (mean scalar gradient along x2 direction). Therefore, for a passive scalar with Sc = 1, we aim to obtain the fully-turbulent scalar field. We need to note that the spatial resolution required for the passive scalars with Sc ≥ 1 is defined based on ηB = η Sc−1/2 [155] and in this example the spatial resolution for the velocity field is sufficient for a well-resolved passive scalar. We manage to resolve the passive scalar field for 25 large-eddy turnover times and the rest of the simulation parameters remain the same as values reported in the Table 2.2. Figure 2.3 shows the records of scalar variance production over dissipation rate, −2 ⟨ϕ u2 ⟩/⟨χ⟩, and the flatness factor for the scalar gradient along the direction of mean scalar gradient, Kϕ,2 . As it is observed, after resolving the passive scalar field for approximately two large-eddy turnover times, −2 ⟨ϕ u2 ⟩/⟨χ⟩ ∼ 1.0 that means the equilibrium state for the passive scalar variance is obtained. Moreover, after approximately three large-eddy turnover times the high-order statistical moments of the scalar gradient reach to a statically stationary state. For instance, Sϕ,2 ∼ 1.4, and Kϕ,2 ∼ 20.8 throughout the time-averaging of these statistical moments when t/Te ≥ 5. By resolving the passive scalar field through AD equation and for long enough time after the equilibrium and stationary state, the fully-turbulent and realistic scalar field is ensured. 2.4 Impact and Applications Current work offers a framework to obtain highly accurate spatio-temporal data for ho- mogeneous turbulent transport with proper statistical testing from the recorded quantities. In turbulent transport research, this provides a great source of high-fidelity data for a va- riety of innovative contributions. For instance in large-eddy simulation (LES) of turbulent transport, novel nonlocal models for the subgrid-scale (SGS) stress and flux terms appearing in the filtered NS and AD equations heavily depend on the DNS data for such transport phenomena to compute the exact values of SGS terms to evaluate the model performance 23 [126, 17, 127, 128, 16]. On the other hand, in the abundance of data and emergence of the data-driven turbulence models [78], current computational platform would be a reliable candidate to generate data for training and testing such models [156, 90, 157, 158, 159, 101]. Moreover, high-Reynolds and well-examined high-fidelity turbulent transport data from the present DNS framework could be directly employed in studying the role of coherent turbulent structures and their effects on the turbulence statistics [160, 161], as well as inves- tigating topological characteristics of turbulent transport [162, 163], and analysis of extreme events as well as the internal intermittency [7, 164, 6]. 24 CHAPTER 3 FRACTIONAL-ORDER SUBGRID-SCALE MODELING FOR LES OF SCALAR TURBULENCE 3.1 Background Traditionally, SGS modeling is categorized into two main branches: (i ) functional model- ing, and (ii ) structural modeling [88]. Functional modeling requires a prior knowledge of the interactions between resolved-scale and subgrid-scale is required so that one can represent the LES closure in terms of a mathematical function of resolved transport variables. Functional models are usually representing the net transfer of turbulent kinetic energy from resolved scales to the subgrid scales. The Smagorinsky model initially conceptualized in [165], and its variations are well-known examples of functional SGS modeling. On the other hand, struc- tural models seek to reconstruct the statics and structure of SGS stresses and fluxes from the resolved-scale variables. For instance, scale-similarity models initially introduced by Bardina et al. [166] are among well-known examples of structural models. Functional models usually are poorly correlated with the true SGS terms a priori and by construction are incapable of reproducing backward transfer of energy (backscattering); however, in an LES setting they have shown to be dissipative enough for solver stability. In contrast, structural models such as scale-similarity type models have been found to be sufficiently correlated with the true SGS terms and fairly capable of following backscattering phenomenon in an a priori sense. Nonetheless, their significant drawback is that in LES they are under-dissipative; hence, the stable time-integration is intractable. As a practical remedy to the mentioned issues, further efforts have been devoted to formulating a mixed representation of functional and structural models [167, 168]. Recently, abundance of the high-fidelity data for the SGS closures mainly available through filtered DNS data, and with the advent of modern machine learning (ML) techniques and their application to fluid mechanics and in particular, turbulence model- ing, [75, 78, 76, 169] have resulted in a wide variety of predictive data-driven SGS models. Among the numerous contributions in ML-based SGS modeling and LES, interested readers 25 are referred to the following notable works [90, 170, 157, 171]. The structure of the rest of this work is organized as follows: in section 3.2, we state the problem and show the governing equations. In section 3.3, we motivate the necessity of our modeling strategy to address nonlocality using statistical measures obtained from the filtered DNS data. In section 3.4, the mathematical framework of our SGS modeling that includes fractional calculus and Boltzmann transport is described and derivation of the SGS model is presented. Afterwards, in section 3.5, a two-stage data-driven calibration procedure is introduced to optimize the model performance. Finally, section 3.6 delivers an a priori testing on the SGS dissipation of the resolved-scale scalar variance. 3.2 Governing Equations Considering flows governed by incompressible Navier-Stokes (NS) equations ∂ui ∂  1 ∂p ∂ 2 ui + ui uj = − +ν + A ui , i, j = 1, 2, 3, (3.1) ∂t ∂xj ρ ∂xi ∂xi ∂xj subject to the continuity, ∇·u = 0, where the velocity and the pressure fields are denoted by u(x, t) = (u1 , u2 , u3 ) and p(x, t) for x = xi and i = 1, 2, 3, respectively. ρ specifies the den- sity and ν represents the kinematic viscosity for a Newtonian fluid. In (3.1), A is a dynamic coefficient associated with the artificial forcing scheme to enforce statistical stationary state on the kinetic energy to reach to a realistic and fully turbulent state. It is worth mentioning that all the values in (3.1) are taken to be zero-mean values, therefore, u(x, t) corresponds to the turbulent fluctuations. In our study, a passive scalar with an imposed mean gradient along the x2 direction is considered to be transported with the described turbulent flow. According to the Reynolds decomposition for the total concentration of the passive scalar, Φ(x, t), one can write that Φ = ⟨Φ⟩ + ϕ. Here, ⟨·⟩ is the ensemble-averaging operator, and ϕ denotes the fluctuating part of the passive scalar concentration. More specifically, the imposed mean scalar gradient is taken to be uniform as ∇⟨Φ⟩ = (0, β, 0), where β is a constant. Therefore, the turbulent scalar concentration obeys an advection-diffusion (AD) 26 equation that is simplified into the following form ∂ϕ ∂ ∂ 2ϕ + (ϕ ui ) = −β u2 + D , i = 1, 2, 3, (3.2) ∂t ∂xi ∂xi ∂xi where D denotes the molecular diffusion coefficient of the passive scalar. Accordingly, the Schmidt number is defined as Sc = ν/D. In the LES of turbulent transport, the fluid and passive scalar motions are resolved down to a prescribed length scale namely as filter width, ∆, which linearly decomposes the veloc- ity and scalar concentration fields into the filtered (resolved) and the residual (unresolved) components. For instance, for the scalar concentration, ϕe and ϕR = ϕ − ϕe represent the filtered and residual fields, respectively. The filtered fields are obtained by a convolution, ϕe = G ∗ ϕ, where G = G(r) denotes the generic spatial filtering kernel [4]. Applying such filtering operation on the governing equations returns the subsequent LES equations. For example, the filtered AD equation is formulated as ∂ ϕe ∂ e  ∂ 2 ϕe ∂qiR + ϕuei = −β u e2 + D − , i = 1, 2, 3, (3.3) ∂t ∂xi ∂xi ∂xi ∂xi where qiR denotes the residual or SGS scalar flux that is defined exactly as qiR = ϕg ui − ϕe u ei . In the LES sense, the SGS scalar flux needs to be closed (modeled) in terms of the resolved- scale (filtered) variables through proper and physically consistent SGS modeling. 3.3 Why SGS Dynamics are Statistically Nonlocal? In an idealistic LES, one of the main elements reflecting the dynamics of turbulent trans- port is capturing the true filtered (resolved-scale) turbulent intensity through robust SGS modeling that is physically and mathematically consistent. In fact, such transport equation includes closure terms that directly link the correct time-evolution of turbulent intensity to the nature of the SGS closure and its modeling. In the LES of scalar turbulence, multiply- ing both sides of the filtered AD equation (3.3) by ϕ, e yields the time evolution of filtered turbulent intensity as 1 ∂ 2 e ∂ e  ∂ 2 ϕe ∂ qR ϕ̃ + ϕ ei = −β ϕe u ϕu e2 + D ϕe − ϕe i . (3.4) 2 ∂t ∂xi ∂xi ∂xi ∂xi 27 100 1.0 100 ∆/η = 8 10−1 ∆/η = 20 0.8 10−1 ∆/η = 41 e k) ∆/η = 53 C(qkR, G 0.6 PDF 10−2 10−2 0.4 10−3 101 102 0.2 10−4 0.0 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0 100 200 300 400 500 e −qR · G/kq R e · Gk r/η (a) (b) Figure 3.1: Statistics of true subgrid-scale contribution to the filtered scalar variance rate. (a) PDF of normalized SGS dissipation of filtered scalar variance, −q R · G,e computed over a sample space of 10 TLE of statically stationary turbulence. (b) Time-averaged two-point correlation function (3.7) between q∥R and G e with r = |r ⊥ |. ∥ Using the continuity equation and chain rule for differentiation, ! 1 ∂ 2 e ∂ ϕe ∂ ∂ e ϕ ∂ ϕe ∂ ϕe ∂  e R ∂ ϕe ei ϕ̃ + ϕ u e = −β ϕ u e2 + D ϕe −D − ϕ qi + qiR . (3.5) 2 ∂t ∂xi ∂xi ∂xi ∂xi ∂xi ∂xi ∂xi Applying the ensemble-averaging operator, ⟨·⟩, on (3.5), returns the transport equation for the filtered scalar variance, ⟨ϕ̃2 ⟩. In this study, we are considering the case of homogeneous D E ∂ turbulent velocity and scalar fields; therefore, ∂x (·) = ∂x ∂ ⟨(·)⟩ = 0. Defining the filtered i i e e scalar gradient as G(x) = ∇ϕ(x), time-evolution of the filtered scalar variance takes the following form 1d 2 ⟨ϕ̃ ⟩ = −Te + P e−χ e + ⟨Π⟩, (3.6) 2 dt D E D E D E Te = ϕe u ei , ei G e e P = −β ϕ u e2 , χ e = D Gi Gi , Π = qiR G e e ei . In (3.6), Te denotes the turbulent transport of filtered scalar variance while P e represents the e is the production of resolved scalar variance by the uniform mean scalar gradient, and χ resolved scalar variance dissipation due to the molecular diffusion. Unlike these three terms, Π (representing the SGS production of resolved scalar variance) is the only contributing term in (3.6) that contains the effects of the SGS scalar flux. Therefore, as pointed out earlier, e is essential for the SGS modeling and understanding the true statistical nature of q R · G 28 precise evaluation of the resolved scalar variance in the LES. This examination of q R · G e might be viewed both from single-point and two-point statistics as discussed in [12] in the context of the LES for homogeneous isotropic turbulent flows. In a recent comprehensive study by Di Leoni et al., effects of the SGS contribution in the evolution of the two-point velocity correlation was explored for the incompressible Navier-Stokes equations using filtered DNS data for HIT and turbulent channel flows at high-Reynolds numbers, and revealed the importance of nonlocal effects in the SGS dynamics [125]. In the present study, we are also focused on the two-point statistics of the SGS production of resolved scalar variance. This quantity is well represented in terms of the following normalized two-point correlation function D E qiR (x) Gei (x + r) C(qiR , G ei ) = D E , (3.7) ei (x) qiR (x) G where r = (r1 , r2 , r3 ) denotes the spatial shift from the location x. Moreover, probability distribution function (PDF) of the SGS production of scalar variance normalized by its L2 - e norm i.e., q R · G/∥q e is another measure to learn about the statistical behavior of Π R · G∥, and have a more comprehensive insight into the SGS modeling. High-Fidelity Database of the SGS Scalar Flux In order to study the statistics of Π, we compute true values of the SGS scalar flux using the box filtering kernel with isotropic filter width ∆ as,     1 , r ≤ ∆/2 ∆ G(r) = (3.8)   0, r ≥ ∆/2. Applying this convolution kernel on a well-resolved DNS database of passive scalar with imposed mean gradient in synthetic (forced) homogeneous isotropic turbulence. To perform the simulation, we employ an open-source parallel statistical-computational platform for turbulent transport equipped with a Fourier pseudo-spectral spatial discretization of the NS and AD equations, fourth-order Runge-Kutta (RK4) time-integration scheme, and an 29 artificial forcing method (to keep the turbulent kinetic energy at low wavenumbers constant) (see Capter 2). Our computational domain is a triply periodic cube of Ω = [0, 2π]3 that is discretized on a uniform Cartesian grid with N = 5203 Fourier collocation points while a constant ∆t = 5 × 10−4 is utilized for the stable time-integration. In construction of this DNS database, the imposed mean scalar gradient is taken as β = 1, and Sc = 1 according to the section 3.2. Letting kmax be the maximum resolved wavenumber in our simulation and η = (ν 3 /ε)1/4 be the Kolmogorov length scale while ε denotes the turbulent dissipation rate, we measure kmax η ≈ 1.5; therefore, one can ensure that the small-scales in the velocity and scalar fields are well-resolved (see Capter 2). Moreover, our records indicate that the Taylor-scale Reynolds number is Reλ = 240 (averaged over 25 large-eddy turnover times, TLE , of resolving the passive scalar field). Statistical Analysis of the SGS Effects in Filtered Scalar Intensity By taking a large sample space over 10 TLE of this stationary process (after resolving the passive scalar field for 15 TLE ), we compute the PDF of the normalized SGS production of filtered scalar variance for four different filter widths, ∆/η = 8, 20, 41, 53. As a result, we observe that as ∆ becomes larger the PDF exhibits broader tails as shown in Figure 3.1(a). Emergence of these heavy PDF tails implies that as we increase the filter width, long-range spatial interactions become stronger and more pronounced [161]. Motivated by this observation, a two-point diagnosis of the SGS scalar production of the filtered variance as defined in equation (3.7) would be another statistical measure shedding light on the long- range interactions in addition to the filter width effects. Considering ∥ as the direction along the imposed mean scalar gradient and ⊥ representing the directions perpendicular to the imposed mean gradient, we are interested in evaluating C(q∥R , G e ). Here, we take ∥ r = (r1 , 0, 0) and r = (0, 0, r3 ) and take the average of the resulting two-point correlation functions. Due to the statistically stationary turbulence, we perform such procedure for 20 data snapshots that are uniformly spaced over 10 TLE (on the same spatio-temporal data we used to compute the PDFs); hence, we obtain the time-averaged value of C(q∥R , G e ). Figure ∥ 30 1.0 Filtered DNS, ∆/η = 8 0.8 Filtered DNS, ∆/η = 53 Local EDM, ∆/η = 8 0.6 e k) C(qkR , G Local EDM, ∆/η = 53 0.4 0.2 0.0 0 100 200 300 400 500 r/η Figure 3.2: Comparison between the true values of two-point correlation function given in (3.7) and the ones obtained from the local eddy-diffusivity modeling of the SGS scalar flux given in (3.9). The evaluations are performed at two filter widths of ∆/η = 8, 53. 3.1(b) illustrates this two-point correlation function extending over a wide range of spatial shift, r = |r|, and evaluated at four filter widths similar to the ones utilized in Figure 3.1(a). This plot quantitatively and qualitatively reveals that as we increase ∆, greater correlation e (x+r) are observed values between the SGS scalar flux q∥R (x), and filtered scalar gradient G∥ at a fixed r. These spatial correlations are significant both in the dissipation and also inertial subranges. This confirms the substantial nonlocal effects in the true SGS dynamics, which needs to be carefully addressed in the SGS modeling for LES. A popular and fairly simple approach for modeling the SGS scalar flux is Eddy-Diffusivity Modeling (EDM). In EDM, the main assumption is that the SGS scalar flux is proportional to the resolved scale scalar gradient as q R (x) ≈ −DED G(x), e (3.9) and DED is the proportionality coefficient. Obviously, EDM is a local modeling approach by e ) while q R is approximated with EDM, one can com- its construction. Computing C(q∥R , G∥ ∥ pare it with its true value as shown in Figure 3.1(b). Figure 3.2 illustrates such comparison for two filter widths, ∆/η = 8, 53, reveals that in both of the cases local EDM substantially 31 fails to predict the conspicuous long-range spatial correlations observed in the true two-point correlation values. This observation is closely similar to the results reported in Di Leoni et al. [125] that showed local eddy-viscosity model is structurally incapable of reproducing the two-point SGS dissipation for the HIT and turbulent channel flows. This concrete ev- idence urges to go beyond the conventional means of SGS modeling for the scalar flux in order to address the matter of nonlocality with more sophisticated mathematical modeling tools. Thus, a nonlocal construction for the EDM would be a fairly relevant remedy to this problem. 3.4 Boltzmann Transport Framework In studying the turbulent transport and mixing, kinetic Boltzmann theory has shown a rich and promising ground based upon principles of statistical mechanics, which by con- struction is well-suited for the stochastic description of turbulence at microscopic level [172]. In the following, the fundamental sources of nonlocal closure and the SGS modeling for the residual passive scalar flux are studied at the kinetic Boltzmann transport framework. Our objective is to derive a nonlocal eddy-diffusivity SGS model at the continuum level. BGK Model and Double Distribution Function Considering classical kinetic theory of gases, we are concerned with the evolution of a single particle distribution function, f , that is governed by the Boltzmann Transport Equation (BTE), ∂f + v · ∇f = C(f ). (3.10) ∂t In (3.10), the probability distribution f = f (t, x, v) is defined such that there exists mass of fluid particles that are located inside the infinitesimal volume element dx centered at x, velocity element dv centered at v, and at time t. In the phase space of particle, x, v, and t are considered as independent variables. The left-hand side of (3.10) represents the streaming of the non-reacting particles that is balanced by the collision operator, C(f ), on the right-hand side. As a widely common model for the collision operator, Bhatnagar–Gross–Krook (BGK) 32 approximation considers scattering of the fluid particle due to collision with another particle. Therefore, the BGK model characterizes C(f ) = CBGK (f ) with a single parameter, that is called the relaxation time, τ [173]. Therefore, the collision operator is written as f − f eq CBGK (f ) = − , (3.11) τ where the local equilibrium distribution function, f eq = f eq (t, x, v) is given by the Maxwell distribution [174], and is parameterized by the locally conserved quantities (density ρ, particle speed v, and temperature T ) as ! ρ −(v − u)2 f eq = exp . (3.12) (2π c2T )d/2 2 c2T p In (3.12), cT = kB T /m is the thermal speed at T in which kB is the Boltzmann constant, and m represents the molecular air weight, while d denotes the spatial dimensions [175]. In order to study the passive scalar transport phenomena in this context, Double Dis- tribution Function (DDF) method has been a successful approach [176]. In the DDF, we consider one distribution function to address the conservation of mass and momentum while another distribution function is taken to represent the conservation of energy. In the case of passive scalar transport, the compressive work and heat dissipation are considered to be negligible in the incompressible limit [177, 178, 179]. Therefore, the extra BTE that governs the energy distribution function, g = g(t, x, v), with the BGK collision model is expressed as ∂g g − g eq + v · ∇g = CBGK (g) = − . (3.13) ∂t τg In (3.13), τg represents the relaxation time, which is the time-scale associated with the colli- sional relaxation to the local energy equilibrium denoted by the Maxwell energy distribution, ! Φ −(v − u) 2 g eq = exp . (3.14) 2 (2π cT ) d/2 2 c2T Defining L = (v − u)2 /c2T and F (L) = exp(−L/2), the Maxwell distribution in (3.14) (for the most general case where d = 3) is reformulated as g eq = Φ F (L). (2π) 3/2 c3T 33 Subsequently, continuum averaging yields the macroscopic flow variables for the incom- pressible flow, ρ = ρ(t, x), as follows: Z ρ = f (t, x, v) dv, (3.15) d ZR ρ u(t, x) = v f (t, x, v) dv, i = 1, 2, 3, (3.16) d ZR Φ(t, x) = g(t, x, v) dv, (3.17) Rd where Φ(t, x) is the total passive scalar concentration field appearing in the AD equation. Let us define L as the macroscopic characteristic length, ls as the microscopic charac- teristic length associated with the smallest length-scale of the passive scalar, and lm as the mean-free path (the average distance traveled by a particle between successive collisions). Considering x′ to be the location of particles before scattering while we characterize their current location with x, one can assume that x′ = x − δx, where δx = (t − t′ ) v. Here we assume that during the time t − t′ , v approximately remains constant [16]. According to Chen et al. [180, 181], the Boltzmann BGK kinetics with “constant” relaxation time, equa- tions (3.10) and (3.13), admit analytical solutions for f (t, x, v) and g(t, x, v) based upon their local equilibrium distribution that is valid in a general flow where the distance from the wall is large compared to lm . Focusing on equation (3.13) and defining s = (t − t′ )/τg , the exact solution to g(t, x, v) would be Z ∞ Z ∞ eq g(t, x, v) = e−s g eq (t − sτ g, x − v sτg , v) ds = e−s gs,s (L) ds, (3.18) 0 0 eq where gs,s (L) = g eq (t − sτg , x − v sτg , v). Filtered BTE, Closure Problem, and Kinetic-Boltzmann Modeling Statistical description of LES is well-represented through incorporating a filtering pro- cedure into the kinetic Boltzmann transport. For the purpose of passive scalar transport, applying a spatially and temporally invariant filtering kernel, G = G(r), onto the distribu- tion function g(t, x, v) linearly decomposes that into the filtered, ge = G ∗ g, and the residual, g ′ = g − ge, components. Therefore, filtering the equation (3.13) results in the following 34 filtered BTE (FBTE) for the passive scalar: ∂e g eq (L) ge − g^ + v · ∇ ge = − . (3.19) ∂t τg As it was elaborated by Girimaji [182], the nonlinear nature of the collision operator, CBGK (g), prohibits the filtering kernel to commute with CBGK (g); thus, it initiates a source of closure at the kinetic level in FBTE (3.19). Defining Le := (v − u e )2 /c2T , this closure problem is manifested in the following inequality, e e eq (L) = Φ exp(−L/2) ̸= Φ exp(−L/2) = g eq (L). g^ e (3.20) (2π)3/2 c3T (2π)3/2 c3T The identified closure requires proper means of modeling so that one can numerically solve the FBTE (3.19). A common practice is to approximate this closure problem with a modified relaxation time approach that is described in detail in [183]. Despite the success of this approach in some applications, it is not physically consistent with the filtered turbulent transport dynamics [182]. Nevertheless, here we manage to adjust this inconsistency by looking at the nonlocal effects arising from filtering the Maxwell distribution function, g eq (L), and model them with proper mathematical tools. Considering the spatial filtering kernel G(r) with the filter-width ∆, and applying it on the Maxwell equilibrium distribution as Z   eq (L) = G ∗ g eq L(t, v, x) = g^ G(r) g eq L(t, v, x − r) dr, (3.21) Rf where Rf = [−∆/2 , ∆/2]3 . The integral form of the convolution (3.21) implies that g^ eq (L) consists of a summation of the exponential functions. Thus, filtering encodes a multi-exponential behavior into the filtered equilibrium distribution that is gets intensified as the filter-width enlarges. Moreover, this multi-exponential structure of the filtered Maxwell distribution induces a heavy-tailed form for the filtered distribution that essentially entails the non-Gaussian behavior and jus- tifies the spatial nonlocality [16]. This statistical rationale strongly indicates that modeling this closure problem with a Gaussian-type distribution is fundamentally insufficient. On the 35 other hand, it is well-known that the statistical behavior of a multi-exponential distribution could be sufficiently approximated with a power-law distribution [184, 16]. Subsequently, by rewriting the right-hand side of the passive scalar FBTE (3.19) into the following form 1   eq (L) = − 1 g   e + 1 g^   − ge − g^ e − g eq (L) eq (L) − g eq (L) e , (3.22) τg τg τg | {z } | {z } closed unclosed the unclosed part is structurally multi-exponentially distributed and maybe approximated by a power-law distribution model as we propose e g^eq (L) − g eq (L) e ≈ g α (L)e = Φ F α (L), e (3.23) c3T where F α (L) e denotes an α-stable Lévy distribution that is mathematically designed based on heavy-tailed stochastic processes and replicate the power-law behavior [185, 105]. Regarding the decomposition given in (3.22), and by applying the filtering kernel on the analytical solution to g(t, x, v) that is given in (3.18), we obtain Z ∞ Z ∞ Z ∞   eq eq e eq eq e ge(t, x, v) = e−s g^ s,s (L) ds = e−s gs,s (L) ds + e−s ^ gs,s (L) − gs,s (L) ds, 0 0 0 (3.24) eq  eq L(t − sτ , x − v sτ , v) , and the second integral represents the closure where g^ s,s (L) = g g g source. Therefore, employing the power-law distribution model in (3.23) returns the following analytic form for ge(t, x, u) Z ∞ Z ∞ eq e ge(t, x, v) = e−s gs,s (L) ds + e−s gs,s α (L) e ds, (3.25) 0 0  wherein, gs,s e := g α L(t α (L) e − sτg , x − v sτg , v) . Fractional-Order Model for the SGS Scalar Flux Similar to the continuum averaging shown in (3.15) to (3.17), the macroscopic continuum variables associated with (3.3), are obtained in terms of the filtered distribution functions, 36 fe and ge, as Z Φe = ge(t, x, v) dv, (3.26) RZd 1 ei = u v fe(t, x, v) dv, i = 1, 2, 3. (3.27) ρ Rd i Multiplying both sides of the passive scalar FBTE by a collisional invariant X = X (v) and then integrating over the kinetic momentum would return Z   Z ! ∂e g g ge − g(L) X + v · ∇ ge dv = X − dv. (3.28) Rd ∂t Rd τg Here, choosing X = 1 would result in recovering the filtered AD equation (3.3). According to the microscopic reversibility of the particles that assumes the collisions occur elastically, the right-hand side of (3.28) equals zero [186]. Therefore, (3.28) reads as e Z ∂Φ +∇· v ge dv = 0. (3.29) ∂t Rd Since we are working with spatial filtering kernels, G = G(r), Z Z Z v ge dv = (v − u e ) ge dv + e ge dv. u (3.30) Rd Rd Rd By plugging (3.30) into (3.29), we obtain that ∂Φe   +∇· Φ eu e = −∇ · q, (3.31) ∂t where Z qi = (vi − u ei ) ge dv. (3.32) Rd Using (3.25), we formulate qi as Z Z ∞ Z Z ∞ −s eq e qi = e (ui − u ei ) gs,s (L) ds dv + e−s (vi − u α (L) ei ) gs,s e ds dv. (3.33) Rd 0 Rd 0 It is straightforward to show that the temporal shift can be removed from (3.33). Moreover, since (vi − u e and (vi − u ei ) g eq (L) e both represent odd functions of vi , thus, ei ) g eq (L) Z Z (vi − u eq e ei ) g (L) dv = (vi − u ei ) g α (L) e dv = 0. (3.34) Rd Rd 37 As a result, qi in (3.33) can be rewritten as Z Z ∞   eq e qi = e−s (vi − u ei ) gs,s (L) − g eq (L) e ds dv + (3.35) d ZR Z0∞   e−s (vi − u α (L) ei ) gs,s e − g α (L) e ds dv. Rd 0 In an LES setting, the first integral on the right-hand side of (3.35) represents the filtered scalar flux, qe, while the second integral aims to model the residual scalar flux, q R , associ- ated with unresolved small scales of turbulent transport. In other words, by assigning the Gaussian distribution g eq (L) e to qei and the isotropic α-stable Lévy distribution, g α (L), e to qiR , the total passive scalar flux, q = q e + q R , in (3.35) may be decomposed as Z ∞Z   eq e e e−s dv ds, qei = (vi − u ei ) gs,s (L) − g eq (L) (3.36) d Z0 ∞ ZR   qiR = ei ) gsα (L) (vi − u e e−s dv ds. e − g α (L) (3.37) 0 Rd In B, the details of derivation of q e and q R in terms of macroscopic transport variables including Φ e and ue are presented. As the result, the filtered passive scalar flux is obtained as e = −D ∇Φ, q e (3.38) and the divergence of residual scalar flux is derived as the fractional Laplacian of the filtered total scalar concentration, ∇ · q R = −Dα (−∆)α Φ, e α ∈ (0, 1], (3.39) Cα (cT τg )2α where Dα := τg (2α + 2) Γ(2α) is a model coefficient with the unit [L2α /T ]. The filtered AD equation for the total passive scalar concentration, developed from the filtered kinetic BTE with an α-stable Lévy distribution model, yields a fractional-order SGS scalar flux model at the continuum level. The aforementioned filtered AD equation reads as ∂Φe ∂ e  e + Dα (−∆)α Φ. e + Φuei = D ∆Φ (3.40) ∂t ∂xi 38 Through a proper choice for the fractional Laplacian order α, the developed model optimally works in an LES setting. Applying the Reynolds decomposition and considering the pas- sive scalar with imposed uniform mean gradient, equation (3.40) fully recovers the filtered e transport equation (3.3) for the transport of the filtered scalar fluctuations, ϕ. In order to explicitly derive the modeled residual scalar flux in terms of the filtered transport fields, from the Fourier definition of fractional Laplacian and the Riesz transform in given in A, one can verify that n o   n o F (−∆)α ϕe = i ξj − i ξj /|ξ| (|ξ|2 )α−1/2 F ϕe , (3.41) which leads to   (−∆)α ϕe = ∇j Rj (−∆)α−1/2 ϕe . (3.42) Therefore, using (3.39) we may write   ∇ · q R = ∇ · −Dα R(−∆)α−1/2 ϕe . (3.43) Finally, from (3.43) one can find the explicit form of the modeled SGS flux as qiR = −Dα Ri (−∆)α−1/2 ϕe + c, (3.44) where c is a real-valued constant. 3.5 Data-driven Nonlocal SGS Modeling Deriving the structure of the residual scalar flux as a nonlocal SGS model, there are two levels of model calibration in order to employ this SGS model in an LES. In fact, this model calibration problem could be viewed as a two-stage procedure where its first part is dealing with estimation of the fractional order, α, and the other stage infers the proportionality coefficient of the model, Dα . Subsequently, we propose a two-stage a priori parameter identification strategy based upon spatio-temporal data for the true q R , obtained from filtering well-resolved DNS of scalar turbulence as described in section 3.3. 39 Capturing Nonlocality with Fractional Modeling of the SGS Scalar Flux This is the first stage of this data-driven model identification, which targets finding an optimal fractional order, αopt . Our ground-truth data comes from exact evaluation of the two-point correlation function, C(q∥R , G e ) as described in section 3.3. In fact, we aim to ∥ capture the spatial nonlocality we showed in the statistics of SGS production of filtered scalar variance (see Figure 3.1b). Since we employ the fluctuating part of q∥R in comput- ing the two-point correlation function, and from the definition C(q∥R , G e ) is normalized by ∥ D E e (x) , finding αopt is essentially independent of the other model parameters ap- q∥R (x) G∥ peared in (3.44). Using the exact values of q∥R from filtered DNS, (3.7) returns the ground- truth two-point correlation function, C True Using the database described in section 3.3, while using the fractional model for SGS scalar return flux C Model as functions of spatial shift, r. In our study, for a fixed filter width, the fractional order that minimizes the mismatch function ∥C True − C Model ∥, simply determines αopt capturing the entire range of spatial nonlocality. By changing 0 < α ≤ 1, we evaluate C Model for four different filter widths, ∆/η = 8, 20, 41, 53. Figure 3.3, shows C True in addition to the variations of C Model with r/η as we change α. We observe that as α decreases, the nonlocal correlations in C True are better approximated over r with the fractional SGS model. According to the minimization of the mismatch function we introduced, αopt for the four values of filter width is reported in Table 3.1. Moreover, given αopt for each filter width, single-point correlation coefficient between   the true and modeled values of the SGS scalar flux, ϱ q∥True , q∥Model , is computed and acceptably good correlation values (in an a priori sense) are reported in Table 3.1. We need to emphasize that the passive scalar transport occurs in a statistically homogeneous medium with a direction of large-scale anisotropy. This source of anisotropy significantly impacts the intensity of nonlocal effects in the SGS dynamics so that the identified fractional-order in the SGS model is found to be less than 0.5. A similar observation in the study by Di Leoni et al. showed that presence of anisotropy effects in the turbulent channel flow (due to the non-zero mean velocity gradient along the stream-wise direction) increases the nonlocality 40 (a) ¢/¥ = 8 (b) ¢/¥ = 20 1.0 100 DNS 1.0 100 DNS Æ = 1.0 Æ = 1.0 0.8 Æ = 0.8 0.8 Æ = 0.8 10°1 Æ = 0.6 Æ = 0.6 0.6 e k) C(qkR , G Æ = 0.5 e k) C(qkR , G 0.6 10°1 Æ = 0.5 Æ = 0.45 Æ = 0.45 10°2 Æ = 0.4 Æ = 0.4 0.4 0 50 100 150 Æ = 0.35 0.4 0 50 100 150 Æ = 0.35 Æ = 0.3 Æ = 0.3 Æ = 0.25 Æ = 0.25 0.2 0.2 0.0 0.0 0 100 200 300 400 500 0 100 200 300 400 500 r/¥ r/¥ (c) ¢/¥ = 41 (d) ¢/¥ = 53 1.0 100 DNS 1.0 100 DNS Æ = 1.0 Æ = 1.0 0.8 Æ = 0.8 0.8 Æ = 0.8 Æ = 0.6 Æ = 0.6 0.6 e k) C(qkR , G Æ = 0.5 e k) C(qkR , G 0.6 Æ = 0.5 Æ = 0.45 Æ = 0.45 10°1 Æ = 0.4 Æ = 0.4 0.4 0 50 100 150 Æ = 0.35 0.4 0 50 100 150 Æ = 0.35 Æ = 0.3 Æ = 0.3 0.2 Æ = 0.25 0.2 Æ = 0.25 0.0 0.0 0 100 200 300 400 500 0 100 200 300 400 500 r/¥ r/¥ Figure 3.3: Variations of the the two-point correlation function given in (3.7) obtained for the modeled SGS flux C Model α is changed from 0 to 1, in addition to the exact evaluation e ) via exact values of the SGS scalar flux from DNS, illustrated for four filter of C(q∥R , G∥ widths (a) ∆/η = 8, (b) ∆/η = 20, (c) ∆/η = 41, and (d) ∆/η = 53. The arrows indicate the increase of α. Insets depict the two-point correlation function values on smaller regions of the spatial shift, r/η < 150, in logarithmic scale. These plots show that the true values of the two-point correlation function over the entire range of spatial shift is well-approximated with finding the αopt in the fractional-order SGS model. Table 3.1: Optimal fractional orders, and their corresponding single-point correlation coeffi- cients between true and modeled SGS scalar fluxes.   ∆/η αopt ϱ q∥True , q∥Model 8 0.40 0.35 20 0.35 0.40 41 0.36 0.44 53 0.37 0.45 in the SGS dynamics in a way that it requires α < 0.5 to properly capture that with the fractional gradient SGS model [125]. 41 Figure 3.4: Regression plots between q∥True and q∥Model for the filter widths, (a) ∆/η = 41, and (b) ∆/η = 53. The corresponding optimal fractional-orders are reported in Table 3.1. Sparse Regression on the Fractional-Order Model After obtaining the αopt for a choice of filter width, we can compute the explicit term X = R(−∆)αopt −1/2 ϕe noting the linear mapping q R = −Dα X +c, in (3.44). Having access to the true values of SGS scalar flux on an extensive spatio-temporal database (described in section 3.3) turns the second stage of our model calibration into a sparse linear regression procedure. Therefore, this procedure leads to learning and inferring of Dα that is appeared in the filtered AD equation (3.40). Similar to Beetham and Capecelatro’s work for sparse regression [187], we employ a regularized linear regression method namely as elastic net that combines the L1 and L2 penalties as its regularizer [188]. Using the implementation of the elastic net method in scikit-learn [189] and assigning equal weights to the L1 and L2 regularizes, we perform the regression and the its quality is examined through scatter plots. As a common practice, and in order to choose proper training data size, we perform cross-validation tests over our spatio-temporal dataset [75]. As a result, Figure 3.4 shows the resulting scatter plots after the regression for two cases with ∆/η = 41, 53. Using the described procedure, the proportionality coefficient for each filter width is achieved. Figure 3.5 illustrates predicted Dα through this regression procedure as a function of chosen filter width, and it is notable that the predicted Dα decreases to lower values as we chose smaller filter widths. This numerical observation is consistent with our theoretical 42 0.10 Regression ∆2 αopt 0.08 Dα 0.06 0.04 0.02 0 10 20 30 40 50 60 ∆/η Figure 3.5: Variation of the proportionality coefficient, Dα , for fractional-order SGS model with filter width, and the scale invariance study. (a) ∆/η = 41 (b) ∆/η = 53 Filtered DNS Local EDM 10−1 10−1 Fractional-order Model 10−2 PDF PDF 10−2 10−3 10−3 10−4 10−4 −0.02 −0.01 0.00 0.01 0.02 0.03 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 e −q · G/kq e · Gk e −q · G/kq e · Gk Figure 3.6: Probability distribution functions of the SGS dissipation of scalar variance, for the exact values from filtered DNS, local eddy-diffusivity model, and fractional-order SGS model at filter widths ∆/η = 41, 53. interpretation of Dα as we pointed out in section 3.4. A vital consideration in developing an SGS model is the concept of scale-invariant closure model, especially within the inertial- convective subrange [190]. As indicated in section 3.4, Dα takes the unit of [L2α /T ]. There- fore, to study the scale invariance property, by choosing the filter width as the length-scale, one can compare the variations of Dα obtained from the sparse regression against ∆2αopt . Figure 3.5 shows that the developed fractional-order SGS model is scale-invariant. 43 3.6 A Priori Testing via SGS Dissipation of the Resolved Scalar Variance We subsequently examine the capability of the optimal fractional SGS model in repro- ducing the PDF of SGS dissipation of scalar variance, q R · G. e Through addressing: • The ability of the SGS model to capture heavy tails in the true PDF, and • If the SGS model is capable of representing the backward scattering of the scalar variance cascade i.e., reproducing the negative values in the PDF. Considering two filter widths of ∆/η = 41, 53, Figure 3.6 shows the PDF of normalized SGS dissipation of filtered scalar variance for the optimal fractional-order model, local EDM, and the true SGS flux. The sample space to compute the PDFs is identical to the one we utilized to obtain the PDFs illustrated in Figure 3.1a as fully described in section 3.3. Here, one can see that for both of the filter widths the fractional-order SGS model successfully captures the broad tail of the PDF in the positive value region for the SGS dissipation, however, the local eddy-diffusivity model fails to completely do that. The positive side of the PDF is associated with the cascade of scalar variance from the resolved scales to the unresolved ones i.e., forward scattering of the scalar variance. On the other hand, this figure remarkably demonstrates that unlike the local EDM, fractional-order SGS model is able to predict the events with the negative SGS dissipation values as observed in the true SGS dissipation PDFs. In fact, our resulting PDFs display that the nonlocal modeling of the SGS scalar flux through fractional-order operator makes it possible to include the backward scattering in the LES of turbulent scalar transport. Similar observation in the context of fractional-order SGS modeling was reported by Di Leoni et al., where their developed fractional SGS model was shown to be able to reproduce the back-scattering of the filtered turbulent kinetic energy [125]. 44 CHAPTER 4 A NONLOCAL SPECTRAL TRANSFER MODEL AND SCALING LAW FOR SCALAR TURBULENCE 4.1 Background According to [132], anomalies due to the effect of large-scale motions in passive scalars occur as the result of turbulent mixing, and arising from rare events in which a parcel of fluid moves a distance much greater than the integral length scale without equilibrating. In the analogy of Lagrangian path integrals, [191] argued that this behavior is identified for a typical fluid path for which the mixing rate is anomalously long rather than for a typical mixing rate but with an atypical path [132]. This interpretation is an evidence for nonlocal interactions at the large scale levels of turbulent motion originating from the presence of anisotropy. According to [132], this behavior is directly linked to the emergence of heavy tails (exponential tails) in the PDF of passive scalar and has been experimentally observed in the turbulent behavior of passive scalars with non-zero mean gradient [192, 193, 194, 195]. A proper approach to account for a mathematical model representing the accumulative source of these nonlocal motions is to revisit the spectral transfer model for the cascade of the passive scalar. In fact, this has been a thriving area of research as reported in different studies such as [196, 197]. A nonlocal spectral transfer model provides a robust link between the large-scale anisotropy at the energy containing range and the universal range throughout the turbulent cascade while accounting for the breakdown of local isotropy at small scales. The rest of this chapter is organized as follows: in section 4.2, we introduce the math- ematical model for the transport of turbulent passive scalars, their spectral transfer view, and the nonlocal modeling for the spectral transfer. In section 4.3, we provide a detailed statistical analysis for the nonlocal and standard models in a DNS setting by comparing the single-point, two-point, and high-order small-scale statistical quantities. In section 4.4, the similarities between the current model and the fractional-order subfilter modeling for large-eddy simulation are reconciled. 45 4.2 Turbulent transport of passive scalars The Navier-Stokes (NS) equations that govern incompressible fluid flow dynamics are given by ∂u 1 + u · ∇u = − ∇p + ν ∆u + F ; ∇ · u = 0, (4.1) ∂t ρ where u is the velocity field, ρ denotes the density of fluid, p is the pressure, and ν indicates the kinematic viscosity. Moreover, F represents an external forcing mechanism and in this setting we take it as A u, where A is a linear indicator function in the spectral domain in order to artificially inject the dissipated TKE into the large scales (low wavenumbers) to generate a statistically stationary isotropic turbulent velocity field (see chapter 2). In order to model the transport of a conserved passive scalar with the diffusivity D in this medium, the advection-diffusion equation is a well-known Fickian mathematical model, which is written in the following form: ∂Φ + u · ∇Φ = D ∆Φ. (4.2) ∂t Reynolds decomposition allows for Φ = ⟨Φ⟩ + ϕ, where ⟨·⟩ denotes the ensemble-averaging operator, and ϕ is the fluctuating part of the scalar field. In the homogeneous turbulence, assumption of a uniform imposed mean gradient, ∇⟨Φ⟩ = G, is a common practice in order to consider a forcing mechanism for the turbulent intensity [154]. As a result, (4.2) is rewritten as ∂ϕ + u · ∇ϕ = −G · u + D ∆ϕ. (4.3) ∂t The governing equations are numerically solved in the pseudo-spectral setting introduced in chapter 2, and the simulation setup is further explained in section 4.2. Considering Eϕ (k, t) as the three-dimensional scalar spectrum, spectral budget of the scalar variance reads as Z ∞ ⟨ϕ2 ⟩ = Eϕ (k, t) dk, (4.4) 0 46 where k indicates the wavenumber. Through the assumption of small-scale isotropy, the spectral budget for the dissipation rate of scalar variance by molecular diffusion, χ, is given as Z ∞ χ = 2D k 2 Eϕ (k, t) dk. (4.5) 0 Depending on the ratio of ν/D, and the dissipation rate of TKE (ε), three main wavenumbers are identified for the scalar spectrum:  ε 1/4  ε 1/4  ε 1/4 kηK ≡ , kηB ≡ , kηOC ≡ , ν3 νD2 D3 associated with Kolmogorov (ηK ), Batchelor (ηB ), and Obukhov-Corrsin (ηOC ) length- scales. Unlike the case with Schmidt number Sc := ν/D ≈ 1 where all of these three wavenumbers are nearly equal, kηB and kηOC encode the different behavior of the scalar spectrum in the presence of viscous-convective subrange (Sc ≫ 1), and inertial-diffusive sub- range (Sc ≪ 1), respectively. It is convenient to differentiate the scales of turbulent cascade by the wavenumbers kEI , kDI , as the k < kEI represents the energy-containing range, and kEI < k indicates the universal equilibrium range [4]. Moreover, the universal equilibrium range is split into inertial-convective (kEI < k < kDI ), and dissipation (kDI < k) subranges for the passive scalars with Sc = 1 [196]. Scalar spectral transfer and modeling Time-evolution for the spectrum of a conserved scalar is governed by (see e.g., [4, 196]) ∂ E (k, t) − T (k, t) = P(k, t) − 2 D k 2 Eϕ (k, t), (4.6) ∂t ϕ where P(k, t) denotes the spectral content for the large-scale production rate of the scalar variance, and T (k, t) is the scalar spectral transfer function. The unknown nature of the T (k, t) causes a closure problem in (4.6); thus, a proper modeling for the spectral transfer function is required. T (k, t) could be defined as the rate of spectral flux function, F (k, t), per unit wavenumber as (see e.g., [134, 196]) ∂F (k, t) T (k, t) ≡ − . (4.7) ∂k 47 100 (a) Time-averaged spectrum 100 (b) Time-averaged spectrum C k −5/3 C k −5/3 10−1 10−1 ε−2/3 E(k) ε1/3 χ−1 Eφ(k) 10−2 10−2 10−3 10−3 10−4 10−4 10−5 10−5 10−6 100 101 102 100 101 102 k k Figure 4.1: Time-averaged 3-D spectra for (a) turbulent kinetic energy (E(k)), and (b) turbulent scalar intensity (Eϕ (k)), obtained from the DNS results described in section 4.2. By integrating (4.7) over the 3-D spectral domain, the spectral flux function is obtained. As a well-known assumption, P(k, t) mainly contributes to the energy-containing scales directly, while it is obvious that k 2 Eϕ (k, t) is mainly considerable in the small scales of the turbulent cascade where diffusion is the dominant transport mechanism [4]. Therefore, integrating (4.6) over the wavenumber space yields the following: Z k Z k ∂ 2 EI DI ⟨ϕ ⟩ = [P(k, t) + T (k, t)] dk + T (k, t) dk ∂t 0 kEI Z ∞ h i + 2 T (k, t) − 2 D k Eϕ (k, t) dk (4.8) kDI In the statistically stationary state, the second and third integrals in (4.8) are approximately zero [4, 196]. In other words, within the inertial-convective subrange ∂F (k)/∂k = 0 by the constant cascade assumption, and for the diffusive subrange ∂F (k)/∂k = −2 D k 2 Eϕ (k). As a result, for the wavenumbers in the inertial-convective subrange F (k) = χ, integrating with respect to k and employing (4.5). Subsequently, (4.8) is rewritten as: ∂ 2 ⟨ϕ ⟩ = P(t) − F (kEI , t) + F (kEI , t) − F (kDI , t) + F (kDI , t) − χ(t) (4.9) ∂t | {z } | {z } | {z } Energy-containing Inertial-Convective Dissipation | {z } Universal Equilibruim Range This picture motivates the concept of modeling for the turbulent cascade process, here specifically for the case of scalar transport. This approach was initially introduced by On- sager in [22, 23], and later was generalized by [134] to the cascade mechanism of other systems 48 such as passive scalars. In the cascade transfer, assuming the geometric doubling at each wavenumber step would approximate the step length as ∆k ≈ k. Therefore, the spectral flux function could be represented in the following form: scalar variance ∆k Eϕ (k) k Eϕ (k) F (k) ≈ = ≈ , (4.10) unit time τ (k) τ (k) where τ (k) is the time-scale associated with the step at wavenumber k [134]. Within the inertial-convective subrange, choosing τ (k) to be the local time-scale associated with the eddies of size ℓ = k −1 , reads as !−1 [k E(k)]1/2 h i−1/2 τ (k) = = k 3 E(k) . (4.11) 1/k In (4.11), E(k) is the spectrum of the turbulent kinetic energy (TKE). Thus, according to the well-known Kolmogorov’s scaling for the inertial subrange, E(k) ∼ ε2/3 k −5/3 , where ε is the dissipation rate of TKE. Plugging (4.11) into (4.10), F (k) = χ ≈ k 5/2 E(k)1/2 Eϕ (k), and according to the scaling of TKE spectrum, the scalar spectrum scales as Eϕ (k) = C χ ε−1/3 k −5/3 , (4.12) where C is the scaling constant. For a comprehensive overview of a variety of the spectral transfer models as well as the analysis of their dynamics and properties, interested readers are referred to [135, Sec. 4.7.1]. Direct numerical simulation (DNS) of the scalar turbulence with a uniform imposed mean gradient, G = (0, 1, 0), advected on statistically stationary homogeneous isotropic turbulent (HIT) flow, provides a proper database to examine the scaling law in (4.12). In this study, a well-resolved DNS at the Taylor-scale Reynolds number Reλ = 240 for the case with Sc = 1 is obtained. Resolving over sufficiently large time period in the statistically stationary state provides a rich sample space to obtain the ensemble-averaged spectra for the TKE and scalar. Figure 4.1 illustrates these time-averaged spectra over approximately 15 large-eddy turnover times, τET . Although the well-known Kolmogorov scaling for the TKE is achieved as shown 49 Time-averaged spectrum −1/2 10 −1 C k −2/3 Cα k 1.3 + k 2 ε1/3 χ−1 Eφ(k) 10−2 10−3 10−4 10−5 100 101 102 k Figure 4.2: A priori identification of the fractional order and Cα , for the modified scaling law introduced in (4.15) based upon the calibration of the scaling law with the large-scale content of the time-averaged 3-D scalar spectrum (k < 10) that induces the nonlocality. The data to compute the scalar spectrum comes from the DNS using standard scalar model. The identified values are α = 0.65, and Cα ≈ 3.9. in Figure 4.1a, Figure 4.1b reveals that the scalar spectrum does not obey the scaling law given in (4.12). Nonlocal modeling of the scalar spectral transfer In Corrsin’s generalization to the Onsager’s cascade model, ∂F/∂k is considered as the rate of gain or loss of the spectral content per unit wavenumber [134]. Moreover, this generalized model could be applied to • non-conservative systems, • systems with different characteristic time-scales, • systems with different cascading mechanisms. In fact, one can realize that the scaling given in (4.12) was obtained based upon this generalization. However, in the case of scalar spectrum, a main assumption that might be questioned is small-scale isotropy. Recently, high-fidelity computational studies showed that effects of large-scale anisotropic forcing do not vanish at the small-scales of turbulent 50 transport of passive scalars [see e.g. 32]. Moreover, these small-scale anisotropic fluctuations are identified to be highly intermittent due to intensified presence of nonlocal interactions in passive scalar turbulence. Anomalous scaling of high-order scalar structure functions is a clear and well-known experimental evidence supporting this significant deviation from local isotropy at small-scales of transport. In fact these effects are available in the inertial- convective subrange. Accordingly, we also observed that the scaling law for the scalar spec- trum has a notable discrepancy with the ensemble-averaged spectrum obtained from DNS. Based on the Corrsin’s generalization, we propose that the effects of the anisotropy could be effectively modeled in the spectral transfer function when the effect of spatial nonlocality in the cascade of scalar variance is properly considered. Inclusion of an added power-law behavior into the eddies of size ℓ mathematically enables modeling of the long-range inter- actions in the physical domain (spatial nonlocality), and would return a modification in the local time-scale given in (4.11). As a result, we propose ℓ = (k 2 + Cα k 2α )−1/2 , α ∈ (0, 1], (4.13) where Cα is a non-negative model coefficient, and Cα = 0 yields the limit case ℓ = k −1 . Consequently, the modified time-scale is derived as h i−1/2 τ (k) = (k 2 + Cα k 2α ) k E(k) . (4.14) Plugging this nonlocal time-scale into (4.10) yields the following modified scaling law Eϕ (k) = C χ ε−1/3 k −2/3 (k 2 + Cα k 2α )−1/2 , (4.15) and C denotes the scaling constant. Testing this modified scaling law on the time-averaged scalar spectrum shows a promising agreement through proper parameterization of α and Cα . In order to parameterize these two values, given the time-averaged 3-D scalar spectrum we obtained from the standard DNS (Figure 4.1b), in the modified scaling law (4.15), we initially set Cα = 1, then vary α ∈ (0, 1]. We observe that α = 0.65 yields the slope of the time-averaged scalar spectrum; however, the level of the spectral content (for the 51 scalar variance) is achieved when 1 < Cα . Therefore, by fixing the α = 0.65 and trying incremented realizations of Cα in (4.15), we realize that with Cα ≈ 3.9 we have reached to the desired parameterization. Accordingly, Figure 4.2 illustrates that with this parameterization (α = 0.65 and Cα ≈ 3.9) the universal scaling observed in the TKE spectrum is achieved for the scalar spectrum with respect to (4.15). The multi-scale nature of the turbulent cascade process implies that the nonlocal trans- port effects induced by the large-scale anisotropy are fundamentally connected to the small- scales of motion through inter-scale interactions. Here, the inertial-convective subrange essentially plays the role of a meso-scale region for the turbulent cascade, where the presence of these nonlocal inter-scale interactions are highly pronounced. In fact, several fundamental studies focused on these nonlocal interactions and tried to unravel their nature by triad and tetrad models in the spectral domain. Therefore, a modification to the dissipation model (4.5) (initiated from the small-scale isotropy hypothesis) would compliment the nonlocal ef- fects observed in larger scales of turbulent cascade. Subsequently, the total dissipation (χT ) is revised into Z ∞ χT = 2 D [k 2 + Cα k 2α ] Eϕ (k, t) dk Z0∞ Z ∞ = 2D 2 k Eϕ (k, t) dk + 2 D Cα k 2α Eϕ (k, t) dk . (4.16) 0 0 | {z } | {z } χ χα Defining Dα := D Cα , χα characterizes the essence of having an underlying nonlocal diffu- sion operator in the advection-diffusion equation governing the turbulent transport of passive scalar. From a mathematical point of view, χα directly stems from a fractional-order Lapla- cian operator, (−∆)α (·), acting on the scalar concentration; thus, the modified transport model reads as ∂ϕ + u · ∇ϕ = −G · u + D ∆ϕ + Dα (−∆)α ϕ. (4.17) ∂t 52 Table 4.1: Time-averaged values of the contributing terms in the time-evolution of scalar variance over the statistically stationary region. T P χ χα Standard -0.00056 1.1146 -1.1077 – Nonlocal -0.00033 1.1427 -0.9261 -0.2462 Numerical discretization and simulation details A standard pseudo-spectral scheme based on the Fourier-collocation method is utilized to discretize and resolve the NS and AD equations. The triply periodic computational domain Ω = [0, 2π]3 is discretized in space by a uniform grid with 5203 grid points. A fourth-order Runge-Kutta (RK4) scheme is employed to perform the time-integration with a fixed ∆t = 8 × 10−4 , while the CFL < 1 condition is always checked; therefore, the numerical stability is ensured. In this study, we select a fully developed HIT field at Reλ = 240 as the initial state of our computational experiment, and investigate the development of the passive scalar concentration under the effect of a large-scale uniform imposed mean-gradient, G = (0, 1, 0). Concentration field, ϕ(x), is initiated from zero and is resolved for approximately 30 large- eddy turnover times under the standard model (4.3), and the introduced nonlocal model (4.17), while Sc = 1. More detailed discussions about the numerical method, computational approach for the simulations, and the flow characteristics of the utilized HIT data have been presented in Chapter 2. 4.3 Statistical analysis of the nonlocal scalar turbulence model Given the fact that turbulent transport is a stochastic process, sophisticated statistical analysis of the mathematical models for such phenomena has been a center of attention in turbulence research. In this study, we consider the single- and two-point statistical quantities of interest in passive scalar transport to examine the performance of our nonlocal modeling within an equilibrium turbulent state against its conventional counterpart. 53 d dt hφ φi P T +P +χ 2 (a) T χ Scalar variance rate 1 0 −1 −2 0 5 10 15 20 25 30 t/τET 2.0 (b) d dt hφ φi P χα T χ T + P + χ + χα 1.5 Scalar variance rate 1.0 0.5 0.0 −0.5 −1.0 −1.5 −2.0 0 5 10 15 20 25 30 t/τET Figure 4.3: Records of the contributing terms in the time-evolution of scalar variance d ⟨ϕϕ⟩) defined in the right-hand side of (4.19), for: (a) standard model, (b) nonlocal ( dt fractional-order model. The time-averaged quantities (over the statistically stationary re- gion of the time-records) are reported in Table 4.1. Transport of the scalar variance Transport of the scalar variance provides an important information about the evolution of the turbulence intensity. In fact, computational fluid dynamics approach make it possible to identify and keep track of the records of different influential mechanisms obtained from the mathematical modeling of the physics. It is well known that the multiplying both sides of the AD equation would yield the time-evolution equation for the turbulent intensity, ϕ2 . Therefore, applying that to equation (4.17), after using the incompressibility condition and 54 −(P + T )/(χ + χα) 102 Standard (χα = 0) Nonlocal 101 100 0 5 10 15 20 25 30 t/τET Figure 4.4: Tracking the record of the balance in scalar variance equation ensuring the equilibrium state in simulations of standard and nonlocal models. chain rule for the spatial derivatives, one can derive that 1 ∂ϕ2 h i 2 = − ∇ · (u ϕ ) − (u ϕ) · ∇ϕ − G · (uϕ) 2 ∂t + D ∇ · (ϕ ∇ϕ) − D ∇ϕ · ∇ϕ (4.18)   + Dα ∇ · ϕ R(−∆) α−1/2 ϕ − Dα ∇ϕ · R(−∆)α−1/2 ϕ, where R(−∆)α−1/2 (·) denotes the fractional-order gradient obtained from the Riesz trans- form [16, 126]. Due to the homogeneity of the scalar fluctuations, averaging over the spatial domain is equivalent to the ensemble-averaging operation ⟨ · ⟩ [4]. Thus, applying this aver- aging operation to (4.18) and considering that homogeneity of the fluctuating fields induces ⟨∇ · (·)⟩ = ∇ · (⟨ · ⟩) = 0, the evolution of scalar variance ⟨ϕ2 ⟩ is obtained as follows: 1d 2 ⟨ϕ ⟩ = ⟨(uϕ) · ∇ϕ⟩ −⟨G · (uϕ)⟩ − D ⟨∇ϕ · ∇ϕ⟩ − Dα ⟨∇ϕ · R(−∆)α−1/2 ϕ⟩ . (4.19) 2 dt | {z }| {z } | {z } | {z } T P χ χα In (4.19), the rate of scalar variance is composed of a balance between the turbulent advection effects (T ), production by the imposed mean-gradient (P), molecular diffusion (χ), and the nonlocal diffusion (χα ). It is clear that for the standard scalar transport model in which Dα = 0, the nonlocal diffusion is consequently zero. According to the simulation considerations described in Section 4.2, we collect the records of the contributing terms in 55 the right-hand side of (4.19), and Figure 4.3 illustrates these time records for the standard (Figure 4.3a) and nonlocal (Figure 4.3b) models. Moreover, during both of the simulations we d ⟨ϕ2 ⟩, using a forward-Euler finite difference scheme, compute the rate of the scalar variance, dt and compare it with the record of the right-hand side of equation (4.19) constructed from the summation of the collected records. For both of the simulations an excellent match between these two computed quantities is observed during the entire simulation time as shown in Figure 4.3. In the current work, we are focused to examine the statistical behavior of the developed nonlocal model at the “turbulence equilibrium” state. In this context, equilibrium is interpreted when for the rate of scalar variance the following condition statistically holds: P +T ∼ 1. (4.20) χ + χα Figure 4.4 shows the displays the record of this quantity for standard and nonlocal models and we notice that after approximately two large-eddy turnover times from resolving the scalar concentration, (P + T )/(χ + χα ) starts to fluctuate around 1. In order to make sure that the transient numerical effects are well past, we continue to simulate up to t/τET = 10, and consider the rest of simulation statistics in the fully developed turbulent equilibrium state. Therefore, the time-averaging operations in our study is performed on a sample space over 10 ≤ t/τET ≤ 30. Accordingly, we can compute the time-averaged values of the contributing terms in evolution of scalar variance given in (4.19) as they are reported in Table 4.1. Comparing the time-averaged values of P from both of the models reveals that the nonlocal model approximately includes 2.5% more production rate of the scalar variance by the large-scale scalar mean-gradient. A reasonable interpretation for this observation is that once the nonlocal transfer of the scalar variance transfer in the cascade process is correctly modeled; therefore, this excessive 2.5% production rate is captured at the equilibrium state for scalar turbulence. In other words, devising a nonlocal turbulent dissipation model (χα ) in the scalar variance cascade mechanism would enable a balance in the equilibrium state so that the nonlocal effects in turbulent transport originating from large-scale “anisotropy” source are better captured throughout the DNS. 56 Time-averaged spectrum −1/2 ε1/3 (χ + χα)−1 Eφ(k) 10−1 C k −2/3 DDα k 2α + k 2 10−2 10−3 10−4 10−5 100 101 102 k Figure 4.5: Time-averaged scalar spectrum computed from the data simulated with the nonlocal model, and evaluation of the identified scaling law in (4.15) for the scalar variance spectrum. Finally, we compute the time-averaged scalar variance spectrum obtained from the scalar field resolved by the nonlocal fractional-order model to examine the modified scaling law (4.15) a posteriori. Figure 4.5 depicts this spectrum and reveals that the scaling (4.15) seamlessly holds. It is worth emphasizing that the total turbulent dissipation is denoted by χ + χα . High-order small-scale statistics of scalar fluctuations It is well-understood that statistics of the turbulence at the small-scales of the transport are represented through the central moments of the gradients of the fluctuating fields. Here, we are interested in discovering the small-scale statistics when the scalar field is sufficiently resolved with the proposed nonlocal scalar transport model. In fact, we seek to understand what would be the prediction of this model for the asymmetric and highly intermittent nature of passive scalar turbulence at the small scales. Therefore, we compute the skewness and flatness factors for the fluctuating concentration gradient, and due to the importance of the statistical behavior along the anisotropy direction ∥, we focus on S ∇ ϕ and K ∇ ϕ . Figure ∥ ∥ 4.6 illustrates the records of these two statistical quantities throughout the entire simulation for the standard and nonlocal models. Over the equilibrium state, explained in section 4.3, 57 2.5 (a) Standard Nonlocal 2.0 S ∇k φ 1.5 1.0 0.5 0.0 0 5 10 15 20 25 30 t/τET (b) 30 K ∇k φ 20 10 Standard Nonlocal 0 0 5 10 15 20 25 30 t/τET Figure 4.6: Time records of (a) skewness, and (b) flatness of the scalar gradient along the anisotropy direction labeled by ∥. The time-averaged values are identified with dashed lines over the statistically stationary state, and their values are reported in Table 4.2. we obtain the mean values of the of these statistical quantities by time-averaging, and their values are reported in Table 4.2. These time-averaged values show that the nonlocal model yields the skewness factor 10% more than the standard model, and the flatness factor is approximately 7% higher in the nonlocal transport model. On the other hand, in a study by [198] on the resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence, for a similar problem setup at Reλ = 240 and Sc = 1, they performed two sets of standard direct numerical simulations with resolutions (a) N = 5123 (kmax ηB ≈ 1.5) and (b) N = 20483 (kmax ηB ≈ 5.1). The resulting high-order statistics we obtained and reported in Table 4.2, are in great agreement with what [198] reported for the same case in Table 4 of their work (µ3 and µ4 ). In fact, our study shows that our nonlocal modeling and the performed DNS based on that, predicts the high-order statistics with about 3% 58 Table 4.2: Time-averaged values of S ∇ ϕ , and K ∇ ϕ over the statistically stationary state ∥ ∥ as illustrated in Figure 4.6. S∇ ϕ K∇ ϕ ∥ ∥ Standard model 1.40 20.8 Nonlocal model 1.54 22.2 error compared to the extra high-resolution simulation (the case with kmax ηB ≈ 5.1) in Donzis & Yeung’s study reports. This comparison essentially proves the effectiveness of our modeling while we reduce the computational cost of the simulation about 64 times [compared to simulation results with N = 20483 spatial resolution in 198]. This evidently implies that an appropriate modeling of the nonlocal turbulent scalar transfer via the fractional-order model properly reflects the well-known statistical features of highly non-Gaussian behavior of the passive scalar turbulence reported in the literature [132, 34]. Two-point statistics and structure functions Structure functions of order n for a turbulent field such as scalar concentration are defined as: ⟨(δr ϕ)n ⟩ = ⟨[ϕ(x + r) − ϕ(x)]n ⟩, n > 1. (4.21) In (4.21), r = re where r is the increment of spatial shift, and e denotes a unit vector along a direction of interest. In fact, the structure functions would provide the nth-order statistics of spatial increments in the fluctuating field, which are interesting metrics in studying the nonlocality. In this study, we are interested in analyzing the behavior of the second- and third order structure function of ϕ along the direction of large-scale anisotropy, i.e. e = (0, 1, 0), and regrading the size of the DNS grid, r = 2 η. Accordingly, Figure 4.7 shows the time- averaged (over the equilibrium turbulent region) structure functions of order 2 and 3 obtained from the simulations from standard and nonlocal scalar transport models. In Figure 4.7(a), one can observe that for r > 40 η the nonlocal second-order structure function starts to exhibit higher values compared to the one computed from the simulations using standard 59 (a) (b) 100 h(δr φ)2i h(δr φ)3i 0 10 10−1 10−1 Standard Nonlocal 101 102 101 102 r/η r/η Figure 4.7: Time-averaged nth-order scalar structure functions obtained from the simulations with standard and nonlocal model, with r = 2 η. (a) n = 2, and (b) n = 3. 101 Standarad ∼ r3 ∼ r2 10 0 Nonlocal ∼ r3 ∼ r −hδr uL (δr φ)2i 10−1 10−2 10−3 10−4 10−5 100 101 102 103 r/η Figure 4.8: Third-order mixed longitudinal structure function, representing the statistics of advective increments. The nonlocal model shows a consistent and extended scaling over universal range. model. For the third-order structure function values, the two models behave similarly up to r/η = 10; however, after that the nonlocal model shows higher values withing the spatial shift domain associated with the inertial-convective and integral-scale domain. It is apparent that the maximum value of the time-averaged ⟨(δr ϕ)3 ⟩ in the nonlocal model is approximately 10 times higher than the standard model both occurring at r/η ≈ 200. As initially introduced in [25], mixed “velocity-scalar” third-order structure function is 60 an importasnt two-point statistical quantity measuring the advective turbulent transport in passive scalars. In the presence of large-scale anisotropy (imposed mean scalar gradient), the longitudinal contribution to this mixed structure function plays the dominant role in the advective transport [199], and its functional form obtained as −⟨δr uL (δr ϕ)2 ⟩. Here, the subscript L indicates the velocity component along the longitudinal direction with respect to spatial shift direction r, where in our computational setup it would be u2 . Similar to the the second- and third-order scalar structure functions, we compute a time-averaged value for the −⟨δr uL (δr ϕ)2 ⟩ over the stationary time domain. Figure 4.8 shows that for the dissipative range (r/η < 6) this structure function scales with r3 in both standard and nonlocal models, while for almost the entire range of 6 < r/η < 200 the mixed structure function obtained from the DNS with the nonlocal transport model scales with r2 . Unlike this universal-range scaling, one can observe that similar behavior is not necessarily seen in the −⟨δr uL (δr ϕ)2 ⟩ when the scalar field is resolved with the standard model. However, in the standard model, a scaling with r could be identified within 20 < r/η < 60. This comparison suggests that the considering the nonlocal effects in the turbulent cascade could result in emergence of more universal behavior in the two-point statistics of the advective transport, which inherently reveal high-order statistics of the nonlocality. 4.4 Reconciliation with the fractional-order SGS modeling for LES Large-eddy simulation is known to be an effective technique in computational turbulence research that reduces the computational cost of the simulations by focusing on resolving the larger scales of the transport while the unresolved scales are modeled from the resolved-scale transport quantities [88]. From a theoretical point of view, the turbulence closure appearing in the LES equations is the result of applying a general filtering operation to the governing equations. In the convolution kernel G∆ (x′ − x) for this filtering operation, ∆ is considered to be the arbitrary filter size. In LES, the common practice is to take ∆ large enough towards the intermediate scales of turbulent transport. However, in theory, the filter size could be considered close to the smaller scales of transport in a way that ∆ → η [190]. 61 10−2 0.08 10−3 0.06 10−4 Dα 1 2 3 4 0.04 0.02 Evaluation from filtered DNS Prediction with GPR DNS-level grid (filter) 0.00 95% confidence interval 100 101 ∆/η Figure 4.9: Reconciliation of the nonlocal model with the fractional-order SGS model devel- oped in [126] when the filter size is assumed to be at the dissipation range of ∆ = 2 η. The Dα is computed from the filtered DNS data for ∆/η = 4, 8, 20, 41, 52 based on the data-driven methodology introduced in [126] and a Gaussian process regressor (GPR) is trained based upon these evaluations, and Dα is approximated based on the trained GPR for ∆/η < 4 The predicted Dα for ∆ = 2 η is found to be in total agreement with the identified one obtained from the scaling analysis a priori in Figure 4.2. With this rationale, the simulations we performed in this study might be seen as an LES when filter size is twice the η. Subsequently, it is interesting and vital to examine if the developed nonlocal transport model is reconciled with a nonlocal functional SGS model in terms of a priori model identification. The fractional-order model for the SGS scalar flux proposed in [126] seems to be the appropriate candidate for this examination. Recalling from section 4.2, we performed an a priori parameter identification that yielded α = 0.65 and Dα /D ≈ 3.9 (see Figure 4.2). Now, in order to fulfill our goal, we need to show that given that the optimal fractional order for the SGS model αopt = 0.65, what would be the value of Dα /D that is obtained from the procedure introduced in [126] that relies on explicitly filtered data and sparse regression. Taking the filtered data from the simulation based on the standard scalar transport model (4.3) with the time-averaged scalar spectrum shown in Figure 4.2, we can obtain the proportionality coefficient for the fractional-order 62 SGS model. Here, we choose a top-hat box filtering kernel and obtain the filtered data for the filter sizes ∆/η = 4, 8, 20, 41, 52. Our goal is to evaluate Dα when ∆/η = 2; however, it is not computationally possible to obtain the filtered data for infer the Dα at . Instead, we manage to employ a feasible machine learning algorithm (ML) to predict the desired Dα while it is trained on the evaluated Dα values from direct filtered data at larger filter sizes. Gaussian Process Regression (GPR) is a known to be a suitable ML algorithm when one seeks to predict a quantity of interest from scarce experimental or computationally expensive high-fidelity data. Using the implementation of GPR in Scikit-Learn package [189], we obtain the predicted value of Dα /D = 3.87 for ∆ = 2 η as illustrated in Figure 4.9. This result shows that the a priori estimates for the proportionality coefficient in the nonlocal scalar transport model is in great agreement with the fractional SGS model when filter size is selected as ∆ = 2 η; therefore, both models are reconciled. It is worth mentioning that the uncertainty in the predictions of the trained GPR for ∆/η < 4 is assessed, and it is observed that the uncertainty is very low and practically negligible (see the 95% confidence interval in Figure 4.9). 63 CHAPTER 5 DEEP LEARNING MODELING FOR SGS FLUX IN THE LES OF SCALAR TURBULENCE AND GENERALIZATION TO OTHER TRANSPORT REGIMES 5.1 Background In the SGS modeling for SGS scalar flux, Vollant et al. [200] managed to represent the closure flux with a deep learning model, and tested their model for homogeneous isotropic turbulent flows. Later, Portwood et al. [157] proposed a deep neural network model for the SGS scalar flux and trained their model on an extensive dataset obtained from filtering the direct numerical simulation of a passive scalar with uniform mean gradient in a HIT flow. They evaluated their model in a priori (off-line) and a posteriori (on-line) tests and showed their model significantly outperforms the functional and structural gradient-type SGS models. In another study, Frezat et al. [103] developed a deep learning SGS scalar flux model using a convolutional neural network assembly, while they enforced multiple transfor- mation invariance properties (such as translation, rotation, and Galilean invariance). They showed that enforcing such physical properties is crucial for improving the performance of SGS model. Looking back to Vollant et al.’s work, they mentioned that were not success- ful in testing their trained SGS model on a different scalar transport regime [200]. In 2-D turbulence settings, Subel et al. [201] and Guan et al. [202] conducted successful generaliza- tion of pre-trained deep learning SGS models to predict on higher Reynolds number regimes using transfer learning technique [203]. Inspired by these pioneering studies, here in the context of 3-D homogeneous turbulence with a uniform scalar gradient, we seek to develop a deep learning SGS scalar flux model that could feasibly be generalized to other turbulent transport regimes. The remaining parts of this chapter is organized as: In section 5.2, we present the method- ologies for numerical simulations and data-driven modeling for base and transfer learned DNN models. In section 5.4, we go over a prioi and a posteriori tests of the trained models 64 and compare them to the conventional gradient-based SGS models for scalar flux. 5.2 Modeling Paradigm Governing Equations Considering the incompressible flow regime and the transport of a conserved passive scalar in that flow medium, Navier-Stokes (NS) and Advection-Diffusion (AD) equations are the set of governing equations constituting the system dynamics [4]. In the LES of turbulent transport, a generic spatial filtering operator, ˜·, is applied to the NS and AD equations returning the LES governing equations (see e.g., [88]) as ∂ ũ 1 + ũ · ∇ũ = − ∇p̃ + ν ∆ũ − ∇ · τ R ; ∇ · ũ = 0, (5.1) ∂t ρ ∂ Φ̃ + ũ · ∇Φ̃ = D ∆Φ̃ − ∇ · q R . (5.2) ∂t In this set of equations, u = (u1 , u2 , u3 ), p, and Φ are the velocity, pressure, and scalar concentration fields, respectively. In (5.1), ρ denotes the fluid density, while ν represents the viscosity of fluid, and in (5.2), D is the diffusivity of the passive scalar field. Moreover, filtering yields sources of closure in the LES governing equations as the divergence of residual stress, τ R = u^ ⊗ u− ũ⊗ ũ, and residual scalar flux, q R = u gΦ− ũ Φ̃. Modeling these residual or subgrid terms using the filtered or resolved flow fields is an essential gateway returning a closed set of equations that are suitable for a predictive and numerically stable LES [190, 88]. The Reynolds decomposition for a general field such as scalar concentration, Φ = ⟨Φ⟩+ϕ, where ⟨Φ⟩ is the ensemble-averaged part of Φ, and ϕ denotes its fluctuating part [4]. In our problem setting, we consider a homogeneous isotropic medium for velocity field; therefore, ⟨u⟩ = 0. For the passive scalar field, we assume the fluctuations are statistically homogeneous while we consider an imposed ensemble-averaged gradient as ∇⟨Φ⟩ = (0, 1, 0) [204, 8]. As a result, the filtered AD equation (5.2) is rewritten as: ∂ ϕ̃ + ũ · ∇ϕ̃ = −ũ2 + D ∆ϕ̃ − ∇ · q R . (5.3) ∂t In (5.3), the residual scalar flux may be restated as: q R = ufϕ − ũ ϕ̃. The goal of our study 65 Figure 5.1: Schematic of the DNN surrogate model for prediction of SGS scalar flux. σ in the units (neurons) indicates the activation function. is to developing a data-driven model for q R using a predictive deep neural network (DNN) architecture. Deep Neural Network (DNN) Closure Modeling Using a deep feed-forward neural network architecture, we model the residual scalar flux. This model has L fully connected layers each with N l neurons in the l-th layer, a weight l l−1 l tensor of W l ∈ RN ×N , and a bias vector bl ∈ RN . For instance, in the l-th layer, ynl formulates the output of the n-th neuron as  l−1  MX ynl = σ l  Wnml y l−1 + bl  , m n (5.4) m where σ l indicates the activation function. In this DNN model, y 0 represents the input features’ layer while y L = q DN N is the layer corresponding to prediction of the targets. Through optimization of the DNN’s parameters (W , b) with FDNS data, q DN N predicts q R = (q1R , q2R , q3R ). The input features are selected based on the gradients of the filtered (resolved) flow fields, i.e. ∇ũ and ∇ϕ̃. Choosing these input features preserves important properties such as frame invariancy and homogeneity of the SGS flux in predictions of the DNN model [157]. 66 5.3 Numerical and Optimization Methods Simulation setup for training/validation data Based on the DNS of passive scalar with a uniform mean gradient of ∇⟨Φ⟩ = (0, 1, 0) in a HIT flow, our required SGS data is constructed. Here, we employ the pseudo-spectral parallel simulation setup elaborated in Chapter 2. Using the mentioned framework we generate a stationary HIT flow using a computational domain as Ω = [0, 2π]3 , which is discretized on a uniform spatial grid of N 3 resolution with N = 520. Accordingly, the size of spatial grid is ∆DNS = 2π/N . A fourth-order Runge-Kutta (RK4) scheme is utilized to perform the time-integration with a fixed ∆t = 8 × 10−4 , while the CFL < 1 condition is checked during the simulation; therefore, the numerical stability is always ensured. Aiming for the Taylor- scale Péclet number (P eλ = Sc Reλ = 240) with the Schmidt number (Sc = ν/D = 1) and = Reλ representing the Taylor-scale Reynolds number, at fully turbulent (statistically stationary) state, we obtain well-resolved velocity and scalar turbulent fields with kmax ηB ≈ 1.5; therefore, we ensure that the small-scales in the velocity and scalar fields are well- resolved (see Chapter 2). In order to reach to this turbulent state, first the NS equations are resolved from a randomly initialized field for approximately 15 τLE while an artificial forcing mechanism is enforced to the low wavenumbers (energy containing range) to maintain the turbulent kinetic energy (see Chapter 2). Afterwards, we start resolving the AD equation from ϕ0 (x) = 0 initial fluctuating concentration while imposing a uniform mean-gradient. By resolving the NS and AD equations for approximately another 15 large-eddy turnover times, the skewness and flatness records of the fluctuating passive scalar gradient reach to a statistically stationary state, specifying the fully-developed turbulent scalar fluctuations. Over a sampling time span in the fully turbulent region, the true values of the residual scalar flux are computed using the box filtering kernel with isotropic filter size ∆ as,    1/∆, |x − x′ | ≤ ∆/2 ′ G∆ (x − x ) = (5.5)   0, Otherwise. 67 (a) Full-size spatial data (input and target features) obtained from FDNS at a sampled time in- stance. For ∇ũ tensor only the component 2,2 , and for ∇ϕ̃ and q R vectors only the second components are shown on the 3-D spatial domain. 50 40 30 x3 20 10 0 50 Sampling domain 40 30 0 10 20 x2 20 30 10 x1 40 50 0 (b) Temporal sampling domain, shown on statistical (c) Spatial index sampling to be applied on records of ∇ϕ. each time-sampled full-size 3-D data. Figure 5.2: Spatio-temporal sampling framework for construction of training/validation dataset from FDNS results for scalar transport regime of P eλ = 240, Sc = 1. Construction of a Spatio-temporal Database In our study, we examine our model in two filter sizes ∆∗ = ∆/∆DNS = 10, 20, where the amount of filtered scalar variance for these two filter sizes are 9% and 18%, respectively. Given the time-stationary statistics of our turbulent flow data set, uniform sampling over fully turbulent time domain is a fair approach. Therefore, we select a sampling window over approximately 11 τLE of fully turbulent state, and uniformly sample 64 time instances. To prepare a randomly generated data set in spatial domain (3-D periodic cube) for each time 68 Figure 5.3: Training/validation of the DNN surrogate model. Training MSE value after each gradient decent iteration for (a) ∆∗ = 10, and (c) ∆∗ = 20. Averaged training MSE over evaluations in each epoch for (b) ∆∗ = 10, and (d) ∆∗ = 20. All of the depicted validation MSEs in (a)-(d), shown with dashed red lines, are averaged MSE values over the validation data batches. sample, we utilize Latin Hypercube Sampling (LHS) to obtain a randomly drawn index set based on a uniform distribution. This index set returns the location of the spatially sampled data points on the original 3-D uniform grid of the filtered DNS solution. For instance, for the case with ∆∗ = 10, a set of 233 spatially scattered data points is sampled from the in the full spatial data set of 523 resolution. This reduction of the data at a sampled time frame helps to construct an unbiased database in space and time (using sampled time snapshots of FDNS data) for training/validation of the DNN model, thus, we avoid an unnecessary large size of the entire database. As a result, this spatio-temporal sampling procedure yields a reference FDNS database of approximately 780,000 samples. We hold 5% of the data out for out-of-sample validation during the optimization procedure. Figure 5.2 illustrates this sampling framework. 69 DNN Specifications and Optimization In this study, the DNN model consists of 6 nonlinear hidden layers with rectified linear unit (ReLU) activation functions as: σ(y) = max(y, 0), (5.6) and each hidden layer contains 32 neurons (units). Unlike the nonlinear hidden layers, the output layer that returns the DNN approximation for q R , is set to be linear. As a result, this DNN model yields 6,000 trainable parameters and these model parameters are initialized based on a uniform distribution prior to model optimization. In order to train the DNN model, the database from spatio-temporal sampling is utilized to minimize the mean-squared error (MSE) loss function N 1 X MSE = (qˆi − q i )2 . (5.7) N i In (5.7), q i denotes the DNN prediction for each sample and qˆi is the reference labeled data from FDNS. Using the AdamW optimizer with a learning rate of lr = 10−4 , we manage to minimizing the MSE loss function with respect to the model parameters (weights and biases) of the DNN as we utilize PyTorch library [205]. Using a batch size of 256 data points while we shuffle the training data in the data loader, we perform model optimization for 150 epochs (with each epoch consisting of 2800 optimizer iterations). At the end of each epoch we perform out-of-sample validations with the updated model. The same batch size is selected for the validation data and the resulting MSE evaluations are averaged over the validation batches. For the DNN model optimized based on data set of ∆∗ = 10, Figure 5.3a shows the MSE evaluations of the training iterations, and Figure 5.3b depicts their averaged value over the iterations in each epoch. In both of these plots, the averaged MSE evaluations over the validation data batches are shown with dashed line. At the end, the convergence of the MSE is achieved, and it is shown that the over-fitting is avoided as the averaged MSE values for training and validation are the same over the entire optimization procedure. We 70 Figure 5.4: Schematic of the transfer learning from a DNN model pre-trained on data for Task A to a DNN model for Task B. In the initial training for Task A (upper model in the plot), all of the model parameters are optimized; while in the training the pre-trained model for Task B, the model parameters in the earlier layers are fixed, and the ones in the final layers are optimized based on the data for Task B (lower model in the plot). observe the same behavior in model training/validation using ∆∗ = 20 filtered data set as shown in Figure 5.3c and Figure 5.3d. Transfer Learning for Generalization to Other Transport Regimes Due to the increasing computational cost of high-resolution simulations as P eλ number increases in scalar transport, it is not always feasible to obtain a sufficiently large dataset to properly train a DNN to reliably model the SGS flux. On the other hand, it is widely recognized that the data-driven surrogate models for physical systems do not generalize well to accurately predict in physical regimes outside of their training data. Transfer Learning (TL) provides a practical approach to this problem, and lets us adjust our pre-trained sur- rogate model (from a physical regime with abundant training data) to work well on another physical regime with less available data. Transfer learning aims to re-train the very final 71 layers of the pre-trained DNN to adjust the model parameters learning specific features that belong to the new physical regime utilizing its available data. The rationale behind this re-training step is that the general features that are common among the initial and new physical regimes are learned in the earlier hidden layers (in the pre-trained model), while the specific features are learned in the later hidden layers. Therefore, one can freeze the model parameters associated with the earlier layers in the re-training procedure (by setting them as non-trainable parameters in the optimization), and only update those associated with the later hidden layers using the data for the new physical regime. Figure 5.4 provides a schematic view of the transfer learning in our DNN modeling for SGS scalar flux. In the previous section, we presented a trained DNN model using FDNS dataset from a turbulent scalar transport at regime P eλ = 240 and Sc = 1 (as described in section 5.3. Taking that pre-trained model, we manage to re-train it for transfer learning to the following scalar transport regimes for the specified filter sizes: • P eλ = 240 and Sc = 4 (∆∗ = 10, 20), • P eλ = 360 and Sc = 1 (∆∗ = 8, 16). The data for re-training these two transport regimes were obtained from filtering DNS results where the small scales of passive scalar are sufficiently resolved (kmax ηB ≈ 1.5). The details of these simulations are similar to the description provided in section 5.3. In the re-training step for each transport regime, we construct training/validation datasets that contain approximately 5% of sampled data points we utilized for the pre-trained model (P eλ = 240 and Sc = 1). The construction of these training/validation datasets are based on the procedure we explained in section 5.3; however, in these datasets, we only use 4 sampled time instances of full-sized 3-D FDNS data for spatial sampling. The batch size, and learning rate are chosen similar to optimization of the pre-trained model, yet in the re-training steps we only aim for 50 epochs. As a result, the number of gradient decent iterations in the re-training steps are approximately 2.5% of iterations in the pre-training. 72 Figure 5.5: Training/validation for the transfer learning from the pre-trained DNN model DNN using data of P eλ = 240, Sc = 1 regime (see Figure 5.3) to a DNN model for P eλ = 360, Sc = 1 regime. Training MSE value after each gradient decent iteration for (a) ∆∗ = 8, and (c) ∆∗ = 16. Averaged training MSE over evaluations in each epoch for (b) ∆∗ = 8, and (d) ∆∗ = 16. All of the depicted validation MSEs in (a)-(d), shown with dashed red lines, are averaged MSE values over the validation data batches. Moreover, in the re-training steps, we optimize the model parameters in the final two layers of the pre-trained DNN model, and the rest of the model parameters in other layers remain fixed. For example, Figure 5.5 shows the records of the loss function for the re-training and validation for the transfer learning to the SGS model at transport regime P eλ = 360, Sc = 1. Figures 5.5a-b show the convergence of the MSE for the model re-trained based on data with ∆∗ = 8 (pre-trained model utilized data with ∆∗ = 10), while Figures 5.5c-d illustrate the converged model when ∆∗ = 16 in the re-training data (pre-trained model utilized data with ∆∗ = 20). Similar to the pre-trained models presented in Figure 5.3, close behavior of the training and validation MSE records during the re-training steps ensures that the transfer learned models are not overfitted. 73 5.4 Model Assessments in the Inference Mode In this section, we manage to examine the performance of the trained DNN models described in the previous section in the inference of the SGS and resolved flow quantities. Here we employ the full-size 3-D time instances of flow fields (ũ, ϕ̃) coming from either the FDNS solution (a priori test) or from an LES that utilizes the SGS model (a posteriori test). In both of these tests, we assess the performance of the DNN model against the two traditional closure models for the SGS scalar flux: • Prandtle-Smagorinsky (PSM) model, • Scalar Asymptotic Gradient (SAG) model. In PSM model, the SGS scalar flux is determined based on the eddy-diffusivity concept, and is written as q PSM = −Dt ∇ϕ̃, (5.8) where Dt = νt /P rt denotes the turbulent diffusivity. The constant turbulent Prandtl num- ber, P rt , relates Dt to the turbulent viscosity, νt = (Cs ∆)2 ∥S̃∥ defined by the Smagorin- sky model. Here Cs is the Smagorinsky constant and for the filtered strain-rate tensor   p 1 T S̃ = 2 ∇ũ + ∇ũ , ∥S̃∥ = 2 S̃ : S̃. For passive scalars in homogeneous turbulence, if the filter size is selected withing the inertial subrange, Cs ≈ 0.17 and P rt ≈ 0.5 [206, 207]. Therefore, we adopt these values in our study. In the SAG model, the residual scalar flux is approximated based on a Taylor series expansion of the filtering operation [208, 209], and is represented as ∆2 q SAG = ∇ϕ̃ · ∇ũ. (5.9) 12 Although this model has shown a good structural performance in reproducing the SGS flux from FDNS flow fields, it has been identified to be significantly under dissipative especially when the filter size is sufficiently large compared to the dissipation scales [10]. Therefore, 74 conducting a stable LES (in terms of numerical time-integration) using SAG model is known to require special modifications such as utilizing a “clipping” technique [210, 211] or com- bining it with an eddy-diffusivity model [208]. A priori Test In the a priori testing of developed DNN models, we employ full-size 3-D FDNS solutions for ϕ̃ and ũ. We utilize 10 time instances sampled from approximately 5 τLE away the time domain we obtained our training/validation datasets. Using these time instances of FDNS data and the SGS model, we predict the SGS flux on the 3-D filtered domain. Considering the mean scalar gradient ∇⟨Φ̃⟩ as the anisotropic source of turbulence production for passive scalar, we categorize the predicted SGS flux terms into perpendicular (q⊥∇⟨Φ̃⟩ = q1R , q3R ) and parallel (q∥∇⟨Φ̃⟩ = q2R ) components with respect to the mean scalar gradient direction. Moreover, we compute the SGS dissipation of filtered scalar variance, Π = −q R · ∇ϕ̃, (5.10) as an indicator for capability of the SGS flux model in reproducing the backward scattering of the scalar variance. Backward scattering phenomenon is the transfer of turbulent intensity from unresolved (SGS) scales to the resolved (filtered) scale, and it is identified with the negative values appearing in the true Π (obtained from FDNS). Once we compute q R and Π for each time instance using the SGS model, we take the mean squared error of each quantity with respect to its FDNS solution. In order to average the computed MSEs over the time instances we sampled, each MSE is normalized by the variance of its FDNS solution. Averaged normalized MSEs for q⊥∇⟨Φ̃⟩ , q∥∇⟨Φ̃⟩ , and Π are presented for each model. For q⊥∇⟨Φ̃⟩ , the averaged normalized MSE is the mean value of MSEs for both perpendicular components. Finally, the probability distribution function (PDF) of reconstructed Π from the SGS model is compared to the true PDF from FDNS. In the following, we present these metrics for the base SGS model we developed in section 5.3 as well as the two transfer learned models to other transport regimes we presented in section 5.3. 75 Figure 5.6: A priori estimation (base model - P eλ = 240, Sc = 1) of the time-averaged normalized mean squared error of (a) q⊥∇⟨Φ̃⟩ , (b) q∥∇⟨Φ̃⟩ , (c) Π for each SGS flux model at each filter size. The normalization is done by the variance of FDNS solution for each quantity at its filter size, and the time-averaging is carried out using 10 sampled time instances. Figure 5.7: A priori assessment of the PDFs of normalized SGS dissipation obtained from the DNN (base model - P eλ = 240, Sc = 1), SAG, and PSM models compared to the PDF computed from FDNS for (a) ∆∗ = 10, (b) ∆∗ = 20. Π for all the PDFs is normalized by the standard deviation of its FDNS solution. Base model – P eλ = 240, Sc = 1 Here we compare the performance of the developed DNN model in section 5.3 to the performance of SAG and PSM models. Our comparison is based on the time-averaged normalized MSE metric for predicted q⊥∇⟨Φ̃⟩ , q∥∇⟨Φ̃⟩ , and Π at ∆∗ = 10, 20 (see Figure 5.6). Figure 5.6a shows that at ∆∗ = 10, the DNN model prediction of q⊥∇⟨Φ̃⟩ has 5.5% less error compared to the SAG model and 76.4% less error compared to the PSM model, while at ∆∗ = 20, its prediction returns 11.6% less error compared to the SAG model and 76 54.3% less than the error in the PSM model’s prediction. In Figure 5.6b, we see that at ∆∗ = 10, the DNN model predicts q∥∇⟨Φ̃⟩ with 9.8% less error compared to the SAG model and with 88.2% less error compared to the PSM model, and when ∆∗ = 20, the DNN model yields 18.9% less error compared to the SAG model and its error is 81.1% less than what we observe in the prediction of PSM model. Moreover, from Figure 5.6a,b one can realize that the error level in prediction of q⊥∇⟨Φ̃⟩ , q∥∇⟨Φ̃⟩ components, the DNN model exhibits a fairly consistent behavior; however, such consistency is not necessarily observed in the error level of the predicted flux components in other two models, especially at ∆∗ = 20. Finally, for the predicted Π, Figure 5.6c illustrates that the DNN model has 25.7% less error compared to the SAG model and it returns 82.9% less error compared to the PSM model. However, at ∆∗ = 20, the predicted Π by the DNN model returns 60% lower error compared to the SAG model and 31.4% less than what the PSM model predicts. Here, it is noticeable that when the filter size increases the error level in predicted Π form SAG model increase while the PSM model exhibits an opposite behavior. Nevertheless, the DNN model exhibits a consistent behavior over the filter sizes, while it always maintains error level in the predicted Π significantly lower than the other two models. Given this observation, Figure 5.7 shows that for both filter sizes the DNN model provides a fairly good prediction for the PDF of SGS dissipation of scalar variance compared to the true PDF from FDNS, while the SAG model under-predicts the events in the right side of the PDF (yielding and under-dissipated SGS model), and the PSM model provides an over-dissipated model since it cannot model the negative-valued events associated with the backward scattering of resolved scalar variance. Transfer learned models In this section, we make a comparison between the time-averaged normalized MSEs of q⊥∇⟨Φ̃⟩ , q∥∇⟨Φ̃⟩ , and Π obtained from a priori testing of the TL models we presented in section 5.3 with the error metrics when we only employ the base (pre-trained) model for the predictions in the new scalar transport regime. Case (I) P eλ = 240, Sc = 4: In Figures 5.8a and 5.8b, we can see that at ∆∗ = 10, the 77 Figure 5.8: A priori estimation of transfer learned model (P eλ = 240, Sc = 4) from the base pre-trained model. The subplots are showing the time-averaged normalized mean squared error of (a) q⊥∇⟨Φ̃⟩ , (b) q∥∇⟨Φ̃⟩ , (c) Π for pre-trained and TL SGS flux model at each filter size. The normalization is done by the variance of FDNS solution for each quantity at its filter size, and the time-averaging is carried out using 10 sampled time instances. Figure 5.9: A priori estimation of transfer learned model (P eλ = 360, Sc = 1) from the base pre-trained model. The subplots are showing the time-averaged normalized mean squared error of (a) q⊥∇⟨Φ̃⟩ , (b) q∥∇⟨Φ̃⟩ , (c) Π for pre-traine and TL SGS flux model at each filter size. The normalization is done by the variance of FDNS solution for each quantity at its filter size, and the time-averaging is carried out using 10 sampled time instances. TL model decreases the error in predictions of q⊥∇⟨Φ̃⟩ and q∥∇⟨Φ̃⟩ by 12% and 46.5%, respectively, compared the predictions of pre-trained model. Similar comparison at ∆∗ = 20 indicates that TL model returns 27.3% less error for predicted q⊥∇⟨Φ̃⟩ , and 93.2% less error for the predicted q∥∇⟨Φ̃⟩ compared to the pre-trained DNN model. Therefore, one can realize that the re-training procedure plays a significant role in improving the prediction of parallel component of the SGS flux. Moreover, Figure 5.8c shows that the predicted Π from pre-trained DNN model has 2.6 times higher error at ∆∗ = 10 and 3.3 times higher error at ∆∗ = 20 compared to the TL model’s predictions. In fact, Figures 5.10a and 5.10b illustrate that these improvements in prediction of Π mainly originate from correcting the 78 Figure 5.10: A priori assessment of the PDFs of normalized SGS dissipation obtained from the transfer learned model (for P eλ = 240, Sc = 4), and pre-trained model (using P eλ = 240, Sc = 1 data) compared to the PDF computed from FDNS for (a) ∆∗ = 10, (b) ∆∗ = 20. Π for all the PDFs is normalized by the standard deviation of its FDNS solution. Figure 5.11: A priori assessment of the PDFs of normalized SGS dissipation obtained from the transfer learned model (for P eλ = 360, Sc = 1), and pre-trained model (using P eλ = 240, Sc = 1 data) compared to the PDF computed from FDNS for (a) ∆∗ = 8, (b) ∆∗ = 16. Π for all the PDFs is normalized by the standard deviation of its FDNS solution. 79 over-prediction of forward scattering of SGS dissipation after re-training the DNN model for higher Sc number regime. Case (II) P eλ = 360, Sc = 1: Considering Figures 5.9a and 5.9b, we observe that at both ∆∗ = 8, 16 the TL model reduces the error level in prediction of q⊥∇⟨Φ̃⟩ by 22 times compared to the pre-trained model. However, in prediction of q∥∇⟨Φ̃⟩ , the error level in pre-trained model’s performance is reduced by 33 times at ∆∗ = 8, and by 37 times when ∆∗ = 16 after re-training. Similar to the Case (I), we recognize that TL is playing a more corrective role in enhancing the performance of the pre-trained DNN model in prediction of q∥∇⟨Φ̃⟩ compared to q⊥∇⟨Φ̃⟩ . Concerning the prediction of Π, Figure 5.9c points out that the pre-trained DNN model has 49.5 times higher error at ∆∗ = 8 and 60.5 times higher error at ∆∗ = 16 compared to the TL model. Looking at Figures 5.11a and 5.11b, it is discernible that in the case of TL from P eλ = 240 regime to P eλ = 360, the re-trained DNN model not only corrects the considerable over-prediction of forward scattering phenomenon appearing in the pre-trained model (see the positive sides of PDFs), it also improves its over-prediction of backward scattering of the SGS dissipation (see the negative sides of PDFs). A Posteriori Test In this section, our main objective is to examine the performance of the developed DNN model when it is utilized as the closure model in an LES setting, a.k.a a posteriori testing. Therefore, we seek to understand if the LES resolved transport fields (in our case ϕ̃) are temporally and spatially following their expected true solution which comes from explicit filtering of the DNS outputs. To fulfill this goal, we focus on testing the base DNN model we developed in section 5.3 on the turbulent regime of P eλ = 240 and Sc = 1. In order to perform the large-eddy simulations on the problem setting introduced in section 5.2, we employ the pseudo-spectral solver developed in Chapter 2. We modify this DNS framework to account for utilization of multiple SGS flux models of interest for q R . In our current study, we utilize the DNN model as well as the PSM and SAG models. Typically, both SGS stresses and fluxes in 80 the filtered NS and AD equations are modeled and utilized in an LES setting. Thus, a specific choice of SGS model for τ R could have substantial effects on the solution resolved scalar concentration field [157]. Similar to earlier works in literature for pure purpose of SGS flux model assessment [157, 212], we choose to freeze these potentially dominant effects originating from using an SGS stress model and instead directly make use of FDNS solution of velocity field. Therefore, we can only focus on the performance of the utilized SGS flux model in (5.3). Accordingly, in our numerical setup, we resolve the NS equations based on the DNS required resolution. By explicit filtering of u after each time-step, we obtain the velocity field over the desired LES grid, and use it to solve the equation (5.3) in the time-integrator steps. The initial condition (ϕ̃0 and ũ0 ) for our LES tests are adopted from explicit filtering of a well-resolved and fully developed DNS solution, at P eλ = 240, Sc = 1 as described in section 5.3. As we mentioned earlier, for the SAG model we require to employ special treatment for nu- merical stabilization in LES. Here, we utilize clipping technique during the time integration, which is based on setting the predicted q SAG = 0 before feeding it to the time-integrator where ever ΠSAG < 0. This essentially helps to improve the extremely under-dissipated behavior of the SAG model only for numerical stability [210, 211]. In our tests, the time- integration of filtered AD equation (5.3) is conducted for 5 τLE . In the following, we look at two significant indicators in performance of the SGS scalar flux model. Temporal records of resolved scalar variance Precise time evolution of the resolved scalar variance, ⟨ϕ̃2 ⟩, is a principal indicator for reliable prediction of the turbulent intensity in LES. In our tests, we manage to record this quantity for the utilized SGS models over the simulation time, while the reference record is obtained from filtering the DNS solution. Here, we compute the relative error between the variance record from the LES with an SGS model ⟨ϕ̃2 ⟩LES and from the FDNS ⟨ϕ̃2 ⟩FDNS . Figure 5.12a shows the recorded errors for LES tests at ∆∗ = 10, and 5.12b reports it for simulations at ∆∗ = 20. At both of the filter sizes, the DNN model has a remarkable 81 Figure 5.12: Relative error (with respect to the FDNS) in the records of resolved-scale scalar variance predicted in the LES tests with the DNN, SAG, and PSM models at (a): ∆∗ = 10, and (b): ∆∗ = 20. The LES tests are conducted at P eλ = 240 and Sc = 1 regime, and the utilized DNN model is trained/validated with FDNS dataset obtaned at the same turbulent regime. performance over time compared to the SAG and PSM models, however, in ∆∗ = 20, slight under-prediction of resolved-scale varinace is observed. It is obvious that in both of the filter sizes, SAG model exhibits quite a poor performance especially as time evolves; and this behavior is worse as filter size enlarges. On the other hand, the dissipative nature of the PSM model seems to be helpful in terms of a better performance of the model in LES; nevertheless, its considerable errors in prediction of Π (as we observed in section 5.4) seems to penalize its accuracy in the a posteriori tests. Resolved-scale scalar structure functions (two-point statistics) Assessment of an SGS model performance in prediction of the nonlocal behavior of tur- bulent transport regime is one of the ultimate tests in evolution of a turbulent field during LES [12]. The structure functions of the resolved scalar field are robust two-point statistical measures that return nth-order statistics of resolved-scale scalar increments at a specific direction where 2 ≤ n [132]. These structure functions of order n are defined as D n E ⟨δr ϕ̃n L⟩ = ϕ̃L (x + reL ) − ϕ̃L (x) ; n = 2, 3, . . . , (5.11) where r is the size of spatial increment, L represents the longitudinal direction (the direction along the imposed uniform mean-gradient) [199, 157, 19], and eL specifies the unit vector 82 Figure 5.13: Time-averaged normalized ⟨δr ϕ̃n L ⟩ from LES with different SGS models with respect to the FDNS result as reference solution. For ∆∗ = 10, (a) second-order and (b) third-order illustrate longitudinal structure functions. For ∆∗ = 20, (c) second-order and (d) third-order show the longitudinal structure functions. The time-averaging of these two-point structure functions is done over 4 < t/τLE < 5. The spatial shift r equals the filter size of the LES. along the longitudinal direction. Normalizing the ⟨δr ϕ̃n L ⟩ computed from LES (carried out with an SGS model) with its FDNS measurement would provide an informative metric for comparing the SGS models. Focusing on the temporal region 4 < t/τLE ≤ 5 that the LES solution has undergone a fairly long time-integration, we select uniformly distributed samples of full-size 3-D ϕ̃(x) in time and compute second- and third-order structure func- tions. Since we are simulating a statistically-stationary problem, we are allowed to take the temporal average of these normalized structure function obtained from the sampled struc- ture functions. Therefore, we have a robust indicator measure to examine the performance of each SGS model in predicting nonlocal and high-order statistics of resolved-scale scalar 83 field in a long-time integrated LES. Figure 5.13 shows this time-averaged ⟨δr ϕ̃n L ⟩ against the normalized spatial shift, r/ηB , for the DNN, SAG, and the PMS models. Figures 5.13a and 5.13c are showing the relative errors is the even-order structure functions, where the DNPS model has a considerably better performance. Figures 5.13a and 5.13c are showing the normalized second-order structure functions, for ∆∗ = 10 and ∆∗ = 20, respectively. Remarkably, we observe that for both of the filter sizes predictions of the DNN model closely follow the temporally-averaged behavior of ⟨δr ϕ̃2L ⟩FDNS over the entire range of spatial shift. On the other hand, we realize that the SAG model has a significant over prediction of ⟨δr ϕ̃2L ⟩ especially for r < 100ηB , which becomes worse for the larger filter size. Even in longer range spatial shifts that we observe improvement in the its predictions, there is an asymptotic 20% over prediction at ∆∗ = 10, and approximately 100% over prediction of ⟨δr ϕ̃2L ⟩ at ∆∗ = 20. Conversely, the PSM model seems to exhibit a better performance in prediction of ⟨δr ϕ̃2L ⟩ compared to the SAG model; However, it has an approximate 20% over prediction of the second-order structure function when ∆∗ = 10, and at ∆∗ = 20, a prolonged 50% over- prediction of ⟨δr ϕ̃n L ⟩ is seen (especially for 100 ηB < r). Moreover, it is well understood that capturing the complex behavior of the third-order scalar structure function is quite a cumbersome task [132]. However, by looking at Figures 5.13b and 5.13d, we notice that the DNN model manifests a great performance in prediction of the third-order longitudinal structure function over almost the entire range of spatial shift at both of the investigated filter sizes. For example, at ∆∗ = 10, we observe that over the whole range of spatial shift the DNN model predictions has an approximate ±5% deviation from the FDNS values; nonetheless, we can observe an over-prediction causing inaccuracies over r in the range of 10% to 40% for the SAG model. Regarding the PSM model, the inaccurate predictions of ⟨δr ϕ̃3L ⟩ varies from -15% deviations from the true value for r < 100 ηB to approximately 20% deviation for 100 ηB < r spatial shift region. Therefore, the DNN model plays a notable role in maintaining the accuracy of the third-order statistics of δr ϕ̃L . This overall comparison demonstrates that a well-trained DNN model such as the one we 84 developed could have consequential improvements in the prediction of two-point statistics of the LES resolved scalar concentration. 85 CHAPTER 6 ANOMALOUS TRANSPORT IN INTERNAL CYLINDER FLOW INSTABILITIES SUBJECT TO UNCERTAINTY 6.1 Background Inspired by the flow dynamics after the emergence of symmetry-breaking factors, we are specifically interested in computational study of the onset of flow instabilities and their long-time effects. To model such symmetry-breaking effects in rotational motion of cylinder, we introduce some featured sources of disturbance in angular velocity, which may be cou- pled by eccentricity in rotation of the system. In reality, these sorts of symmetry-breaking noises could be a direct result of unexpected failure in the electro-mechanical rotational sys- tem/fixture, which may be accompanied by secondary inertial disturbances that intensify the instability and transition of the flow regime. From a mathematical modeling and sim- ulation point of view, a deterministic view would inevitably fail to reflect the true physics of such highly complex phenomenon, which is involved with numerous sources of stochas- ticity, (i.e., sources of disturbance). This urges for another level of modeling and investi- gation, which respects the random nature of the problem and is capable of addressing the effects of such sources of randomness in the response of system. In general, these sources of randomness could be categorized into either aleatory or epistemic model uncertainties. Aleatory uncertainty affects the quantities of interest (QoI) by the natural variations of the model inputs and usually are hard to be reduced; nevertheless, epistemic uncertainty mostly comes from our limited knowledge on what we are modeling and could be stochas- tically modeled once we obtain additional information about the system [213]. Uncertainty in modeling procedure and also inaccuracy of the measured data are two main factors in arising epistemic uncertainty. The uncertainty in modeling could be the result of a vari- ety of possibilities including the effects of geometry [214, 215, 216, 217], constitutive laws [218, 219, 220, 221, 222, 223, 224, 225, 226, 227], rheological models [228, 229, 230], low- fidelity and reduced-order modeling [231, 232, 233, 234, 235, 236, 237], and random forcing 86 sources in addition to the random field boundary/initial conditions [238, 239, 115]. The structure of the rest of this chapter is outlined as follows: In section 6.3, we formulate the stochastic version of the Navier-Stokes equations for incompressible flows and develop our stochastic modeling procedure. In section 6.4, we elaborate on the numerical methods we employ in our deterministic solver and generation of a proper grid and later on we introduce the our stochastic discretization approach followed by a discussion on how we study the significance of each source of stochasticty in a global sense. In section 6.6, we show the stochastic convergence, quantification of uncertainty in kinetic energy as QoI and we perform the global sensitivity analysis. Using the expected velocity and vorticity fields we computed from our stochastic computational framework, we obtain the fluctuating responses for a deterministic simulation and study their statistics in a qualitative and quantitative sense. Furthermore, we compute the enstrophy record associated with the fluctuating field and study its time-scaling that unravels a tied link between the observed highly non-Gaussian features and memory effects induced by long-lived coherent vortex structures. 6.2 Problem Statement The present study aims to investigate the fluid dynamics inside a rotating cylinder with radius R that is assumed to be fully-filled with the Newtonian fluid with the kinematic viscosity, ν. Here, the fluid is initially considered to be at solid-body rotation state with an angular velocity of θ̇0 = dθ/dt|t=0 that is enforced by the rotational speed of the cylinder wall. Considering the configuration of the problem, the rotational Reynolds number is defined as Re = R2 θ̇0 /ν. As a common practice and for the sake of generality in comparisons, we manage to formulate our problem in a dimensionless format so that R and θ̇0 are taken to be unity. Consequently, the Reynolds number is simply computed as inverse of the kinematic viscosity of the fluid, Re = ν −1 . The solid-body rotation state represents a laminar flow regime that we take as the initial stage of the flow where right away it encounters a mixture of symmetry-breaking disturbances in the rotational motion of the cylinder. These sources of rotational disturbance include an angular velocity for cylinder with oscillatory and decaying 87 amplitude that is assumed to be accompanied by an eccentric rotation of the cylinder (i.e., resulting from the rotation of cylinder about an off-centered axis). The combination of these factors would make a strong symmetry-breaking effect on the flow that we study the dynamics of the resulting flow instability. In such study, the rotating cylinder is assumed to be long and we consider a 2-D representation of the flow as an acceptably good approximation. In the following section, we go through a detailed mathematical model for the described sources of disturbance considering their randomness that requires a stochastic modeling procedure for the fluid dynamics study. Therefore, we proceed with presenting the stochastic governing equations. 6.3 Stochastic Navier-Stokes Equations Let Ω ⊂ R2 be our bounded convex 2-D spatial domain with boundaries ∂Ω. Moreover, let (Ωs , F, P) be a complete probability space, where Ωs is the space of events, F ⊂ 2Ωs denotes the σ-algebra of sets in Ωs , and P is the probability measure. Then, the govern- ing stochastic incompressible 2-D Navier-Stokes (NS) equations subject to the continuity equation, ∇ · V = 0, for Newtonian viscous fluids ∂V + V · ∇V = −∇p + ν ∇2 V , ∀(x, t; ω) ∈ Ω × (0, T ] × Ωs , (6.1) ∂t V (x, t; ω) = V ∂Ω , ∀(x, t; ω) ∈ ∂Ω × (0, T ] × Ωs , V (x, 0; ω) = V 0 , ∀(x; ω) ∈ Ω × Ωs , hold P-almost surely subject to the corresponding proper initial and boundary conditions, introduced and modeled below. Here, V (x, t; ω) represents vector of the velocity field for the fluid, p(x, t; ω) denotes the specific pressure (including the density). Stochastic Modeling We are interested in learning how the symmetry-breaking factors would affect the onset of flow instability. In our modeling, these factors are reflected in terms of stochastic initial and boundary conditions, subsequently, the rest of possible random effects are treated deter- ministically. Accordingly, these symmetry-breaking effects are modeled through imposing a 88 time-dependent wall angular velocity, θ̇(t; ω) = cos (α(ω)t) e−λ(ω)t , ∀(t; ω) ∈ (0, T ] × Ωs , (6.2) while we consider an off-centered rotation with a radial eccentricity of ϵ(ω), ∀ω ∈ Ωs , with respect to the geometric centroid of the cylinder. In our model, α(ω) and λ(ω) denote the frequency of oscillations and the decay rate appearing in the angular velocity model, respec- tively. In other words, no-slip boundary condition at the wall is imposed by the proposed wall velocity for which the initial condition is a solid-body and off-centered rotation. Re- calling that in our non-dimensional mathematical setup, the initial angular velocity, θ̇(0; ω), and the radius of the cylinder, R, are both taken to be unity; therefore, the stochastic wall velocity field is expressed as   V ∂Ω (x, t; ω) = x − rϵ(ω) θ̇(t; ω), ∀(x, t; ω) ∈ ∂Ω × (0, T ] × Ωs , (6.3) ∥x∥2 = 1, ∥rϵ(ω) ∥2 = ϵ(ω), where ∥ · ∥2 denotes the L2 norm. Moreover, according to (6.2) θ̇(t = 0; ω) = 1, so that the initial condition is denoted as V 0 (x; ω) = x−rϵ(ω) , ∀(x; ω) ∈ Ω×Ωs where ∥rϵ(ω) ∥2 = ϵ(ω). Parameterization of Random Space Let Y : Ωs → R3 be the set of independent random parameters, given as Y (ω) = {Yi }3i=1 = {λ(ω), α(ω), ϵ(ω)}, ∀ω ∈ Ωs , (6.4) with probability density functions (PDF) of each random parameter being ρi : Ψi → R, i = 1, 2, 3, where Ψi ≡ Yi (Ωs ) represent their images that are bounded intervals in R. Q By independence, the joint PDF, ρ(ξ) = 3i=1 ρi (Yi ), ∀ξ ∈ Ψ, with the support Ψ = Q3 3 i=1 Ψi ⊂ R form a mapping of the random sample space Ωs onto the target space Ψ. Thus, an arbitrary point in the parametric space is denoted by ξ = {ξ 1 , ξ 2 , ξ 3 } ∈ Ψ. According to the Doob-Dynkin lemma [240], we are allowed to represent the velocity field V (x, t; ω) as V (x, t; ξ), therefore, instead of working with the abstract sample space, we 89 rather work in the target space. Finally, the formulation of stochastic governing equations in (6.1) subject to the boundary/initial conditions in equation (6.3) can be posed as: Find V (x, t; ξ) : Ω × (0, T ] × Ψ → R such that ∂V + V · ∇V = −∇p + ν∇2 V , (6.5) ∂t V (x, t; ξ) = V ∂Ω , ∀(x, t; ξ) ∈ ∂Ω × (0, T ] × Ψ, V (x, 0; ξ) = V 0 , ∀(x; ξ) ∈ Ω × Ψ, hold ρ-almost surely for ξ(ω) ∈ Ψ and ∀(x, t) ∈ Ω × (0, T ] subject to the incompressibility condition, ∇ · V = 0. 6.4 Stochastic Computational Fluid Dynamics Framework Discretization of Physical Domain and Time-Integration Spectral/hp element method [241] is a high-order numerical method to discretize the governing equations (6.1) in the deterministic physical domain Ω. In particular, SEM is a proper candidate to achieve a high-order accuracy discretization close to the wall boundaries. SNel e In SEM, we partition the spatial domain, Ω, into non-overlapping elements as Ω = e=1 Ω , where Nel denotes the total number of elements in Ω. In practice, a standard element, Ωst , is constructed in a way that its local coordinate, ζ ∈ Ωst , is mapped to the global coordinate for any elemental domain, x ∈ Ωe . This mapping is performed through an iso-parametric transformation, x = χe (ζ). Within the standard element, a polynomial expansion of order P is employed to represent the approximate solution, V δ , as N P Ndof X el X X V δ (x) = V̂je Φej (ζ) = V̂i Φi (x), (6.6) e=1 j=1 i=1 where Ndof indicates the total degrees of freedom (DoF) i.e., the modal coefficients in the solution expansion. Moreover, Φej (ζ) are the local expansion modes, while Φi (x) are the global modes that are obtained from the global assembly procedure of the local modes [241]. NEKTAR++ [242, 243], a parallel open-source numerical framework, provides a seamless platform offering efficient implementation of multiple SEM-based solvers in addition to the 90 10−1 Error 10−2 10−3 0 200000 400000 600000 800000 1000000 1200000 DoF (a) (b) Figure 6.1: (a) Constructed structured grid with transitional h-refinement. (b) Grid conver- gence study based on the error in kinetic energy. pre-/post-processing tools. In our study, we employ its incompressible Navier-Stokes solver namely as IncNavierStokesSolver. Here, the velocity correction scheme along with the C 0 -continuous Galerkin projection are utilized as splitting/projection method in order to decouple the velocity and the pressure fields [242]. We use P th-order polynomial expansions i.e., the modified Legendre basis functions while we vary P for elements at different spatial regions (see section 6.4). Moreover, a second-order implicit-explicit (IMEX) time-integration scheme is used while the time-step is fixed during the time-stepping. The spectral vanishing viscosity (SVV) technique [244, 241] is also used to ensure a stabilized numerical solution from spectral/hp element method. Grid Generation A 2-D structured grid is generated with quadrilateral elements considering h-type refine- ment technique to attain proper grid resolution near the wall. We utilize the open-source finite element grid generator, Gmsh [245], to construct the geometry and then the h-refined grid. The generated grid is illustrated in Figure 6.1a, which shows elemental nodes and h-refinement near the wall. For this h-refined grid, we employ a spatially-variable polyno- mial expansion [243] so that we gain high-accuracy close to wall, while avoiding unnecessary computational cost away from the wall. In order to ensure that our solution is independent of the grid resolution for the Reynolds number that is fixed at Re = 106 , we carry out a grid 91 convergence study based on the error we obtain from the difference of the time-integrated kinetic energy between the solutions after varying the grid resolution and a reference solution with ∼ 2.1 × 106 total DoF. As shown in Figure 6.1b, the total DoF of ∼ 7.5 × 105 gives us a sufficient grid resolution ensuring that the numerical solution is independent of grid resolu- tion. In the applied IMEX time-integration scheme, the time-step is fixed at ∆t = 4 × 10−5 while the numerical stability is always checked during the simulations by ensuring that CFL number being less than unity. In particular, our SEM grid is achieved by utilizing 9th-order polynomial expansions for the elements in the near the wall region and 7th-order polynomial expansions for the elements in the cylinder’s core region. In other words, due to this spatial p-refinement procedure, the near-wall elements would consist of 64 rectangular sub-elements (P = 9) and, the elements in the core region will be finer 36 times (P = 7). For flow at moderately low Reynolds numbers, we verify the resulting solutions from our numerical setup through a comparison with analytical solutions (see Appendix C). Stochastic Discretization Sampling from the parametric random space introduced in section 6.3 is a non-intrusive approach for stochastic discretization. Monte Carlo (MC) sampling method is the most conventional way to perform such task, however, the large number of required realizations of random space is its bottleneck, which prohibits utilizing MC for computationally demanding problems. In our study, we employ probabilistic collocation method (PCM) [246, 247, 248], which is a non-intrusive scheme and has shown affordable efficiency by providing fairly fast convergence rate for statistical moments. In PCM, a set of collocation points {q j }J j=1 is prescribed in parametric random space Ψ, where J denotes the number of collocation points. As a common practice to construct a stable basis, {q j }J j=1 are taken to be the points of a suitable cubature rule on Ψ with integration weights, {wj }J j=1 . In this work, we employ a fast algorithm proposed by Glaser et al. [249] to compute the collocation points based on Gauss quadrature rule. Therefore, let the solution V in the parametric random space be collocated on the set of points {q j }Jj=1 . In other words, we use the SEM setup 92 described in section 6.4 to solve a set of deterministic problems in which the wall velocity field V ∂Ω (x, t; ξ) in equation (6.5) is replaced with its deterministic realization V ∂Ω (x, t; q j ). In order to construct the approximate stochastic solution V̂ (x, t; ξ) from a set of deterministic solutions {V (x, t; q j )}J j=1 , we employ Li (ξ), the Lagrange interpolation polynomials of order i. Let I represent the approximation operator, therefore, the approximate stochastic solution is written as X J V̂ (x, t; ξ) = IV (x, t; ξ) = V (x, t; q j )Lj (ξ). (6.7) j=1 We choose the approximation operator I to be the full tensor product of the Lagrange interpolants in each dimension of parametric random space. Defining the PDF ρ(ξ) over the parametric random space and using the approximate solution, the expectation of V is computed as Z E [V (x, t; ξ)] = V̂ (x, t; ξ)ρ(ξ)dξ. (6.8) Ψ This integral would be approximated using a proper quadrature rule. Letting the set of interpolation/collocation points {q j }J j=1 obtained from Glaser et al. [249] coincide these quadrature points with associated integration weights {wj }J j=1 , one can efficiently compute the approximation to the integral in equation (6.8). Applying the Kronecker delta property of Lagrange interpolants, this integral is approximated as XJ E [V (x, t; ξ)] ≈ wj ρ(q j ) J V (x, t; q j ). (6.9) j=1 In equation (6.9), J represents the Jacobian associated with an affine mapping from standard to the real integration domain regarding the applied quadrature rule. In our study, we utilize uniformly distributed random variables to represent symmetry-breaking effects, hence, J ρ(q j ) yields a constant. In the case of our problem with three stochastic dimensions, J ρ(q j ) = ( 21 )3 . Consequently, the approximate computation of the expectation integral 93 (6.8) is simplified to J 1X E[V (x, t; ξ)] ≈ wj V (x, t; q j ). (6.10) 8 j=1 Similar to the MC approach and using (6.10), the standard deviation in our problem is approximated as v u u X J  2 u1 σ [V (x, t; ξ)] ≈ t wj V (x, t; q j ) − E[V (x, t; ξ)] . (6.11) 8 j=1 6.5 Variance-Based Sensitivity Analysis Grasping knowledge on the significance of sources of randomness in a stochastic model- ing procedure could be very helpful in terms of reducing the computational cost and also decision making during stochastic modeling. Variance-based sensitivity analysis is a well- known technique to assess the relative effect of randomness in each stochastic dimension on the total variance of any QoI, U , as the output of a stochastic model in a global sense [250, 251]. In practice, sensitivity of the QoI to each stochastic parameter is measured by the conditional variance in the QoI, which is caused by that specific parameter. In general, for a k-dimensional stochastic space, ξ, a QoI may be represented as a square-integrable function of the stochastic parameters U = f (ξ). Using Hoeffding decomposition of f [252], and also   the conditional expectation of the stochastic model, E U |ξ i (i=1,...,k) , the total variance of U can be decomposed as X XX V (U ) = Vi + Vij + · · · + V12...k , (6.12) i i i 0.8 at all recorded times, while the effects of the other parameters are always less than 0.2. In particular, by focusing on t < 0.75, we realize that oscillatory effect of the angular velocity model embodied in ξ 2 , is the second dominant source of randomness propagated in the kinetic energy of the entire system, nevertheless, after t = 0.75 as the dynamics of instability evolves with time, the effect of oscillations in the angular velocity decreases. In fact, when 0.75 < t the eccentric rotation is the only effective mechanism appearing in the uncertainty of kinetic energy. On the other hand, by following the summation of the first-order sensitivity indices P depicted in Figure 6.6, we observe that i S i > 0.95, which reveals that the joint interactions of the stochastic parameters on the total variance of kinetic energy are negligible. However, presence of these joint interactions is slightly realized close to the onset of the instability when t < 0.75. Statistical Analysis of Fluctuating Flow Fields Emergence of fluctuating flow velocity field plays a key role in the dynamics of flow instabilities. For instance, Ostilla et al. [257] studied the behavior time-averaged root- 100 (a) (b) Figure 6.7: Snapshots of velocity fluctuations at ti = 0.025, 0.2, 0.375, 0.75 for i = 1, . . . , 4. (a) Radial velocity fluctuations, u′r (x, t), (b) azimuth velocity fluctuations, u′θ (x, t). mean-square (r.m.s.) of the velocity fluctuations to study the dynamics of boundary layer in different regimes of Taylor-Couette flow. In another study, Grossmann et al. [258] examined the behavior of velocity fluctuations profile in a strong turbulent regime of Taylor-Couette problem. In this regard, here we seek to shed light on the mechanism of initiating the flow instability from a statistical perspective through studying the behavior of the fluctuations. In principle, any instantaneous field variable such as velocity, V , which contains a fluctuating part could be decomposed into V = V + V ′, (6.18) where V ′ represents the fluctuations of V and V denotes its ensemble average. Unlike the applied approach in [257, 258] that approximates the ensemble average by time-averaging over a time period on developed flow, here we are not allowed to exploit time-averaging close to the onset of the instability, which essentially takes place in a short period of time. However, our stochastic modeling and CFD platform enables us to properly approximate the ensemble-averaged velocity field with reasonable computational cost. Hence, having the knowledge of ensemble mean velocity field gives us the fluctuating response of the flow field variables. The fluctuations are appeared in the flow at the existence of stochasticity and disturbance in the system. In fact, the ensemble mean is nothing but finding the mathematical expectation of the field variable over the entire sample space that contains 101 large enough number of realizations. Thus, what we obtain as the result of equation (6.10) is the representation of ensemble mean in a PCM setting [259]. The stochastic convergence analysis we performed in section 6.6 ensures that the expectation we compute from PCM is independent of the stochastic discretization, therefore, we are allowed to claim that the expected velocity field on the sufficiently converged PCM is a robust approximation of its ensemble average with large enough number of independent samples. As a result, we can write V = E [V (x, t; ξ)] . (6.19) According to the sensitivity analysis we performed in section 6.6, we are allowed to obtain the ensemble-averaged field by performing a uni-variate PCM on the most sensitive stochastic parameter, ξ 3 = ϵ, while we fix the other two random parameters of the wall velocity model to their mean values as reported in the Table 6.1. Since the uni-variate PCM requires much less realizations evaluated at collocation points, it is computationally feasible to discretize the dominant random direction even beyond the stochastic convergence resolution. Here we proceed with taking 30 collocation/integration points providing a high-resolution expected solution in the stochastic space essentially returning a seamless evaluation of V . According to Table 6.1 and as a physically reasonable assumption, the rotational eccentricity is initially taken to be varying up to 5% of the cylinder radius as ϵ ∼ U (0.0, 0.05). For a randomly drawn realization of the sample space that fixes eccentricity value at ϵ = 0.0263, we evaluate the fluctuating velocity field according to equation (6.18). The procedure of computing the fluctuations from SEM-based realizations is briefly explained in Appendix D. Figure 6.7 shows the resulting velocity fluctuations in polar coordinate system at four snapshots of time illustrating the onset of flow instability. Emergence of Non-Gaussian Statistics in Velocity Fluctuations Tracking the probability density function (PDF) of velocity fluctuations with time ren- ders qualitative statistical information, which characterizes the impacts of the evolution of fluctuations on the dynamics. PDF of fluctuating fields can simply show us the departure 102 t = 0.125 10−1 10−1 t = 0.15 t = 0.175 10−2 10−2 t = 0.2 Gaussian PDF PDF 10−3 10−3 t = 0.025 10−4 10−4 t = 0.05 t = 0.075 10−5 10−5 t = 0.1 Gaussian 10−6 −3 −2 −1 0 1 2 3 4 −30 −20 −10 0 10 u0r /σ u0r /σ (a) (b) t = 0.125 10−1 10−1 t = 0.15 t = 0.175 10−2 10−2 t = 0.2 Gaussian PDF PDF 10−3 10−3 t = 0.025 10−4 10−4 t = 0.05 t = 0.075 10−5 t = 0.1 10−5 Gaussian 10−6 −3 −2 −1 0 1 2 3 −4 −2 0 2 4 6 u0θ /σ u0θ /σ (c) (d) Figure 6.8: Time evolution of PDFs of components of the velocity fluctuations at eight instances of time close to the flow instability onset. Here the PDFs are obtained on the entire computational domain while the fluctuations are normalized by their own standard deviations, σ, and they are all compared with the standard Gaussian PDF, N (0, 1). Radial velocity fluctuations: (a) 0 < t ≤ 0.1, (b) 0.1 < t ≤ 0.2. Azimuth velocity fluctuations: (c) 0 < t ≤ 0.1, (d) 0.1 < t ≤ 0.2. from Gaussian statistical behavior that essentially plays an important role in leading to a chaotic flow dynamic state. Here, we compute the velocity fluctuations’ PDFs over the com- putational domain for the radial and azimuth components, and plot them at eight different time states close to the initiation of the flow instability (see Figure 6.8). All of these PDFs are computed for the velocity fluctuations that are normalized by their standard deviation so that the comparison with the standard Gaussian PDF, drawn from N (0, 1), is readily possible through eyeball measure. Here, Figures 6.8a and 6.8c are depicting the PDFs of normalized radial and azimuth components of velocity fluctuations for 0 < t ≤ 0.1, respectively. For both of the radial and azimuth velocity components the PDFs are showing sub-Gaussian 103 2.5 16 t= 0.1 t= 0.1 2.0 t= 0.125 14 t= 0.125 t= 0.15 t= 0.15 1.5 hu0r 3i/hu0r 2i3/2 t= 0.175 12 t= 0.175 hu0r 4i/hu0r 2i2 1.0 t= 0.2 t= 0.2 10 0.5 8 0.0 6 −0.5 4 −1.0 2 −1.5 10−3 10−2 10−1 10−3 10−2 10−1 r r (a) (b) 0.4 t= 0.1 6 0.2 t= 0.125 t= 0.15 hu0θ 3i/hu0θ 2i3/2 0.0 5 t= 0.175 hu0θ 4i/hu0θ 2i2 t= 0.2 −0.2 4 −0.4 t= 0.1 3 −0.6 t= 0.125 −0.8 t= 0.15 2 t= 0.175 −1.0 t= 0.2 10−3 10−2 10−1 10−3 10−2 10−1 r r (c) (d) Figure 6.9: High-order moments of velocity fluctuations, V ′ = (u′r , u′θ ), as a function of radial distance from the wall, r, where r = 0 indicates the wall. (a) Skewness factor for u′r , (b) Flatness factor for u′r , (c) Skewness factor for u′θ , (d) Flatness factor for u′θ . In (b) and (d), the black-colored dashed lines indicate the flatness factor associated with the standard Gaussian distribution. 100 t = 0.25 t = 0.25 10−1 t = 0.5 t = 0.5 10−1 t = 0.75 t = 0.75 Gaussian 10−2 Gaussian 10−2 PDF PDF 10−3 10−3 10−4 10−4 10−5 10−5 10−6 10−6 −20 −10 0 10 20 −4 −2 0 2 4 6 u0r /σ u0θ /σ (a) (b) Figure 6.10: Comparison between the standard Gaussian PDF and PDFs of the velocity fluctuations at t = 0.25, 0.5, 0.75. (a) Radial velocity fluctuations, (b) azimuth velocity fluctuations. 104 behavior that is commonly expected given the laminar initial state of the flow, however, the former rapidly tends to show broader tails compared to the latter with time. Moreover, we can observe that the onset of the flow instability causes noticeable deviations from symmetry in the PDF of radial velocity fluctuations. By tracking the PDFs of velocity fluctuations at further times, i.e. 0.1 < t ≤ 0.2, one can clearly observe that emergence of broad PDF tails and asymmetries quickly leads to a highly non-Gaussian statistical behavior (see Figures 6.8b and 6.8d and compare with the standard Gaussian PDF). More specifically, Figure 6.8b shows that the velocity fluctuations in the radial direction are essentially the main source of this non-Gaussianity as the heavy-tailed PDF accompanied with intermittent events dis- tributed at the PDF tails are arising (see 0.15 ≤ t ≤ 0.2). On the other hand, a noticeable skewness towards the negative-valued fluctuations of the radial velocity component tends to grow with time as shown in Figure 6.8b. Comparing the radial and azimuth components of velocity fluctuations qualitatively show that emerging the aforementioned features that are essentially the fingerprints of non-Gaussian statistics is much milder and at slower rates for the azimuth component, u′θ . In order to obtain a quantitative measure on the non-Gaussian statistics of the velocity fluctuations, we manage to compute their skewness and flatness factors as a function of radial distance from the wall, r. This effectively helps to understand how the non-Gaussian behavior evolves through time as we move away from the wall towards the center. Our approach involves uniformly sampling the velocity values on the circular stripes with a thickness of δr where their radial distance from the wall is r. Once we performed such sampling, we 3 2 4 2 can simply attain the skewness and flatness factors as ⟨V ′ ⟩/⟨V ′ ⟩3/2 and ⟨V ′ ⟩/⟨V ′ ⟩2 , respectively. In our measurements, we took δr = 2 × 10−4 and ⟨·⟩ denotes spatial averaging over the uniformly sampled velocity space on each circular stripe with radial distance r from wall. As a result, Figure 6.9 illustrates such radial skewness and flatness factors for both components of velocity fluctuations at five instances of time for 0.1 ≤ t ≤ 0.2. The resulting measures for u′r depicted in Figures 6.9a and 6.9b show that the non-zero skewness 105 factor and flatness factor greater than 3 (measures associated with standard Gaussian) are appearing for 0.15 ≤ t. This record is in total agreement with what we observe in their non-Gaussian PDFs in Figure 6.8b. For u′θ , Figure 6.9a illustrates non-zero skewness factor values close to the wall at all the recorded times and Figure 6.9b shows that for a narrow region close to the wall the flatness factor exceeds 3 for 0.15 < t. Again, these observations are in complete agreement with the behavior we observe in PDFs of u′θ shown in Figure 6.8d. More specifically on the heavy-tailed velocity fluctuations PDFs, one can link the radial records of flatness factor in both components u′r and u′θ as shown in Figures 6.9b and 6.9d, respectively. In radial velocity fluctuations, it is clearly seen that as time passes the flatness factor increases for the closest radial distances to wall, i.e. r < 10−3 , and in farther distances from the wall, a span of radial region of high flatness factor that essentially contributes to the rare events occurring at the PDF tails (for 0.15 ≤ t) is observed. As we pointed out, this high flatness factor span is expanding towards the center of the cylinder as flow instability evolves in time. Although such behavior is also seen for the azimuth component of velocity fluctuations, its intensity is much milder compared to u′r . In fact, our records show that for u′θ the flatness factor rarely exceeds 3 (see Figure 6.9d). Finally, by comparing the PDFs of velocity fluctuations for 0.2 < t with the one associated with standard Gaussian (see Figure 6.10), we recognize that the statistical features such as non-symmetric distributions and heavy PDF tails with high intermittency are remarkably discernible. However, as illustrated for the prior times closer to the flow instability initiation, these features seem to be manifested more prominently in the radial component of velocity fluctuations. Memory Effects in Vorticity Dynamics and Anomalous Time-Scaling of Enstro- phy Although early theories of Batchelor [260] assumed that for decaying two-dimensional turbulence it is only kinetic energy that is mainly remembered for a long time, later it has been shown that vorticity field plays a key role in the flow dynamics, which was initially 106 100 t = 0.25 t = 0.5 10−1 t = 0.75 Gaussian 10−2 PDF 10−3 10−4 10−5 10−6 −60 −40 −20 0 20 40 60 ωz0 /σ (a) t = 0.25 t = 0.25 10 t = 0.5 400 t = 0.5 t = 0.75 t = 0.75 hωz0 3i/hωz0 2i3/2 5 hωz0 4i/hωz0 2i2 300 0 200 −5 −10 100 −15 0 10−2 10−1 100 10−2 10−1 100 r r (b) (c) Figure 6.11: (a) Comparison between the standard Gaussian PDF and normalized vorticity fluctuations’ PDFs at t = 0.25, 0.5, 0.75, (b) skewness factor for ωz′ and, (c) flatness factor for ωz′ . The dashed lines indicate the measures associated with Gaussian behavior. failed to be addressed by Batchelor [261]. Here, while the filamentation of the vorticity field is occurring, there exist small yet sufficiently strong patches of vorticity surviving the filamentation process and comprise coherent vortices that somehow live even longer than many large-eddy turnover times [2]. These coherent vortices are interacting with each other quite similar to a collection of point vortices. On some occasions, these coherent vortices could approach each other and merge into larger ones. Therefore, the number of coherent vortices decreases while their average size increases as flow evolves. On the other hand, given the discussion on non-Gaussian behavior velocity fluctuations, one can make a con- 107 nection between the statistical behavior of the vorticity field and generation and intensity of coherent vortices resulting from the flow instability. Thus, similar to the procedure in the previous section, we compute the vorticity PDFs in addition to the radial skewness and flatness factors for the same realization of the fluctuating flow field we considered. Figure 6.11 provides this statistical information at t = 0.25, 0.5, and 0.75. Comparing the vorticity PDFs shown in Figure 6.11a to the standard Gaussian PDF makes it evident that finger- prints of non-Gaussian statistics, i.e. non-symmetric probability distributions in addition to broad and intermittent PDF tails, are immensely evolving in vorticity field. Moreover, the radial skewness and flatness factors obtained for these three time instances quantitatively demonstrate that such intense non-Gaussian statistical behavior is swiftly extending towards the center of cylinder (see the radial region of 0.01 < r < 0.2 at Figures 6.11b and 6.11c). Given the discussion on the generation and evolution of the coherent vortices, and our quantitative/qualitative study on the emergence of strong non-Gaussian statistical behavior for velocity and vorticity fluctuations, one can argue that such statistics are closely tied to and in other words, the direct result of generation and growth of coherent vortical structures due to the effect of the rotational symmetry-breaking factors. In prior studies, such connection was investigated and partially addressed in the contexts of planar mixing and free shear layers [160, 262, 263], subgrid-scale (SGS) motions and their nonlocal modeling for wall- bounded turbulent flows, boundary layer flows [264, 265], and turbulent flows interacting with wavy-like moving/actuated surfaces (with application to reduction and control of flow separation) [266, 267]. Here, an interesting yet, practical question that could be raised is that if such “intensified” coherent vortical structures induced by the symmetry-breaking parameters in the rotational motion are capable of incorporating more memory effects into the dynamics of vorticity field. This potentially could lead to the engineering means to increase effective chaotic mixing in rotating systems by introducing factors that initiate deviation from symmetry in rotation. In a two-dimensional turbulent/chaotic flow, the very presence of “long-lived” coherent vortices 108 105 E(t) ∼ t−2 (I) ∼ t−1 ∼ t−1/2 Enstrophy 104 (III) 103 (II) 102 100 101 t (a) (b) Figure 6.12: Time-scaling of enstrophy record and its link to evolution of coherent vortical structures. (a) Enstrophy record, E(t), and its early-time (I), transient-time (II), and long- time (III) scaling affected by the imposed symmetry-breaking disturbances on the rotational motion of cylinder, (b) snapshots of instantaneous vorticity field, ωz (x, t), (left) t = 7 and (right) t = 24 showing the structure and growth of coherent vortical regions attached to the cylinder wall. normally cause the time-scaling of enstrophy record at long-time to be close to t−1 , however, it initially is scaled with t−2 at the early stages of flow which is also what Batchelor’s theory envisions [2]. Therefore, a relevant approach to seek an answer to this question is to study the long-time behavior of enstrophy record that contains the spatially integrated information in the vortical motions over the entire domain and also is a representative for the dissipation dynamics. Similar to the kinetic energy (6.16), we define the enstrophy, E(t), in our problem setting as Z 1 2 E(t) = ωz′ (x, t) dΩ. (6.20) µ(Ω) Ω By computing the record of enstrophy for relatively long times (obtained from the same flow 109 realization we studied its fluctuating velocity and vorticity behavior), studying the early- /long-time scaling trend of enstrophy would be possible. To perform this very study, the validity and stability of long-time evaluation of QoIs for stochastic mathematical models is of crucial importance to be considered and it has been addressed in multiple prior studies. For instance, Xiu and Karniadakis [238] used generalized polynomial chaos (gPC) with relatively high resolutions in order to study the long-time behavior of vorticity field for the flow past a cylinder under the uncertain inflow boundary conditions. In another study, Xiu and Hesthaven [246] employed high-order stochastic collocation methods to achieve stable second- order moment response to the stochastic differential equations at the long times. Moreover, Foo et al. [248] utilized multi-element probabilistic collocation method (ME-PCM) with high resolution in random space to compute stable long-time flow records. Therefore, maintaining sufficiently high resolutions in discretization of random space is a key point. In our study, the high-resolution uni-variate PCM we employed to obtain the fluctuating flow fields (as described in section 6.6) essentially guarantees the validity and statistical stability of our evaluations for the long-time fluctuating vortictiy field and computing the enstrophy record as illustrated in Figure 6.12a. This plot shows that in terms of enstrophy time-scaling, we observe three stages of time. Here at stage (I), enstrophy behaves as E ∼ t−2 (for t < 2.5), however, after a transition period, stage (II), it persistently follows E ∼ t−1/2 time-scaling in stage (III). At the third stage, this “anomalous” long-time scaling with t−1/2 rather than the expected t−1 scaling could essentially be interpreted as the result of an “intensified” mechanism for birth and growth of coherent vortices that live for effectively long periods of time during the evolution of this internal flow right after the occurrence of the flow instability. Figure 6.12b portrays two snapshots of instantaneous vorticity field, ωz , on a segment of cylinder close to the wall to show the evolution and form of these coherent vortex structures survived the vortex filamentation process. We emphasize that the long life of the mature and relatively large coherent vortical zones (clearly visible and attached to the cylinder’s wall) is the main reason of the anomalous enstrophy time-scaling we observe at 110 stage (III) in Figure 6.12a. As we mentioned earlier, this phenomenon could potentially be a practical engineering candidate to enhance and reinforce the effective chaotic/turbulent mixing qualities by inducing more memory effects resulted from a symmetry-breaking flow instability. 111 CHAPTER 7 SUMMARY AND FUTURE WORKS The objective of this dissertation was to develop robust modeling strategies utilizing data- driven methods, nonlocal mathematical operators, and high-performance computing for sim- ulations and analysis of turbulent mixing phenomena. In the following sections, we point out the concluding remarks for each chapter of this study, and later we provide some outlook to the future directions. 7.1 Concluding Remarks • In chapter 2, we presented a computational platform for DNS of homogeneous turbulent flow and passive scalar transport. This open-source software works based upon a pseudo- spectral representation of the NS and AD equations on a triply cubic computational domain with periodic boundary conditions for fluctuating fields. Using Fourier collocation method, the governing equations are discretized in space and by employing RK4 scheme the time- stepping is performed. The software provides a pre-processing step to construct homogeneous and isotropic divergence-free velocity IC based on prescribed energy spectrum and decom- pose it into user-defined partitions. Using artificial forcing scheme, the dissipated energy is injected to the low wavenumbers so that after long-time integration, the statistically station- ary state is achieved. In order to examine and identify if the fully-developed turbulent flow is obtained, small-scale statistical quantities of turbulence in addition to the central moments of VGT components are computed and recorded. Once the realistic turbulent velocity field is obtained, the user is able to start resolving a passive scalar that is transported on the HIT flow while a uniform mean scalar gradient is imposed. Resolving the scalar fluctuations for long enough time after reaching to the equilibrium and stationary state provides the fully- developed turbulent scalar field. Statistical records of scalar gradients in addition to the records of production and dissipation of scalar variance helps the user to properly identify when the fully-developed scalar turbulence is achieved. • In chapter 3, we developed a new data-driven nonlocal/fractional SGS model for the 112 LES of passive scalars transported in the homogeneous isotropic turbulent flow. The main focus of our work was on obtaining an SGS model that is structurally designed based on the nonlocal nature of the SGS scalar flux. Therefore, we first managed to present a through statistical interpretation of nonlocality in the SGS dynamics using the single- and two- point statistics of the SGS scalar dissipation. Using a rich dataset of high-fidelity data for the SGS flux obtained from direct filtering of DNS results, we illustrated the statistical nonlocality embedded in the SGS dynamics and showed that it amplifies as the filter-width increases. Moreover, we showed that the conventional means of SGS modeling originate from a local statistical representation for the SGS dynamics and are intrinsically incapable of predicting the statistical nonlocality. As a robust starting point for our mathematical modeling, we started from Boltzmann-BGK kinetics as the microscopic transport framework for passive scalars in homogeneous turbulence and considered the closure problem manifested in filtering the transport equations. By revisiting the kinetic-level strategy for the LES modeling taking into account the consistency of the the model for the filtered equilibrium distribution with its macroscopic representation at the continuum level, we proposed to proceed with closure modeling using α-stable Lévy distribution to address the nonlocal and non-Gaussian behavior of the closure at the kinetic level. In order to derive a macroscopic representation of such model to employ in the filtered AD equation, we used continuum averaging and obtained the filtered and residual (modeled) passive scalar flux components that essentially return the filtered AD equation. Throughout this procedure, the up-scaled model for the divergence of the residual flux takes the form of a fractional Laplacian acting on the filtered scalar concentration with a model-specific proportionality coefficient. Next, we managed to calibrate the fractional-order model in two separate data-driven stages. First we targeted identification of the optimal fractional order using two-point statistics data for the normalized SGS dissipation function obtained from the DNS and minimizing the mismatch function with its counterpart in the fractional-order SGS model. This procedure returned the optimal fractional order that minimizes the single-point correlation between 113 the modeled and true SGS scalar flux. Afterwards, following an sparse regression strategy over the spatio-temporal data for the SGS scalar flux in a statistically-stationary turbulent scalar field, we obtained the proportionality coefficient of the model. Moreover, we showed the consistency of the derived model in terms of the relationship between the obtained proportionality coefficient and decreasing the filter-width. Finally, in an a priori test, we showed that the identified model is capable of capturing the PDF tail associated with the forward scattering of the filtered scalar variance and illustrated that our model has the capability to partially reproduce the backward scattering phenomenon. • In chapter 4, we proposed a modification to the spectral transfer model for the turbulent cascade of passive scalars under the effect of large-scale anisotropy. Employing the Corrsin’s generalization to Onsager’s turbulent cascade model, our modified model introduced an ad- ditional power-law term in the definition of local time-scale, τ (k), in order to account for the induced nonlocal contributions originated from the anisotropy sources in the energy contain- ing range. Subsequently, our approach yielded a modified scaling law for the passive scalar spectrum, Eϕ (k). This modified scaling showed a great match with the time-averaged 3-D scalar spectrum obtained from a well-resolved standard DNS after the parameter identifi- cation procedure. Using the integrated equation for the evolution of scalar spectrum, we revised the total scalar dissipation definition, which introduced an additional term into the total scalar dissipation representing the integrated effects of nonlocal turbulent dissipation cascade. This modification to the scalar dissipation returned that a fractional-order Lapla- cian acting on the scalar concentration is required in the AD equation. Using this revised AD equation, we performed a DNS study to analyze different quantities of turbulent transport compared to the DNS results obtained from convectional model at the statistical equilibrium state. Our analysis on the rate of the scalar variance showed that considering the effects of nonlocality in the scalar dissipation results in pronounced prediction of the production rate of scalar variance by the imposed mean-gradient (large-scale anisotropy), which could be interpreted in consistency with the breakdown of local isotropy in small scales of passive 114 scalars. On the other hand, we showed that incorporation of the nonlocality effects in the scalar dissipation improves the accuracy of predicted time-averaged records of the skewness and flatness factors for the ∇∥ ϕ, confirming the essence of devising a proper modeling mech- anism for cascade of the anisotropy effects from large to small scales. Moreover, a two-point statistical analysis for the advective scalar increments (the third-order mixed longitudinal structure function) revealed that the DNS results obtained from nonlocal model provides a long-range scaling with r2 . This observation on long-range scaling suggested that the inclu- sion of nonlocal cascading mechanism in the presence of large-scale anisotropy could result in prediction of more universal behavior over a wide span of scales in turbulent scalar trans- port. Finally, we showed an accurate consistency between the developed spectral transfer model and the fractional-order SGS modeling with ∆ ≈ 2 η, after employing a well-trained GPR model. • In chapter 5, we developed a data-driven surrogate model for prediction of the closure flux term appearing in the LES equations of a passive scalar. The developed model was a fully connected deep neural network architecture that takes gradient of the filtered/resolved scale transport fields as input features. In order to provide a rich set of training/validation data from FDNS, we introduced a spatio-temporal procedure to sufficiently sample from a in fully-turbulent time domain and also the homogeneous spatial turbulent domain at each sampled time frame. Using dataset, we managed to optimize the model parameters of our specified DNN model with respect to the MSE loss function. Eventually, we presented a proper optimized DNN model where training and validation loss records have shown perfectly matching behavior, therefore, the any over-fitting was avoided. Using the notion of transfer learning, we efficiently (in terms of data requirement and Computational cost) generalized our initially trained model into a case for a higher Schmidt number, and another case for a higher Péclet number. In training the TL models, we ensured that we are avoiding an over- fitting after the re-training procedure. Afterwards, we managed to test the performance of the developed DNN models in the inference mode. Therefore, in an a prori assessment, 115 we compared the performance of the base DNN model with the PSM and SAG models at two filter scales, and realized that the DNN model not only reduces the level of error in predictions of SGS flux and dissipation, but also improves the consistency in prediction of SGS flux with respect to the direction of large-scale anisotropy, ∇⟨Φ̃⟩ compared to the traditional SGS models. For the TL models, we tested the pre-trained model’s performance in prediction of a new scalar transport regime against the re-trained models and found out that re-training significantly enhances the over-prediction of the forward scattering phenomenon (and even the backward scattering of filtered variance at higher Péclet number regime). Finally, along the SAG and PSM model, we tested our base DNN model in an LES setting. We compared the records of resolved-scale scalar variance from the employed SGS model as well as the second-order and third-order structure functions of the resolved scalar. Compared to the traditional models we tested, we realized that a well-trained DNN model can result in striking error reduction in the recorded scalar variance as well as significant improvement in prediction of the two-point structure functions. • Finally, in chapter 6, our study leverages the outcome of stochastic modeling and simula- tions to carry out a thorough analysis on the initiation of flow instabilities within high-speed rotating cylinders. Considering the random nature of the problem, a detailed mathematical representation of the stochastic incompressible Navier-Stokes equations was presented. Fur- ther, a high-fidelity stochastic CFD framework was introduced, which employs spectral/hp element method in the forward solver and later on the stochastic space was numerically handled by probabilistic collocation method. Detailed grid generation steps and required convergence studies for the deterministic solver were obtained and stochastic discretization convergence were studied for the solutions of first and second moments. The time-evolution of expected kinetic energy of the flow in addition to its variance were computed and the uncertainty bounds propagated in the solution were identified with time. A variance-based sensitivity analysis of the random parameters of the model were conducted to globally char- acterize the most effective stochastic factor on the total variance of kinetic energy, conse- 116 quently, the “eccentric rotation” was learned to be the dominant source of stochasticity. Later on, the expected solution from a very fine uni-variate PCM discretization on the dom- inant random parameter was utilized to compute the fluctuating velocity and vorticity fields for a randomly drawn realization of the sample space. These fluctuations were statistically analyzed through the time-evolution of their PDFs for radial and azimuth components in a qualitative manner while comparing to the standard Gaussian PDF. Statistical features such as appearance of intermittent and rare events in terms of heavy-tailed PDFs in addition to observing asymmetries in velocity and vorticity PDFs were spotted out. In particular, very close to the flow instability onset, these non-Gaussian statistical features were found to quickly get intensified especially for the radial velocity fluctuations and therefore fluctuating vorticity field as the flow evolves in time. Moreover, the statistics of flow fields were quan- titatively measured through computing the skewness and flatness factors on narrow radial stripes extending from the wall to the cylinder’s center. These records closely supported our qualitative findings from studying the PDFs of fluctuations and identified that in velocity field we quickly face regions with skewness factor of O(1) and flatness factor of O(10) while for the vorticity field these factors were recorded with about one order of magnitude higher than their velocity counterparts emphasizing on significantly high non-Gaussian vorticity induced by cylinder rotation affected by symmetry-breaking factors. Here, we need to em- phasize that the modeling strategy in this work was conducted based upon an approximate 2-D representation of the NS equations for long rotating cylinders fully-filled with fluid. It properly provided a valuable insight into the dominant effect that breaks the flow symmetry and the dynamics of the further fluctuations in the flow variables. Certainly, a similar study focusing on the short-height rotating flow systems requires three-dimensional modeling and simulations. In fact, since the top roof and bottom base surfaces of the cylinder additionally have major contributions to the symmetry-breaking instabilities, the axial propagation of uncertainty in the flow would be inevitable and must be taken into consideration. Motivated by the observed strong non-Gaussianity in flow fluctuations, we sought to study the effects of 117 coherent vortical structures essentially inducing memory effects into the vorticity dynamics. Thus, we managed to compute the time-scaling of the enstrophy record. Interestingly, we learned that unlike the early stages of flow after introduction symmetry-breaking rotational effects, enstrophy is scaled as t−1/2 at long-time. This anomalous time-scaling essentially reveals the very existence of long-lasting and growing coherent vortical regions initially gen- erated due to the non-symmetric rotation of the cylinder wall. This mechanism seems to be a promising engineering strategy to increase the chaotic/turbulent mixing time and quality for the rotating hydrodynamic systems. 7.2 Future Directions and Outlook Some open topics remain in this study to be addressed in the future works. Here, we list some of them as follows: • Testing the SGS model on other types of turbulent flow: investigating the charac- teristics of the fractional SGS model in different types of turbulent flows such as shear flows and wall-bounded ones to uncover the trend in the parameter estimation required to capture the nonlocality in each regime, followed by a comprehensive statistical comparison between model performance at each regime. • Other deep learning modeling strategies: developing other deep learning SGS models for the turbulent transport by embedding the relevant features in the surrogate model, and different model architectures such as convolutional neural networks, transformer networks with attention mechanisms, and reinforcement learning approaches. • Learning the SGS dynamics from a Lagrangian point of view: employing a parti- cle tracking approach over the fully-developed turbulent fields, the nonlocal behavior of the subgrid-scale fluctuations is aimed to be learned from analysis of time-series through sophis- ticated deep learning algorithms. After integration of the learned statistical information, discovery of the proper differential operators to govern the evolution of the SGS dynamics is pursued. 118 BIBLIOGRAPHY [1] R. Klages, G. Radons, and I. M Sokolov. Anomalous transport: foundations and applications. John Wiley & Sons, 2008. [2] Peter Alan Davidson. Turbulence: an introduction for scientists and engineers. Oxford university press, 2015. [3] Peter William Egolf and Kolumban Hutter. Nonlinear, Nonlocal and Fractional Tur- bulence: Alternative Recipes for the Modeling of Turbulence. Springer Nature, 2019. [4] Stephen B Pope. Turbulent flows, 2001. [5] Pierre Sagaut and Claude Cambon. Homogeneous turbulence dynamics, volume 10. Springer, 2008. [6] Themistoklis P Sapsis. Statistics of extreme events in fluid flows and waves. Annual Review of Fluid Mechanics, 53. [7] PK Yeung, XM Zhai, and Katepalli R Sreenivasan. Extreme events in computational turbulence. Proceedings of the National Academy of Sciences, 112(41):12633–12638, 2015. [8] Ali Akhavan-Safaei and Mohsen Zayernouri. A Parallel Integrated Computational- Statistical Platform for Turbulent Transport Phenomena. arXiv preprint arXiv:2012.04838, 2020. [9] Lin Fu, Sanjeeb Bose, and Parviz Moin. Heat transfer in three-dimensional intersect- ing shock-wave/turbulent boundary-layer interactions with wall-modeled large-eddy simulations. arXiv preprint arXiv:2009.02411, 2020. [10] Athony Leonard. Energy cascade in large-eddy simulations of turbulent fluid flows. In Advances in geophysics, volume 18, pages 237–248. Elsevier, 1975. [11] Massimo Germano. Turbulence: the filtering approach. Journal of Fluid Mechanics, 238:325–336, 1992. [12] Charles Meneveau. Statistics of turbulence subgrid-scale stresses: Necessary conditions and experimental tests. Physics of Fluids, 6(2):815–833, 1994. [13] Robert D Moser, Sigfried W Haering, and Gopal R Yalla. Statistical properties of subgrid-scale turbulence models. Annual Review of Fluid Mechanics, 53, 2021. [14] Fabian Waleffe. The nature of triad interactions in homogeneous turbulence. Physics of Fluids A: Fluid Dynamics, 4(2):350–363, 1992. [15] Peter E Hamlington, Jörg Schumacher, and Werner JA Dahm. Local and nonlo- cal strain rate fields and vorticity alignment in turbulent flows. Physical Review E, 77(2):026303, 2008. 119 [16] Mehdi Samiee, Ali Akhavan-Safaei, and Mohsen Zayernouri. A fractional subgrid- scale model for turbulent flows: Theoretical formulation and a priori study. Physics of Fluids, 32(5):055102, 2020. [17] Mehdi Samiee, Ali Akhavan-Safaei, and Mohsen Zayernouri. Tempered Fractional LES Modeling. Journal of Fluid Mechanics, 932, 2022. [18] Mehdi Samiee. Data-Infused Fractional Modeling and Spectral Numerical Analysis for Anomalous Transport and Turbulence. Michigan State University, 2021. [19] Ali Akhavan-Safaei and Mohsen Zayernouri. A Nonlocal Spectral Transfer Model and New Scaling Law for Scalar Turbulence. arXiv preprint arXiv:2111.06540, 2021. [20] A. N. Kolmogorov. Energy dissipation in locally isotropic turbulence. In Dokl. Akad. Nauk. SSSR, volume 32, pages 19–21, 1941. [21] A. N. Kolmogorov. Local structure of turbulence in incompressible fluid under very high values of reynolds numbers. Reports of AS USSR, 30(4):299–303, 1941. [22] L. Onsager. The distribution of energy in turbulence. In Physical Review, volume 68, pages 286–286, 1945. [23] L. Onsager. Statistical hydrodynamics. Il Nuovo Cimento, 6(2):279–287, March 1949. [24] A. M. Obukhov. Structure of temperature field in turbulent flow. In Izv. Akad. Nauk. SSSR, Georg. i Geofiz., volume 13, pages 58–69, 1949. [25] A. M. Yaglom. On the local structure of a temperature field in a turbulent flow. In Dokl. Akad. Nauk. SSSR, volume 69, 1949. [26] Stanley Corrsin. On the spectrum of isotropic temperature fluctuations in an isotropic turbulence. Journal of Applied Physics, 22(4):469–473, 1951. [27] George K Batchelor. Small-scale variation of convected quantities like temperature in turbulent fluid Part 1. General discussion and the case of small conductivity. Journal of Fluid Mechanics, 5(1):113–133, 1959. [28] GK Batchelor, ID Howells, and AA Townsend. Small-scale variation of convected quantities like temperature in turbulent fluid Part 2. The case of large conductivity. Journal of Fluid Mechanics, 5(1):134–139, 1959. [29] Andrey Nikolaevich Kolmogorov. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number. Journal of Fluid Mechanics, 13(1):82–85, 1962. [30] A. M. Oboukhov. Some specific features of atmospheric tubulence. Journal of Fluid Mechanics, 13(1):77–81, 1962. [31] A. S. Monin and A. M. Yaglom. Statistical fluid mechanics, volume II: mechanics of turbulence, volume 2. MIT press, 1975. 120 [32] Dhawal Buaria, Matthew P Clay, Katepalli R Sreenivasan, and PK Yeung. Small- scale isotropy and ramp-cliff structures in scalar turbulence. Physical Review Letters, 126(3):034504, 2021. [33] Alain Pumir and Boris I. Shraiman. Persistent Small Scale Anisotropy in Homogeneous Shear Flows. Phys. Rev. Lett., 75:3114–3117, Oct 1995. [34] Katepalli R Sreenivasan. Turbulent mixing: A perspective. Proceedings of the National Academy of Sciences, 116(37):18175–18183, 2019. [35] Kerstin Avila, David Moxey, Alberto de Lozar, Marc Avila, Dwight Barkley, and Björn Hof. The onset of turbulence in pipe flow. Science, 333(6039):192–196, 2011. [36] JH Gerrard. The mechanics of the formation region of vortices behind bluff bodies. Journal of fluid mechanics, 25(2):401–413, 1966. [37] GH Koopmann. The vortex wakes of vibrating cylinders at low reynolds numbers. Journal of Fluid Mechanics, 28(3):501–512, 1967. [38] Owen M Griffin and Steven E Ramberg. Vortex shedding from a cylinder vibrating in line with an incident uniform flow. Journal of Fluid Mechanics, 75(2):257–271, 1976. [39] C Barbi, DP Favier, CA Maresca, and DP Telionis. Vortex shedding and lock-on of a circular cylinder in oscillatory flow. Journal of Fluid Mechanics, 170:527–544, 1986. [40] A Ongoren and D Rockwell. Flow structure from an oscillating cylinder part 1. mech- anisms of phase shift and recovery in the near wake. Journal of fluid Mechanics, 191:197–223, 1988. [41] A Ongoren and D Rockwell. Flow structure from an oscillating cylinder part 2. mode competition in the near wake. Journal of Fluid Mechanics, 191:225–245, 1988. [42] SCR Dennis, P Nguyen, and Serpil Kocabiyik. The flow induced by a rotationally oscillating and translating circular cylinder. Journal of Fluid Mechanics, 407:123–144, 2000. [43] MRH Nobari and H Naderan. A numerical study of flow past a cylinder with cross flow and inline oscillation. Computers & Fluids, 35(4):393–415, 2006. [44] E Crespo Del Arco, Eric Serre, Patrick Bontoux, and BE Launder. Stability, transition and turbulence in rotating cavities. Advances in Fluid Mechanics, 41:141–196, 2005. [45] Marcello Lappa. Rotating thermal flows in natural and industrial processes. John Wiley & Sons, 2012. [46] Min Chan Kim and Chang Kyun Choi. The onset of Taylor–Görtler vortices during impulsive spin-down to rest. Chemical engineering science, 61(19):6478–6485, 2006. [47] Alban Sauret, David Cébron, Michael Le Bars, and Stéphane Le Dizès. Fluid flows in a librating cylinder. Physics of Fluids, 24(2):026603, 2012. 121 [48] Frieder Kaiser, Bettina Frohnapfel, Rodolfo Ostilla-Mónico, Jochen Kriegseis, David E. Rival, and Davide Gatti. On the stages of vortex decay in an impulsively stopped, rotating cylinder. Journal of Fluid Mechanics, 885:A6, 2020. [49] J. P. Gollub and Harry L. Swinney. Onset of turbulence in a rotating fluid. Phys. Rev. Lett., 35:927–930, Oct 1975. [50] B. Dubrulle, O. Dauchot, F. Daviaud, P.-Y. Longaretti, D. Richard, and J.-P. Zahn. Stability and turbulent transport in Taylor–Couette flow from analysis of experimental data. Physics of Fluids, 17(9):095103, 2005. [51] Yufeng Lin, Jerome Noir, and Andrew Jackson. Experimental study of fluid flows in a precessing cylindrical annulus. Physics of Fluids, 26(4):046604, 2014. [52] Romain Lagrange, Christophe Eloy, François Nadal, and Patrice Meunier. Instability of a fluid inside a precessing cylinder. Physics of Fluids, 20(8):081701, 2008. [53] Rodolfo Ostilla, Richard J. A. M. Stevens, Siegfried Grossmann, Roberto Verzicco, and Detlef Lohse. Optimal Taylor–Couette flow: direct numerical simulations. Journal of Fluid Mechanics, 719:14–46, 2013. [54] Hao Teng, Nansheng Liu, Xiyun Lu, and Bamin Khomami. Direct numerical simulation of Taylor-Couette flow subjected to a radial temperature gradient. Physics of Fluids, 27(12):125101, 2015. [55] Christopher J. Crowley, Michael C. Krygier, Daniel Borrero-Echeverry, Roman O. Grigoriev, and Michael F. Schatz. A novel subcritical transition to turbulence in Taylor–Couette flow with counter-rotating cylinders. Journal of Fluid Mechanics, 892:A12, 2020. [56] JM Lopez, JE Hart, F Marques, S Kittelman, and J Shen. Instability and mode interactions in a differentially driven rotating cylinder. Journal of Fluid Mechanics, 462:383, 2002. [57] J. M. Lopez and F. Marques. Instabilities and inertial waves generated in a librating cylinder. Journal of Fluid Mechanics, 687:171–193, 2011. [58] Juan M. Lopez and Francisco Marques. Sidewall boundary layer instabilities in a rapidly rotating cylinder driven by a differentially corotating lid. Physics of Fluids, 22(11):114109, 2010. [59] Juan M. Lopez and Francisco Marques. Rapidly rotating precessing cylinder flows: forced triadic resonances. Journal of Fluid Mechanics, 839:239–270, 2018. [60] Karianne J Bergen, Paul A Johnson, Maarten V de Hoop, and Gregory C Beroza. Machine learning for data-driven discovery in solid earth geoscience. Science, 363(6433):eaau0323, 2019. 122 [61] Noah D Brenowitz and Christopher S Bretherton. Prognostic validation of a neural network unified physics parameterization. Geophysical Research Letters, 45(12):6289– 6298, 2018. [62] Noah D Brenowitz and Christopher S Bretherton. Spatially extended tests of a neural network parametrization trained by coarse-graining. Journal of Advances in Modeling Earth Systems, 11(8):2728–2744, 2019. [63] J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor. Dynamic mode decomposition: data-driven modeling of complex systems. SIAM, 2016. [64] Paul Raccuglia, Katherine C Elbert, Philip DF Adler, Casey Falk, Malia B Wenny, Aurelio Mollo, Matthias Zeller, Sorelle A Friedler, Joshua Schrier, and Alexander J Norquist. Machine-learning-assisted materials discovery using failed experiments. Na- ture, 533(7601):73–76, 2016. [65] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neu- ral networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. [66] Siddhant Kumar, Stephanie Tan, Li Zheng, and Dennis M Kochmann. Inverse-designed spinodoid metamaterials. npj Computational Materials, 6(1):1–10, 2020. [67] Amirhossein Arzani, Jian-Xun Wang, and Roshan M D’Souza. Uncovering near-wall blood flow from sparse data with physics-informed neural networks. Physics of Fluids, 33(7):071905, 2021. [68] Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differen- tial equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018. [69] Justin Sirignano and Konstantinos Spiliopoulos. DGM: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339–1364, 2018. [70] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equa- tions from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932–3937, 2016. [71] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Data-driven discovery of partial differential equations. Science advances, 3(4):e1602614, 2017. [72] Kai Fukami, Koji Fukagata, and Kunihiko Taira. Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics, 870:106–120, 2019. [73] Jin-Long Wu, Karthik Kashinath, Adrian Albert, Dragos Chirila, Heng Xiao, et al. Enforcing statistical constraints in generative adversarial networks for modeling chaotic dynamical systems. Journal of Computational Physics, 406:109209, 2020. 123 [74] Matthew Li and Christopher McComb. Using Physics-Informed Generative Adversarial Networks to Perform Super-Resolution for Multiphase Fluid Simulations. Journal of Computing and Information Science in Engineering, 22(4):044501, 2022. [75] J Nathan Kutz. Deep learning in fluid dynamics. Journal of Fluid Mechanics, 814:1–4, 2017. [76] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos. Machine learning for fluid mechanics. Annual Review of Fluid Mechanics, 52:477–508, 2020. [77] Jin-Long Wu, Heng Xiao, and Eric Paterson. Physics-informed machine learning ap- proach for augmenting turbulence models: A comprehensive framework. Physical Re- view Fluids, 3(7):074602, 2018. [78] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao. Turbulence modeling in the age of data. Annual Review of Fluid Mechanics, 51:357–377, 2019. [79] Karthik Duraisamy. Perspectives on machine learning-augmented reynolds-averaged and large eddy simulation models of turbulence. Physical Review Fluids, 6(5):050504, 2021. [80] Andrea Beck and Marius Kurz. A perspective on machine learning methods in turbu- lence modeling. GAMM-Mitteilungen, 44(1):e202100002, 2021. [81] Claudia Drygala, Benjamin Winhart, Francesca di Mare, and Hanno Gottschalk. Gen- erative modeling of turbulence. Physics of Fluids, 34(3):035114, 2022. [82] Wenhui Peng, Zelong Yuan, and Jianchun Wang. Attention-enhanced neural network models for turbulence simulation. Physics of Fluids, 34(2):025111, 2022. [83] Julia Ling, Andrew Kurzawski, and Jeremy Templeton. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics, 807:155–166, 2016. [84] Jian-Xun Wang, Jin-Long Wu, and Heng Xiao. Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data. Physical Review Fluids, 2(3):034603, 2017. [85] Rozie Zangeneh. Data-driven model for improving wall-modeled large-eddy simulation of supersonic turbulent flows with separation. Physics of Fluids, 33(12):126103, 2021. [86] Jin-Ping Li, Deng-Gao Tang, Chen Yi, and Chao Yan. Data-augmented turbulence modeling by reconstructing reynolds stress discrepancies for adverse-pressure-gradient flows. Physics of Fluids, 34(4):045110, 2022. [87] Alexis-Tzianni G Charalampopoulos and Themistoklis P Sapsis. Machine-learning energy-preserving nonlocal closures for turbulent fluid flows and inertial tracers. Phys- ical Review Fluids, 7(2):024305, 2022. 124 [88] Pierre Sagaut. Large eddy simulation for incompressible flows: an introduction. Springer Science & Business Media, 2006. [89] Fabrizio Sarghini, G De Felice, and Stefania Santini. Neural networks based subgrid scale modeling in large eddy simulations. Computers & fluids, 32(1):97–108, 2003. [90] Andrea Beck, David Flad, and Claus-Dieter Munz. Deep neural networks for data- driven LES closure models. Journal of Computational Physics, 398:108910, 2019. [91] Chenyue Xie, Jianchun Wang, and E Weinan. Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence. Physical Review Fluids, 5(5):054606, 2020. [92] Zelong Yuan, Chenyue Xie, and Jianchun Wang. Deconvolutional artificial neural net- work models for large eddy simulation of turbulence. Physics of Fluids, 32(11):115106, 2020. [93] Yunpeng Wang, Zelong Yuan, Chenyue Xie, and Jianchun Wang. Artificial neural network-based spatial gradient models for large-eddy simulation of turbulence. AIP Advances, 11(5):055216, 2021. [94] Josh Williams, Uwe Wolfram, and Ali Ozel. Neural stochastic differential equations for particle dispersion in large-eddy simulations of homogeneous isotropic turbulence. Physics of Fluids, 2022. [95] Masataka Gamahara and Yuji Hattori. Searching for turbulence models by artificial neural network. Physical Review Fluids, 2(5):054604, 2017. [96] Jonghwan Park and Haecheon Choi. Toward neural-network-based large eddy simula- tion: Application to turbulent channel flow. Journal of Fluid Mechanics, 914, 2021. [97] Bo Liu, Huiyang Yu, Haibo Huang, Nansheng Liu, and Xiyun Lu. Investigation of non- local data-driven methods for subgrid-scale stress modeling in large eddy simulation. AIP Advances, 12(6):065129, 2022. [98] H Jane Bae and Petros Koumoutsakos. Scientific multi-agent reinforcement learning for wall-models of turbulent flows. Nature Communications, 13(1):1–9, 2022. [99] Junhyuk Kim, Hyojin Kim, Jiyeon Kim, and Changhoon Lee. Deep reinforcement learning for large-eddy simulation modeling in wall-bounded turbulence. arXiv preprint arXiv:2201.09505, 2022. [100] XIA Yang, S Zafar, J-X Wang, and H Xiao. Predictive large-eddy-simulation wall modeling via physics-informed neural networks. Physical Review Fluids, 4(3):034602, 2019. [101] Justin Sirignano, Jonathan F MacArt, and Jonathan B Freund. DPM: A deep learning PDE augmentation method with application to large-eddy simulation. Journal of Computational Physics, 423:109811, 2020. 125 [102] Jonathan F MacArt, Justin Sirignano, and Jonathan B Freund. Embedded training of neural-network subgrid-scale turbulence models. Physical Review Fluids, 6(5):050502, 2021. [103] Hugo Frezat, Guillaume Balarac, Julien Le Sommer, Ronan Fablet, and Redouane Lguensat. Physical invariance in neural networks for subgrid-scale scalar flux modeling. Physical Review Fluids, 6(2):024607, 2021. [104] Igor Podlubny. Fractional differential equations: an introduction to fractional deriva- tives, fractional differential equations, to methods of their solution and some of their applications, volume 198. Academic press, 1998. [105] Mark M Meerschaert and Alla Sikorskii. Stochastic models for fractional calculus, volume 43. Walter de Gruyter, 2011. [106] Mehdi Samiee, Ehsan Kharazmi, and Mohsen Zayernouri. Fast spectral methods for temporally-distributed fractional PDEs. In Spectral and High Order Methods for Par- tial Differential Equations ICOSAHOM 2016, pages 651–667. Springer, 2017. [107] Mehdi Samiee, Mohsen Zayernouri, and Mark M Meerschaert. A unified spectral method for FPDEs with two-sided derivatives; part I: a fast solver. Journal of Com- putational Physics, 385:225–243, 2019. [108] Mehdi Samiee, Mohsen Zayernouri, and Mark M Meerschaert. A unified spectral method for FPDEs with two-sided derivatives; part II: Stability, and error analysis. Journal of Computational Physics, 385:244–261, 2019. [109] Anna Lischke, Mohsen Zayernouri, and Zhongqiang Zhang. Spectral and spectral ele- ment methods for fractional advection–diffusion–reaction equations. Numerical Meth- ods, page 157, 2019. [110] Mehdi Samiee, Ehsan Kharazmi, Mark M Meerschaert, and Mohsen Zayernouri. A Uni- fied Petrov–Galerkin Spectral Method and Fast Solver for Distributed-Order Partial Differential Equations. Communications on Applied Mathematics and Computation, pages 1–30, 2020. [111] Yongtao Zhou, Jorge L Suzuki, Chengjian Zhang, and Mohsen Zayernouri. Implicit- explicit time integration of nonlinear fractional differential equations. Applied Numer- ical Mathematics, 165:555–583, 2020. [112] Ning Du, Xu Guo, and Hong Wang. Fast upwind and Eulerian-Lagrangian control volume schemes for time-dependent directional space-fractional advection-dispersion equations. Journal of Computational Physics, 405:109127, 2020. [113] Jorge L Suzuki and Mohsen Zayernouri. A self-singularity-capturing scheme for frac- tional differential equations. International Journal of Computer Mathematics, pages 1–28, 2020. 126 [114] Ehsan Kharazmi and Mohsen Zayernouri. Fractional sensitivity equation method: Ap- plication to fractional model construction. Journal of Scientific Computing, 80(1):110– 140, 2019. [115] Ehsan Kharazmi and Mohsen Zayernouri. Operator-based uncertainty quantification of stochastic fractional partial differential equations. Journal of Verification, Validation and Uncertainty Quantification, 4(4), 2019. [116] Jorge L Suzuki, Mamikon Gulian, Mohsen Zayernouri, and Marta D’Elia. Fractional modeling in action: A survey of nonlocal models for subsurface transport, turbulent flows, and anomalous materials. Journal of Peridynamics and Nonlocal Modeling, pages 1–68, 2022. [117] JL Suzuki, M Zayernouri, ML Bittencourt, and GE Karniadakis. Fractional-order uni- axial visco-elasto-plastic models for structural analysis. Computer Methods in Applied Mechanics and Engineering, 308:443–467, 2016. [118] Giuseppe Failla and Massimiliano Zingales. Advanced materials modelling via frac- tional calculus: challenges and perspectives, 2020. [119] Jorge Suzuki, Yongtao Zhou, Marta D’Elia, and Mohsen Zayernouri. A thermodynam- ically consistent fractional visco-elasto-plastic model with memory-dependent damage for anomalous materials. Computer Methods in Applied Mechanics and Engineering, 373:113494, 2021. [120] Jorge L Suzuki, Maryam Naghibolhosseini, and Mohsen Zayernouri. A General Return-Mapping Framework for Fractional Visco-Elasto-Plasticity. arXiv preprint arXiv:2210.01308, 2022. [121] Jorge L Suzuki, Ehsan Kharazmi, Pegah Varghaei, Maryam Naghibolhosseini, and Mohsen Zayernouri. Anomalous nonlinear dynamics behavior of fractional viscoelastic beams. Journal of Computational and Nonlinear Dynamics, 16(11):111005, 2021. [122] Maryam Naghibolhosseini. Estimation of outer-middle ear transmission using DPOAEs and fractional-order modeling of human middle ear. PhD thesis, City Uni- versity of New York, 2015. [123] Maryam Naghibolhosseini and Glenis R Long. Fractional-order modelling and simu- lation of human ear. International Journal of Computer Mathematics, 95(6-7):1257– 1273, 2018. [124] Jorge L Suzuki, Tyler G Tuttle, Sara Roccabianca, and Mohsen Zayernouri. A data- driven memory-dependent modeling framework for anomalous rheology: Application to urinary bladder tissue. Fractal and Fractional, 5(4):223, 2021. [125] Patricio Clark Di Leoni, Tamer A Zaki, George Karniadakis, and Charles Meneveau. Two-point stress–strain-rate correlation structure and non-local eddy viscosity in tur- bulent flows. Journal of Fluid Mechanics, 914, 2021. 127 [126] Ali Akhavan-Safaei, Mehdi Samiee, and Mohsen Zayernouri. Data-driven fractional subgrid-scale modeling for scalar turbulence: A nonlocal LES approach. Journal of Computational Physics, 446:110571, 2021. [127] S Hadi Seyedi and Mohsen Zayernouri. A data-driven dynamic nonlocal subgrid-scale model for turbulent flows. Physics of Fluids, 34(3):035104, 2022. [128] S. Hadi Seyedi, Ali Akhavan-Safaei, and Mohsen Zayernouri. Dynamic nonlocal passive scalar subgrid-scale turbulence modeling. Physics of Fluids, 34(10):105122, 2022. [129] Brendan Keith, Ustim Khristenko, and Barbara Wohlmuth. A fractional PDE model for turbulent velocity fields near solid walls. Journal of Fluid Mechanics, 916, 2021. [130] Pavan Pranjivan Mehta, Guofei Pang, Fangying Song, and George Em Karniadakis. Discovering a universal variable-order fractional model for turbulent couette flow us- ing a physics-informed neural network. Fractional Calculus and Applied Analysis, 22(6):1675–1688, 2019. [131] Fangying Song and George Em Karniadakis. Variable-Order Fractional Models for Wall-Bounded Turbulent Flows. Entropy, 23(6):782, 2021. [132] Zellman Warhaft. Passive scalars in turbulent flows. Annual Review of Fluid Mechan- ics, 32(1):203–240, 2000. [133] Boris I Shraiman and Eric D Siggia. Scalar turbulence. Nature, 405(6787):639–646, 2000. [134] Stanley Corrsin. Further Generalization of Onsager’s Cascade Model for Turbulent Spectra. The Physics of Fluids, 7(8):1156–1159, 1964. [135] Pierre Sagaut and Claude Cambon. Homogeneous Turbulence Dynamics. Springer, 2018. [136] Parviz Moin and Krishnan Mahesh. Direct Numerical Simulation: A tool in turbulence research. Annual Review of Fluid Mechanics, 30(1):539–578, 1998. [137] Chris D Cantwell, David Moxey, Andrew Comerford, Alessandro Bolis, Gabriele Rocco, Gianmarco Mengaldo, Daniele De Grazia, Sergey Yakovlev, J-E Lombard, Dirk Ekelschot, et al. Nektar++: An open-source spectral/hp element framework. Computer physics communications, 192:205–219, 2015. [138] David Moxey, Chris D Cantwell, Yan Bao, Andrea Cassinelli, Giacomo Castiglioni, Sehun Chun, Emilia Juda, Ehsan Kazemi, Kilian Lackhove, Julian Marcon, et al. Nek- tar++: Enhancing the capability and application of high-fidelity spectral/hp element methods. Computer Physics Communications, 249:107110, 2020. [139] Ping He. A high order finite difference solver for massively parallel simulations of stably stratified turbulent channel flows. Computers & Fluids, 127:161–173, 2016. 128 [140] Paul Bartholomew, Georgios Deskos, Ricardo AS Frantz, Felipe N Schuch, Eric Lam- ballais, and Sylvain Laizet. Xcompact3D: An open-source framework for solving tur- bulence problems on a cartesian mesh. SoftwareX, 12:100550, 2020. [141] Paul T Bauman and Roy H Stogner. GRINS: a multiphysics framework based on the libmesh finite element library. SIAM Journal on Scientific Computing, 38(5):S78–S100, 2016. [142] Mikael Mortensen and Hans Petter Langtangen. High performance python for direct numerical simulations of turbulent flows. Computer Physics Communications, 203:53 – 65, 2016. [143] The OpenFOAM Foundation. Openfoam v8 user guide. [144] R.S. Rogallo. Numerical Experiments in Homogeneous Turbulence. NASA technical memorandum. National Aeronautics and Space Administration, 1981. [145] A. G. Lamorgese, D. A. Caughey, and S. B. Pope. Direct numerical simulation of homogeneous turbulence with hyperviscosity. Physics of Fluids, 17(1):015106, 2005. [146] Lisandro Dalcı́n, Rodrigo Paz, and Mario Storti. MPI for python. Journal of Parallel and Distributed Computing, 65(9):1108–1115, 2005. [147] Lisandro Dalcı́n, Rodrigo Paz, Mario Storti, and Jorge D’Elı́a. MPI for python: Per- formance improvements and MPI-2 extensions. Journal of Parallel and Distributed Computing, 68(5):655–662, 2008. [148] Lisandro D Dalcin, Rodrigo R Paz, Pablo A Kler, and Alejandro Cosimo. Parallel distributed computing using python. Advances in Water Resources, 34(9):1124–1139, 2011. [149] G. S. Patterson and Steven A. Orszag. Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions. The Physics of Fluids, 14(11):2538–2541, 1971. [150] Neal P. Sullivan, Shankar Mahalingam, and Robert M. Kerr. Deterministic forcing of homogeneous, isotropic turbulence. Physics of Fluids, 6(4):1612–1614, 1994. [151] V. Eswaran and S.B. Pope. An examination of forcing in direct numerical simulations of turbulence. Computers & Fluids, 16(3):257 – 278, 1988. [152] K. Alvelius. Random forcing of three-dimensional homogeneous turbulence. Physics of Fluids, 11(7):1880–1889, 1999. [153] A.S. Monin and A.M. Yaglom. Statistical Fluid Mechanics, Volume II: Mechanics of Turbulence. Dover Books on Physics. Dover Publications, 2013. [154] M. R. Overholt and S. B. Pope. Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence. Physics of Fluids, 8(11):3128–3148, 1996. 129 [155] George K Batchelor. Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity. Journal of Fluid Mechanics, 5(1):113–133, 1959. [156] Zhuo Wang, Kun Luo, Dong Li, Junhua Tan, and Jianren Fan. Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation. Physics of Fluids, 30(12):125101, 2018. [157] G. D. Portwood, B. T. Nadiga, J. A. Saenz, and D. Livescu. Interpreting neural network models of residual scalar flux. Journal of Fluid Mechanics, 907:A23, 2021. [158] Arvind T. Mohan, Dima Tretiak, Misha Chertkov, and Daniel Livescu. Spatio-temporal deep learning models of 3D turbulence with physics informed diagnostics. Journal of Turbulence, 0(0):1–41, 2020. [159] Arvind T Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov. Embed- ding hard physical constraints in neural network coarse-graining of 3D turbulence. arXiv preprint arXiv:2002.00021, 2020. [160] M Zayernouri and M Metzger. Coherent features in the sensitivity field of a planar mixing layer. Physics of Fluids, 23(2):025105, 2011. [161] Ali Akhavan-Safaei, S Hadi Seyedi, and Mohsen Zayernouri. Anomalous features in internal cylinder flow instabilities subject to uncertain rotational effects. Physics of Fluids, 32(9):094107, 2020. [162] M Gonzalez. Kinematic properties of passive scalar gradient predicted by a stochastic lagrangian model. Physics of Fluids, 21(5):055104, 2009. [163] Charles Meneveau. Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annual Review of Fluid Mechanics, 43:219–245, 2011. [164] Mohammad Farazmand and Themistoklis P Sapsis. A variational approach to probing extreme events in turbulent dynamical systems. Science advances, 3(9):e1701533, 2017. [165] Joseph Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiment. Monthly weather review, 91(3):99–164, 1963. [166] Jorge Bardina, J Ferziger, and W Reynolds. Improved subgrid-scale models for large- eddy simulation. In 13th fluid and plasmadynamics conference, page 1357. [167] Yan Zang, Robert L Street, and Jeffrey R Koseff. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows. Physics of Fluids A: Fluid Dynamics, 5(12):3186–3196, 1993. [168] Shewen Liu, Charles Meneveau, and Joseph Katz. On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet. Journal of Fluid Mechanics, 275:83–119, 1994. 130 [169] Andrea Beck and Marius Kurz. A Perspective on Machine Learning Methods in Tur- bulence Modelling. arXiv preprint arXiv:2010.12226, 2020. [170] Marius Kurz and Andrea Beck. A machine learning framework for LES closure terms. arXiv preprint arXiv:2010.03030, 2020. [171] Justin Sirignano, Jonathan F. MacArt, and Jonathan B. Freund. DPM: A deep learning PDE augmentation method with application to large-eddy simulation. Journal of Computational Physics, page 109811, 2020. [172] Stewart Harris. An introduction to the theory of the Boltzmann equation. Courier Corporation, 2004. [173] P. L. Bhatnagar, E. P. Gross, and M. Krook. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev., 94:511–525, May 1954. [174] Yoshio Sone. Kinetic theory and fluid dynamics. Springer Science & Business Media, 2012. [175] K. Huang. Statistical Mechanics. Wiley, 1987. [176] Keerti Vardhan Sharma, Robert Straka, and Frederico Wanderley Tavares. Current status of Lattice Boltzmann Methods applied to aerodynamic, aeroacoustic, and ther- mal flows. Progress in Aerospace Sciences, 115:100616, 2020. [177] A Bartoloni, C Battista, S Cabasino, PS Paolucci, J Pech, R Sarno, GM Todesco, M Torelli, W Tross, P Vicini, et al. LBE simulations of Rayleigh-Benard convection on the APE100 parallel processor. International Journal of Modern Physics C, 4(05):993– 1006, 1993. [178] JGM Eggels and JA Somers. Numerical simulation of free convective flow using the lattice-Boltzmann scheme. International Journal of Heat and Fluid Flow, 16(5):357– 364, 1995. [179] Xiaowen Shan. Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method. Physical review E, 55(3):2780, 1997. [180] Hudong Chen, Steven A Orszag, and Ilya Staroselsky. Macroscopic description of arbitrary Knudsen number flow using Boltzmann-BGK kinetic theory. Journal of Fluid Mechanics, 574:495, 2007. [181] Hudong Chen, Steven A Orszag, and Ilya Staroselsky. Macroscopic description of arbitrary Knudsen number flow using Boltzmann-BGK kinetic theory. Part 2. Journal of Fluid Mechanics, 658:294, 2010. [182] Sharath S Girimaji. Boltzmann kinetic equation for filtered fluid turbulence. Physical review letters, 99(3):034501, 2007. 131 [183] Pierre Sagaut. Toward advanced subgrid models for lattice-boltzmann-based large- eddy simulation: theoretical formulations. Computers & Mathematics with Applica- tions, 59(7):2194–2199, 2010. [184] Jesse Chu-Shore, M Brandon Westover, and Matt T Bianchi. Power law versus expo- nential state transition dynamics: application to sleep-wake architecture. PLoS One, 5(12):e14204, 2010. [185] David Applebaum. Lévy processes and stochastic calculus. Cambridge university press, 2009. [186] Laure Saint-Raymond. Hydrodynamic limits of the Boltzmann equation. Number 1971. Springer Science & Business Media, 2009. [187] S. Beetham and J. Capecelatro. Formulating turbulence closures using sparse regres- sion with embedded form invariance. Phys. Rev. Fluids, 5:084611, Aug 2020. [188] Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. Journal of the royal statistical society: series B (statistical methodology), 67(2):301– 320, 2005. [189] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blon- del, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011. [190] Charles Meneveau and Joseph Katz. Scale-invariance and turbulence models for large- eddy simulation. Annual Review of Fluid Mechanics, 32(1):1–32, 2000. [191] Boris I. Shraiman and Eric D. Siggia. Lagrangian path integrals and fluctuations in random flow. Phys. Rev. E, 49:2912–2927, Apr 1994. [192] J. P. Gollub, J. Clarke, M. Gharib, B. Lane, and O. N. Mesquita. Fluctuations and transport in a stirred fluid with a mean gradient. Phys. Rev. Lett., 67:3507–3510, Dec 1991. [193] Jayesh and Z. Warhaft. Probability distribution of a passive scalar in grid-generated turbulence. Phys. Rev. Lett., 67:3503–3506, Dec 1991. [194] Jayesh and Z. Warhaft. Probability distribution, conditional dissipation, and transport of passive temperature fluctuations in grid-generated turbulence. Physics of Fluids A: Fluid Dynamics, 4(10):2292–2307, 1992. [195] B. R. Lane, O. N. Mesquita, S. R. Meyers, and J. P. Gollub. Probability distributions and thermal transport in a turbulent grid flow. Physics of Fluids A: Fluid Dynamics, 5(9):2255–2263, 1993. [196] Reginald J Hill. Models of the scalar spectrum for turbulent advection. Journal of Fluid Mechanics, 88(3):541–562, 1978. 132 [197] Katepalli R Sreenivasan. The passive scalar spectrum and the Obukhov–Corrsin con- stant. Physics of Fluids, 8(1):189–196, 1996. [198] DA Donzis and PK Yeung. Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence. Physica D: Nonlinear Phenomena, 239(14):1278– 1287, 2010. [199] KP Iyer and PK Yeung. Structure functions and applicability of Yaglom’s relation in passive-scalar turbulent mixing at low Schmidt numbers with uniform mean gradient. Physics of Fluids, 26(8):085107, 2014. [200] Antoine Vollant, Guillaume Balarac, and C Corre. Subgrid-scale scalar flux modelling based on optimal estimation theory and machine-learning procedures. Journal of Tur- bulence, 18(9):854–878, 2017. [201] Adam Subel, Ashesh Chattopadhyay, Yifei Guan, and Pedram Hassanzadeh. Data- driven subgrid-scale modeling of forced burgers turbulence using deep learning with generalization to higher reynolds numbers via transfer learning. Physics of Fluids, 33(3):031702, 2021. [202] Yifei Guan, Ashesh Chattopadhyay, Adam Subel, and Pedram Hassanzadeh. Stable a posteriori LES of 2D turbulence using convolutional neural networks: Backscattering analysis and generalization to higher Re via transfer learning. Journal of Computa- tional Physics, 458:111090, 2022. [203] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson. How transferable are features in deep neural networks? Advances in neural information processing systems, 27, 2014. [204] MR Overholt and SB Pope. Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence. Physics of Fluids, 8(11):3128–3148, 1996. [205] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. 2017. [206] D. K. Lilly. The representation of small-scale turbulence in numerical simulation ex- periments. in Proceedings of the IBM Computer Science Symposium on Environmental Sciences, page p. 167, Yorktown Heights, NY, 1967. [207] M Antonopoulos-Domis. Large-eddy simulation of a passive scalar in isotropic turbu- lence. Journal of Fluid Mechanics, 104:55–79, 1981. [208] Robert A Clark, Joel H Ferziger, and William Craig Reynolds. Evaluation of subgrid- scale models using an accurately simulated turbulent flow. Journal of fluid mechanics, 91(1):1–16, 1979. 133 [209] KW Bedford and WK Yeo. Conjunctive filtering procedures in surface water flow and transport. Large eddy simulation of complex engineering and geophysical flows, pages 513–537, 1993. [210] Bert Vreman, Bernard Geurts, and Hans Kuerten. Large-eddy simulation of the tur- bulent mixing layer. Journal of fluid mechanics, 339:357–390, 1997. [211] Hao Lu and Fernando Porté-Agel. A modulated gradient model for scalar trans- port in large-eddy simulation of the atmospheric boundary layer. Physics of Fluids, 25(1):015110, 2013. [212] Antoine Vollant, Guillaume Balarac, and Christophe Corre. A dynamic regularized gradient model of the subgrid-scale stress tensor for large-eddy simulation. Physics of Fluids, 28(2):025114, 2016. [213] Armen Der Kiureghian and Ove Ditlevsen. Aleatory or epistemic? does it matter? Structural Safety, 31(2):105 – 112, 2009. Risk Acceptance and Risk Communication. [214] Daniel M. Tartakovsky and Dongbin Xiu. Stochastic analysis of transport in tubes with rough walls. Journal of Computational Physics, 217(1):248 – 259, 2006. Uncertainty Quantification in Simulation Science. [215] Mohsen Zayernouri, Sang-Woo Park, Daniel M. Tartakovsky, and George Em Karni- adakis. Stochastic smoothed profile method for modeling random roughness in flow problems. Computer Methods in Applied Mechanics and Engineering, 263:99 – 112, 2013. [216] Stefano Berrone, Claudio Canuto, Sandra Pieraccini, and Stefano Scialò. Uncertainty quantification in discrete fracture network models: Stochastic geometry. Water Re- sources Research, 54(2):1338–1352, 2018. [217] Chunsong Kwon and Daniel M. Tartakovsky. Modified immersed boundary method for flows over randomly rough surfaces. Journal of Computational Physics, 406:109195, 2020. [218] N. Vu-Bac, M. Silani, T. Lahmer, X. Zhuang, and T. Rabczuk. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Com- putational Materials Science, 96:520 – 535, 2015. Special Issue Polymeric Composites. [219] E. E. Prudencio, P. T. Bauman, D. Faghihi, K. Ravi-Chandar, and J. T. Oden. A computational framework for dynamic data-driven material damage control, based on Bayesian inference and model selection. International Journal for Numerical Methods in Engineering, 102(3-4):379–403, 2015. [220] Kathryn Farrell, J. Tinsley Oden, and Danial Faghihi. A Bayesian framework for adaptive selection, calibration, and validation of coarse-grained models of atomistic systems. Journal of Computational Physics, 295:189 – 208, 2015. 134 [221] Rebecca E. Morrison, Todd A. Oliver, and Robert D. Moser. Representing model inadequacy: A stochastic operator approach. SIAM/ASA Journal on Uncertainty Quantification, 6(2):457–496, 2018. [222] S Hadi Seyedi, Behzad Nemati Saray, and Ali Ramazani. High-accuracy multiscale simulation of three-dimensional squeezing carbon nanotube-based flow inside a rotating stretching channel. Mathematical Problems in Engineering, 2019, 2019. [223] Guangping Li, Jalil Manafian, Mirmehdi Seyyedi, R Sivaraman, and Subhiya M Zey- nalli. Periodic, Cross-Kink, and Interaction between Stripe and Periodic Wave Solu- tions for Generalized Hietarinta Equation: Prospects for Applications in Environmen- tal Engineering. Advances in Mathematical Physics, 2022, 2022. [224] Eduardo A Barros de Moraes, Mohsen Zayernouri, and Mark M Meerschaert. An integrated sensitivity-uncertainty quantification framework for stochastic phase-field modeling of material damage. International Journal for Numerical Methods in Engi- neering, 122(5):1352–1377, 2021. [225] Eduardo A Barros de Moraes, Hadi Salehi, and Mohsen Zayernouri. Data-driven failure prediction in brittle materials: A phase field-based machine learning framework. Journal of Machine Learning for Modeling and Computing, 2(1), 2021. [226] Eduardo A Barros de Moraes, Jorge L Suzuki, and Mohsen Zayernouri. Atomistic-to- meso multi-scale data-driven graph surrogate modeling of dislocation glide. Computa- tional Materials Science, 197:110569, 2021. [227] Eduardo Augusto Barros de Moraes, Marta D’Elia, and Mohsen Zayernouri. Nonlocal machine learning of micro-structural defect evolutions in crystalline materials. arXiv preprint arXiv:2205.05729, 2022. [228] Ali A. Safaei. Analysis and implementation of multiple models and multi-models for shallow-water type models of large mass flows. Master’s thesis, State University of New York at Buffalo, 2018. [229] Abani Patra, Andrea Bevilacqua, Ali Akhavan-Safaei, E. Bruce Pitman, Marcus Bur- sik, and David Hyman. Comparative analysis of the structures and outcomes of geo- physical flow models and modeling assumptions using uncertainty quantification. Fron- tiers in Earth Science, 8:275, 2020. [230] Abani K Patra, Andrea Bevilacqua, and Ali Akhavan Safaei. Analyzing complex mod- els using data and statistics. In International Conference on Computational Science, pages 724–736. Springer, 2018. [231] Dongbin Xiu and Spencer J Sherwin. Parametric uncertainty analysis of pulse wave propagation in a model of a human arterial network. Journal of Computational Physics, 226(2):1385–1407, 2007. 135 [232] C Gorlé and G Iaccarino. A framework for epistemic uncertainty quantification of turbulent scalar flux models for reynolds-averaged navier-stokes simulations. Physics of Fluids, 25(5):055105, 2013. [233] Mohammad Khalil, Guilhem Lacaze, Joseph C. Oefelein, and Habib N. Najm. Uncer- tainty quantification in LES of a turbulent bluff-body stabilized flame. Proceedings of the Combustion Institute, 35(2):1147 – 1156, 2015. [234] Catherine Gorlé, Stéphanie Zeoli, Michael Emory, J Larsson, and G Iaccarino. Epis- temic uncertainty quantification for reynolds-averaged navier-stokes modeling of sepa- rated flows over streamlined surfaces. Physics of Fluids, 31(3):035101, 2019. [235] Andrea F. Cortesi, Paul G. Constantine, Thierry E. Magin, and Pietro M. Congedo. Forward and backward uncertainty quantification with active subspaces: Application to hypersonic flows around a cylinder. Journal of Computational Physics, 407:109079, 2020. [236] Zengrong Hao and Catherine Gorlé. Quantifying turbulence model uncertainty in reynolds-averaged navier-stokes simulations of a pin-fin array. part 1: flow field. Com- puters & Fluids, page 104641, 2020. [237] Zengrong Hao and Catherine Gorlé. Quantifying turbulence model uncertainty in reynolds-averaged navier-stokes simulations of a pin-fin array. part 2: scalar transport. Computers & Fluids, page 104642, 2020. [238] Dongbin Xiu and George Em Karniadakis. Modeling uncertainty in flow simulations via generalized polynomial chaos. Journal of computational physics, 187(1):137–167, 2003. [239] Mohammad Khalil and Abhijit Sarkar. Independent component analysis to enhance performances of karhunen–loeve expansions for non-gaussian stochastic processes: Ap- plication to uncertain systems. Journal of Sound and Vibration, 333(21):5600–5613, 2014. [240] Bernt Oksendal. Stochastic differential equations: an introduction with applications. Springer Science & Business Media, 2013. [241] George Karniadakis and Spencer Sherwin. Spectral/hp element methods for computa- tional fluid dynamics. Oxford University Press, 2013. [242] C.D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. De Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, B. Jordi, H. Xu, Y. Mohamied, C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R.M. Kirby, and S.J. Sherwin. Nektar++: An open-source spectral/hp element framework. Computer Physics Communications, 192:205 – 219, 2015. [243] David Moxey, Chris D. Cantwell, Yan Bao, Andrea Cassinelli, Giacomo Castiglioni, Sehun Chun, Emilia Juda, Ehsan Kazemi, Kilian Lackhove, Julian Marcon, Gianmarco Mengaldo, Douglas Serson, Michael Turner, Hui Xu, Joaquim Peiró, Robert M. Kirby, 136 and Spencer J. Sherwin. Nektar++: Enhancing the capability and application of high- fidelity spectral/hp element methods. Computer Physics Communications, 249:107110, 2020. [244] Robert M. Kirby and Spencer J. Sherwin. Stabilisation of spectral/hp element meth- ods through spectral vanishing viscosity: Application to fluid mechanics modelling. Computer Methods in Applied Mechanics and Engineering, 195(23):3128 – 3144, 2006. Incompressible CFD. [245] Christophe Geuzaine and Jean-François Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11):1309–1331, 2009. [246] Dongbin Xiu and Jan S Hesthaven. High-order collocation methods for differential equations with random inputs. SIAM Journal on Scientific Computing, 27(3):1118– 1139, 2005. [247] Fabio Nobile, Raúl Tempone, and Clayton G Webster. A sparse grid stochastic colloca- tion method for partial differential equations with random input data. SIAM Journal on Numerical Analysis, 46(5):2309–2345, 2008. [248] Jasmine Foo, Xiaoliang Wan, and George Em Karniadakis. The multi-element prob- abilistic collocation method (ME-PCM): Error analysis and applications. Journal of Computational Physics, 227(22):9572–9595, 2008. [249] Andreas Glaser, Xiangtao Liu, and Vladimir Rokhlin. A fast algorithm for the cal- culation of the roots of special functions. SIAM Journal on Scientific Computing, 29(4):1420–1438, 2007. [250] I.M Sobol’. Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. Mathematics and Computers in Simulation, 55(1):271 – 280, 2001. The Second IMACS Seminar on Monte Carlo Methods. [251] Andrea Saltelli, Paola Annoni, Ivano Azzini, Francesca Campolongo, Marco Ratto, and Stefano Tarantola. Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index. Computer Physics Communications, 181(2):259 – 270, 2010. [252] I.M Sobol’. On sensitivity estimation for nonlinear mathematical models. Matem. Mod., 2(1):112 – 118, 1990. [253] Olivier Le Maı̂tre and Omar M Knio. Spectral methods for uncertainty quantification: with applications to computational fluid dynamics, 2010. [254] Hyung Jun Yang, Francesca Boso, Hamdi A Tchelepi, and Daniel M Tartakovsky. Method of distributions for quantification of geologic uncertainty in flow simulations. Water Resources Research, page e2020WR027643. 137 [255] Khader M. Hamdia, Mohammed A. Msekh, Mohammad Silani, Nam Vu-Bac, Xiaoy- ing Zhuang, Trung Nguyen-Thoi, and Timon Rabczuk. Uncertainty quantification of the fracture properties of polymeric nanocomposites based on phase field modeling. Composite Structures, 133:1177 – 1190, 2015. [256] D. Faghihi, S. Sarkar, M. Naderi, J.E. Rankin, L. Hackel, and N. Iyyer. A probabilistic design method for fatigue life of metallic component. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part B: Mechanical Engineering, 4(3), 2018. [257] Rodolfo Ostilla-Mónico, Erwin P. van der Poel, Roberto Verzicco, Siegfried Grossmann, and Detlef Lohse. Boundary layer dynamics at the transition between the classical and the ultimate regime of Taylor-Couette flow. Physics of Fluids, 26(1):015114, 2014. [258] Siegfried Grossmann, Detlef Lohse, and Chao Sun. Velocity profiles in strongly turbu- lent Taylor-Couette flow. Physics of Fluids, 26(2):025114, 2014. [259] Dongbin Xiu and Daniel M Tartakovsky. Numerical methods for differential equations in random domains. SIAM Journal on Scientific Computing, 28(3):1167–1185, 2006. [260] George K Batchelor. Computation of the energy spectrum in homogeneous two- dimensional turbulence. The Physics of Fluids, 12(12):II–233, 1969. [261] R Benzi, S Patarnello, and P Santangelo. Self-similar coherent structures in two- dimensional decaying turbulence. Journal of Physics A: Mathematical and General, 21(5):1221, 1988. [262] M Fathali and M Khoshnami Deshiri. Sensitivity of the two-dimensional shearless mixing layer to the initial turbulent kinetic energy and integral length scale. Physical Review E, 93(4):043122, 2016. [263] Siwei Dong, Yongxiang Huang, Xianxu Yuan, and Adrián Lozano-Durán. The coherent structure of the kinetic energy transfer in shear turbulence. Journal of Fluid Mechanics, 892, 2020. [264] Xiang IA Yang, Sergio Pirozzoli, and Mahdi Abkar. Scaling of velocity fluctuations in statistically unstable boundary-layer flows. Journal of Fluid Mechanics, 886, 2020. [265] Marco Atzori, Ricardo Vinuesa, Adrián Lozano-Durán, and Philipp Schlatter. Co- herent structures in turbulent boundary layers over an airfoil. In Journal of Physics: Conference Series, volume 1522, page 012020. IOP Publishing, 2020. [266] Amir M. Akbarzadeh and Iman Borazjani. Large eddy simulations of a turbulent channel flow with a deforming wall undergoing high steepness traveling waves. Physics of Fluids, 31(12):125107, 2019. [267] AM Akbarzadeh and I Borazjani. Controlling flow separation on a thick airfoil using backward traveling waves. AIAA Journal, pages 1–8, 2020. 138 [268] Anna Lischke, Guofei Pang, Mamikon Gulian, Fangying Song, Christian Glusa, Xi- aoning Zheng, Zhiping Mao, Wei Cai, Mark M Meerschaert, Mark Ainsworth, et al. What is the fractional Laplacian? A comparative review with new results. Journal of Computational Physics, 404:109009, 2020. [269] Elias M Stein. Singular integrals and differentiability properties of functions (PMS-30), volume 30. Princeton university press, 2016. [270] GP Neitzel and Stephen H Davis. Energy stability theory of decelerating swirl flows. The Physics of Fluids, 23(3):432–437, 1980. [271] Scott Callaghan, Gideon Juve, Karan Vahi, Philip J Maechling, Thomas H Jordan, and Ewa Deelman. rvGAHP: push-based job submission using reverse SSH connections. In Proceedings of the 12th Workshop on Workflows in Support of Large-Scale Science, pages 1–8, 2017. [272] Nikolay A Simakov, Renette L Jones-Ivey, Ali Akhavan-Safaei, Hossein Aghakhani, Matthew D Jones, and Abani K Patra. Modernizing Titan2D, a parallel AMR geo- physical flow code to support multiple rheologies and extendability. In International Conference on High Performance Computing, pages 101–112. Springer, 2019. 139 APPENDIX A FRACTIONAL-ORDER DIFFERENTIAL OPERATORS According to Lischke et al. [268], the fractional Laplacian operator, denoted by (−∆)α with 0 < α ≤ 1, is defined as Z (−∆)α u(x) = 1 |ξ|2α u, e−iξ·x  eiξ·x dξ (2π)d Rd n  o = F −1 2α |ξ| F u (ξ) , (A.1) where F and F −1 represent the Fourier and inverse Fourier transforms for a real-valued √ vector ξ = ξj , j = 1, 2, 3, respectively, and i = −1. Moreover, (· , ·) specifies the L2 -inner product on Rd , d = 1, 2, 3. Therefore, the Fourier transform of the fractional Laplacian is then obtained as n o  F (−∆)α u(x) = |ξ|2α F u (ξ), (A.2) where α = 1 recovers the integer-order Laplacian. Considering the definition of α-Riesz potential as Z u(x) − u(s) Iα u(x) = Cd,−α d−2α ds, (A.3) Rd |x − s| the fractional Laplacian can also be expressed in the integral form as Z u(x) − u(s) (−∆)α u(x) = Cd,α 2α+d ds, (A.4) Rd |x − s| 22α Γ(α+d/2) where Cd,α = d/2 for α ∈ (0, 1] and Γ(·) represents Gamma function [268]. The π Γ(−α) α-Riesz potential is also formulated [269] as n  o Iα u(x) = (−∆)−α u(x) = F −1 |ξ| −2α F u (ξ) . (A.5) Considering (A.5), the Riesz transform is then given by n i ξj  o Rj u(x) = ∇j I1 u(x) = F −1 − F u (ξ) , (A.6) |ξ| which is utilized in formulating the SGS scalar flux. 140 APPENDIX B DERIVATION OF PASSIVE SCALAR FLUX Regarding the filtered passive scalar flux given in (3.36), one can write that Z ∞Z   qei = (vi − u ei ) Φ g s,s F ( g Ls,s ) − e Φ F ( e L) e−s dv ds, (B.1) 0 Rd where using the Taylor expansions of Φ gs,s and F (L g s,s ) about their not shifted values and later on by utilizing the incompressibility constraint, one arrives at the following Z ∞Z e  ∂Φ qei = (vi − u ei ) −vj sτg F (L)e e−s dv ds (B.2) 0 R d ∂x j e  Z ∞  Z τg ∂Φ ≃ − 3 s e−s ds (vj − u ej )(vj − u e dv. ej ) F (L) cT ∂xi 0 Rd R Knowing that 0∞ s e−s ds = 1, the diffusivity coefficient of the passive scalar, D, would be expressed as Z τg e dv. D := 3 (vj − u ej )(vj − u ej ) F (L) (B.3) cT Rd As a result, the filtered (resolved) passive scalar flux, q e, and its divergence appearing in the right-hand side of (3.31) could be written as ∂Φe qei = −D =⇒ ∇·q e = −D ∆ Φ. e (B.4) ∂xi On the other hand, the integral form of the modeled SGS flux in (3.37) can be written in the following form Z ∞Z   qiR = (vi − u ei ) Φ gs,s F α (Lg s,s ) − e Φ F α (L) e e−s dv ds, (B.5) 0 Rd   By adding and subtracting Φ gs,s F e α (L) to Φg s,s F α (Lg s,s ) − e Φ F e α (L) , one can rewrite (B.5) as Z ∞Z   qiR = (vi − u ei ) Φ g s,s − e Φ F α (L) e e−s dv ds (B.6) 0Z Rd ∞Z   + (vi − u g ei ) Φ s,s F α (Lg s,s ) − F α (L) e e−s dv ds, 0 Rd 141 where the second integral is approximated as zero. Therefore, the modeled subgrid-scale scalar flux is simplified into Z ∞Z   qiR = (vi − u ei ) Φg s,s − e Φ F α (L) e e−s dv ds. (B.7) 0 Rd Considering vi = (x′i − xi )/sτ and approximating vi − uei ≃ vi , one can obtain that dv = dx′ /(sτg )3 [16]. As a result, Le = (v − u e)2 /c2T ≈ v 2 /c2T = (x′ − x)2 /(s τg cT )2 . According to the definition of isotropic α-stable Lévy distribution, F α (L) e = Cα /Le(2α+3)/2 , where Cα is a real-valued constant. Consequently, (B.7) may be reformulated as Z ∞Z  ′    1 xi − x i  g e  Cα qiR = 3 Φs,s − Φ e−s dx′ ds (B.8) 0 3 R cT s τ g d 3 s τg e L (2α+3)/2   2α Z ∞  Z (x′ − x ) Φ(x e e ′ ) − Φ(x) Cα (cT τg ) i i = − e−s s2α−1 ds ′ 2α+3 dx′ τg 0 d |x − x| R  Z ′ e (xi − xi ) Φ(x ) − Φ(x)′ e Cα (cT τg )2α = − Γ(2α) ′ − x|2α+3 dx′ , τg R d |x By taking the divergence of the modeled SGS scalar flux in (B.8), we obtain   2α Z (x′ − x ) Φ(x e e ′ ) − Φ(x) Cα (cT τg ) i i (∇ · q R )i = − Γ(2α) ∇ · ′ 2α+3 dx′ (B.9) τg R d |x − x| " Z e ′ Z # Cα (cT τg )2α Φ(x ) − Φ(x) ′ e e ∂ Φ/∂xi = − Γ(2α) (2α + 2) ′ 2α+3 dx − ′ 2α+2 dx′ . τg Rd |x − x| Rd |x − x| Due to symmetry, the second integral inside the large brace is zero and the other integral is nothing but the definition of fractional Laplacian of the filtered passive scalar concentration, e Thus, (B.9) takes the following compact form Φ. Cα (cT τg )2α ∇ · qR = − (2α + 2) Γ(2α) (−∆)α Φ, e α ∈ (0, 1]. (B.10) 2τg 142 APPENDIX C VALIDATION OF NUMERICAL SETUP This appendix provides a comparison study between the analytic and numerical solutions for specific cases of impulsive and exponential spin-decay at low-Reynolds numbers in order to validate our CFD results. Simplifying the governing equations in cylindrical coordinate system, (r, θ, z), for a non-stationary 2-D viscous incompressible flow, gives ! u2θ ∂p ρ − =− , (C.1) r ∂r    2  ∂uθ ∂ uθ 1 1 ∂uθ ρ =µ − u + . ∂t ∂r2 r2 θ r ∂r Here, the first and second equations represent the momentum equation in r and θ directions, respectively. By considering no-slip boundary conditions on the wall and taking the initial condition as V (r, 0) = rθ̇ (rigid-body rotation), equation (C.1) can be solved through the Laplace transform on the variable t [270, 46]. If the length is scaled by the radius of cylinder, r, time is scaled by r2 /ν, velocity in the sudden stop case is scaled by rθ̇, and velocity in the exponential decay case by λr3 θ̇/ν, the resulting solution would be dimensionless. Therefore, the exact solutions for the complete sudden stop and exponential decay cases at low-Reynolds numbers are obtained as X∞ J1 (βn r) Vs (r, t) = −2 exp(−βn2 t), (C.2) βn J0 (βn ) n=1 √ ∞ R− J1 (r B) exp(−Bt) X J1 (βn r) exp(−βn2 t) Ve (r, t) = √ +2 , Re J1 ( B) βn (β 2 − β)J0 (βn ) n=1 where Vs (r, t) indicates the azimuth velocity for sudden stop case, Ve (r, t) is the azimuth velocity for the exponential decay case, J is the Bessel function of the first kind, and βn denotes the positive roots of J1 (βn ) = 0. Also R− = r2 θ̇/ν shows the Reynolds number corresponding to the initial state and Re = r4 θ̇λ/ν 2 denotes the Reynolds number for the spin-decay period (see [270] and [46] for derivations). Using equation (C.2) and implementing the same initial and boundary conditions in the numerical setup for a low Re number, a 143 comparison in different times was made (see Figure C.1). These comparisons are obtained for Re = 1/ν = 100 and Re /R− = 20, while we consider the mentioned dimensionless solution and the physical parameters. Comparing the analytic and the CFD results clearly validates our numerical implementation and procedure. It should be mentioned that the analytic solutions are only valid at the low-Re number regime where no flow instability is created during these processes. Analytic CFD (t = 10−3) Analytic CFD (t = 0.01) 0.8 Analytic 0.04 CFD (t = 10−2) Analytic CFD (t = 0.02) Analytic CFD (t = 10−1) Analytic CFD (t = 0.05) 0.6 0.03 V (r, t) 0.4 V (r, t) 0.02 0.2 0.01 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r r (a) (b) Figure C.1: Comparison between the velocity, V (r, t), obtained from CFD and analytical solution for flow at Re = 100. (a) Complete sudden stop, (b) exponential decay with Re /R− = 20. 144 APPENDIX D COMPUTATIONAL WORKFLOW Performing numerous amount of forward simulations for discretization of random space urges the design of a proper workflow in high-performance computing (HPC) environment [271, 272]. In this work, we are dealing with a forward solver with requires input session files in the xml format, which contain information about the grid and each forward simulation’s conditions. Using parallel computing on O(100) processes is inevitably demanded for each one of these forward simulations. Indeed, the number of simulations addressed in this work, could not be achieved by manually generation of input session files that are fed by realiza- tions of stochastic parameter space. Hence, a Python program is prepared to construct the parameter space realizations (either from MC approach or PCM) and assign them to separate xml scripts that are placed in a directory associated with each forward simulation. Moreover, it enables automation of job submission step in the HPC environment. The statistical so- lutions (i.e., expected fields and their standard deviation) are computed by post-processing through Paraview toolkit. In particular, we exploit Paraview’s Python scripting (executed by pvpython) to extract the flow field variables from xml field files at the SEM integration points and perform required computations on them to obtain the expectation and standard deviation of field variables. Similar procedure is carried out to compute the velocity and vorticity fluctuation fields. 145