DATA-DRIVEN DYNAMIC NONLOCAL SUBGRID-SCALE MODELING FOR THE LARGE EDDY SIMULATION OF TURBULENT FLOWS By Seyedhadi Seyedi A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2022 ABSTRACT This study aims to propose novel solutions to the complex problem of turbulent flows using data-driven statistical and mathematical models. The proposed models reduce the huge computational cost of the direct numerical simulations and make them tractable while maintaining the important statistical features of the chaotic flows. Unlike the conventional models in the literature, the new proposed dynamic models take into account the inherent nonlocality of turbulence and predict the final statistical quantities with higher accuracy and correlations. First, we developed a novel autonomously dynamic nonlocal turbulence model for the large and very large eddy simulation (LES, VLES) of the homogeneous isotropic turbulent flows (HIT). The model is based on a generalized (integer-to-noninteger) order Laplacian of the filtered velocity field, and a novel dynamic model has been formulated to avoid the need for tuning the model constant. Three data-driven approaches were introduced for the determination of the fractional-order to have a model which is totally free of any tuning parameter. Our analysis includes both the a priori and the a posteriori tests. In the former test, using a high-fidelity and well-resolved dataset from direct numerical simulations (DNS), we computed the correlation coefficients for the stress components of the subgrid- scale (SGS) stress tensor and the one we get directly from the DNS results. Moreover, we compared the probability density function of the ensemble-averaged SGS forces for different filter sizes. In the latter, we employed our new model along with other conventional models including static and dynamic Smagorinsky into our pseudo-spectral solver and tested the final predicted quantities. The results of the newly developed model exhibit an expressive agreement with the ground-truth DNS results in all components of the SGS stress and forces. Also, the model exhibits promising results in the VLES region as well as the LES region, which could be remarkably important for the cost-efficient nonlocal turbulence modeling e.g., in meteorological and environmental applications. Afterwards, we extend the same dynamic nonlocal idea to the scalar turbulence. To this end, we formulate the underlying nonlocal model starting from the filtered Boltzmann kinetic transport equation, where the divergence of subgrid-scale scalar fluxes emerges as a fractional-order Laplacian term in the filtered advection-diffusion model, coding the corre- sponding super-diffusive nature of scalar turbulence. Subsequently, we develop a robust data- driven algorithm for estimation of the fractional (non-integer) Laplacian exponent, where we on-the-fly calculate the corresponding model coefficient employing a new dynamic procedure. Our a priori tests show that our new dynamically nonlocal LES paradigm provides better agreements with the ground-truth filtered DNS data in comparison to the conventional static and dynamic Prandtl-Smagorisnky models. Moreover, in order to analyze the numerical sta- bility and assessing the model’s performance, we carry out a comprehensive a posteriori tests. They unanimously illustrate that our new model considerably outperforms other ex- isting functional models, correctly predicting the backscattering phenomena at the same time and providing higher correlations at small-to-large filter sizes. We conclude that our proposed nonlocal subgrid-scale model for scalar turbulence is amenable for coarse LES and VLES frameworks even with strong anisotropies, applicable to environmental applications. Finally, we developed a new dynamic tempered fractional subgrid-scale model, DTF, for the large and very large eddy simulation of turbulent flows. The nonlocality of the turbulent flows is the innate feature that can be seen in the non-Gaussian statistics of the velocity increments and can be addressed properly by the nonlocal models in terms of the fractional operators. Using kinetic transport, we developed a dynamic tempered fractional model that encompasses the three main characteristics of an ideal turbulence model: (i) nonlocal nature, (ii) dynamic model constant computations, and (iii) tempered and finite variance property. Several simulations of forced homogeneous isotropic and multi-layer temporal shear layer turbulent flows have been done in the a priori and a posteriori analyses. The results show that the new model is not only numerically stable and can maintain low- and high-order structures in long-range simulations, but it also provides better predictions than local models and nontempered models. Copyright by SEYEDHADI SEYEDI 2022 ACKNOWLEDGEMENTS I would like to express my sincere gratitude to my advisor, Dr. Mohsen Zayernouri, for his continual unfailing support and guidance throughout my research. The consistent encouragement and kindness he showed me will never be forgotten. I would like to thank the respected committee members, Dr. Yuan, Dr. Gao, and Dr. Lei for all the research discussions and helpful feedback that have helped to improve the quality of this dissertation. I would also like to thank all current and past FMATH group members: Ali, Jorge, Eduardo, Pegah, Mehdi, Vikram, Ehsan and Yongtao. As always, I express my gratitude to my mom, dad, and brother for their endless and unconditional love and support, devotion, and encouragement. I would not have made it this far without them. Additionally, I thank my dear wife for her companionship, patience, and emotional support throughout these years. This work was supported by the department of mechanical engineering and was finan- cially funded by the MURI/ARO grant (W911NF-15-1-0562), the ARO Young Investiga- tor Program (YIP) award (W911NF-19-1-0444), the AFOSR Young Investigator Program (YIP) (FA9550-17-1-0150), and partially by the National Science Foundation award (DMS- 1923201). The high-performance computing resources and services were provided by the Institute for Cyber-Enabled Research (ICER) at Michigan State University. v TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline of the Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 CHAPTER 2 ANOMALOUS FEATURES IN INTERNAL CYLINDER FLOW . . . 8 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Stochastic Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Stochastic Computational Fluid Dynamics Framework . . . . . . . . . . . . 14 2.5 Variance-Based Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 CHAPTER 3 A DATA-DRIVEN DYNAMIC NONLOCAL SUBGRID-SCALE MODEL FOR TURBULENT FLOWS . . . . . . . . . . . . . . . . . . . . . . 36 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 A Priori Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4 A Posteriori Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 CHAPTER 4 DYNAMIC NONLOCAL CLOSURE MODELING FOR SCALAR TURBULENCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2 Development of the new model (DNPS) . . . . . . . . . . . . . . . . . . . . . 60 4.3 Data-driven Identification of Optimum Fractional Order, αopt . . . . . . . . 64 4.4 A Priori Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.5 A Posteriori Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 CHAPTER 5 A DYNAMIC TEMPERED NONLOCAL TURBULENCE MODEL FOR LARGE EDDY SIMULATION OF CHAOTIC FLOWS . . . . . 76 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.2 Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3 A Priori Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.4 A Posteriori Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 CHAPTER 6 SUMMARY AND FUTURE WORKS . . . . . . . . . . . . . . . . . 97 6.1 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 APPENDIX A VALIDATION OF NUMERICAL SETUP . . . . . . . . . . . . . 118 APPENDIX B COMPUTATIONAL WORKFLOW . . . . . . . . . . . . . . . . 120 vi APPENDIX C FRACTIONAL-ORDER DIFFERENTIAL OPERATORS . . . . 121 vii CHAPTER 1 INTRODUCTION 1.1 Motivation The excessively high computational cost of direct numerical simulations (DNS) of realis- tic turbulent flows has inspired many researchers to use coarse-grained techniques including Reynolds-averaged Navier–Stokes (RANS) and large-eddy simulation (LES) methods. Us- ing Reynolds averaging approach in the Navier–Stokes (N-S) equations provide temporally averaged quantities, however, in LES approaches we are specifying a model, which accounts for the effects of the finer scales. In other words, we only solve the N-S equations for the large-scales and model the effects of the smallest scales on the larger scales. Turbulence experimental and DNS features have confirmed that the turbulence is intrin- sically nonlocal and its statistics are non-Gaussian, which means velocity increments have sharp peaks, heavy-skirts and also skewed[1, 2], however, the most of the turbulence models have been built based on the Boussinesq’s turbulent viscosity concept, in which one assumes turbulent stress tensor is proportional to the local mean velocity gradient at any point, and the proportionality coefficient is set to the turbulent viscosity. Prandtl in 1942 aimed to disregard this local constraint by introducing the extended mixing length concept for the first time. The new model was a migration from locality to nonlocality, however, the model and its implementation were not remarkably successful. Afterward, Prandtl parametrized the model in a way that the mixing length was taken to be bigger than the differential length. This strategy was the same as adding a weak nonlocal concept to the model,hence called weak nonlocal in the sense that the stress term is still in the form of Boussinesq and a collinear relation exists with the strain rate tensor at the same point. Bradshaw [3] in 1973 showed that Boussinesq’s hypothesis fails over curved surfaces and noted that form of the stress-strain relations is responsible for this failure. It should be mentioned that there were some important developments mostly based on polynomial series compared to the Boussi- nesq type modeling including the works done by Spencer and Rivlin [4, 5], Lumley [6] and 1 Pope[7]; however, they all lack the accuracy that a “true” physical modeling should provide especially for the second-order and higher tensor series development. Using fractional derivatives is a relatively convenient approach to bring in nonlocal- ity concept from mathematical point of view. Fractional operators represent the heavy- tailed stochastic processes, which can be properly utilized in incorporating the long-range interactions in various mathematical models including but are not limited to beam vibra- tion analysis and damage modeling considering memory effects [8] and visco-elasto-plastic models [9]. Moreover, Harnessing the fractional models capabilities can be obtained prop- erly using the highly accurate numerical schemes for integer and fractional order PDEs [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],which is also an active research topic. In a pioneer work by Hinze et al. [21] in 1974, the authors described the memory effect in a turbulent boundary layer flow. They used the experimental data, produced downstream of a hemispherical cap attached to the lower wall of channel geometry and illustrated that when one computes eddy-viscosity using Boussinesq’s theory in the lateral gradient of the mean flow and turbulence shear-stresses, there is a significant non-uniform distribution that also exists in the outer region of the boundary layer. Interestingly, a nonlocal expression for the gradient of the transported field was proposed in a novel approach by Kraichnan in the same year [22] for the scalar quantity transport. Afterward, fractional-order models based on the RANS approach were offered in [23, 24, 25, 26, 27]. One of the main con- tributions in the development of nonlocal RANS models is by Egolf and Hutter [25, 28]. They started from Lévy flight statistics and generalized the zero-equation local Reynolds shear stress expression to a nonlocal and fractional type. The method is based on Kraich- nanian convolution-integral approach and utilizing different weighting functions. Using the mentioned weighting functions, one can make a bridge between the first-order gradient of the common eddy diffusivity models and the mean velocity difference term. Their proposed model is based on the four distinct steps that can be followed to change a local operator to a nonlocal one. In reality, the final proposed model is a more general and extended version 2 of Prandtl’s zero-equation mixing length and shear-layer turbulence models. The proposed model is called the Difference-Quotient Turbulence Model (DQTM). There is also emerging attention for the nonlocal LES-based models. Samiee et al. pro- posed a new model for homogeneous isotropic turbulent (HIT) flows based on the fractional Laplacian by employing Lévy stable distribution and tempered Lévy stable distribution in kinetic level [29, 30]. They showed that the new nonlocal proposed models potentially re- cover the non-Gaussian statistics of subgrid-scale stress motions. Laval et al. [31] analyzed the effects of the local and nonlocal interactions on the intermittency corrections in the scal- ing properties. They observed that nonlocal interactions are responsible for the creation of the intense vortices and on the other hand, local interactions are trying to dissipate them. Inspired by the mentioned observations, they came up with a new turbulence model that accounts for both the local and nonlocal interactions for the study of intermittency. In their proposed model, the large and small scales are being coupled by nonlocal interactions using a multiplicative process and additive noise along with a turbulent viscosity model for the local interactions. The results of the new model qualitatively cover the previously observed anomaly and intermittency aspects. Akhavan-Safaei et al. [32] proposed a fractional LES approach for the subgrid-scale modelings of the scalar turbulence. They utilized the two- point statistics for defining the optimal fractional order of the new nonlocal model and by using a priori assessment they showed that there is better agreement between the probability distribution function (PDF) of the SGS dissipation and the one that comes from the filtered DNS data. Turbulence fractional modeling in wall-bounded turbulent flows has been done by Keith et al. [33]. They started with consideration of the spectral velocity tensor based on ex- perimental data sets. Considering different scenarios, different shapes for different types of boundary conditions and energy spectra including the energy-containing and dissipative ranges can be imagined. They tuned their model parameters for the shear-free boundary layer and then the calibrated model was utilized to rendering the velocity field for a uniform 3 shear boundary layer. Leoni et al. [34] proposed a new nonlocal eddy viscosity-based model that can be applied in both the isotropic and anisotropic turbulent flows. They assessed model performance using two-point statistics of the filtered quantities for the homogeneous isotropic turbulence and channel flow canonical test cases. There are also some other re- lated studies that one can consult with, including hybrid nonlocal model in the case of magnetically confined plasma [35], generalization of a deconvolution model with fractional regularization for the rotational Navier-Stokes equations [36], fractional-Laplacian closure and its connection to Richardson pair dispersion [37]. Considering the nonlocal models in the literature, there are some important imperfec- tions including high sensitivity to the model constants, relatively low correlation coefficients between the stresses or fluxes of the model and the ground-truth ones, and no back-scatter prediction of kinetic energy from small scales to large scales. However, the conventional and mainly utilized local LES turbulent models are being improved over time to be free from mentioned deficits. One of the methods in the improvement process is using the dynamic procedure for the determination of the model constant [38]. To fill the gap in the litera- ture and provide an applicable and relatively easy to implement nonlocal LES model, we have designed a new dynamic nonlocal model that accounts for the all mentioned downsides. In the new dynamic fractional subgrid-scale models (D-FSGS, DNPS), two nonlocality and dynamic features have been leveraged together for the first time. This coupling between two important features provides a unique and higher performance than the local dynamic or static nonlocal models. Interestingly, the analysis showed that in the new models we have remarkably less sensitivity to the fractional-order, which is needed to be specified in nonlocal models. This relative freedom is obtained thanks to the novel coupling between the dynamic procedure and the nonlocal nature of the base model. We derived and implemented the D-FSGS model in both a priori and a posteriori stringent tests and compared the results with the conventional local models including Smagorisnky (SMG) and dynamic Smagorinsky (D-SMG) models. 4 1.2 Outline of the Study As the overall objective of this research, first, we are showing the emergence of nonlocal features in the transport phenomena especially in the turbulent flows. Second, we are devel- oping new generations of models that incorporate these inherent features. In the following, we concisely introduce these contributions: Chapter 2: We study the stochastic modeling of a rotating internal cylinder flow subject to uncer- tain rotational effects. In this part, we used high-fidelity simulation techniques to precisely address the anomalous transport features in such a system through statistical analysis of the expected field variables. The main sections of this study can be classified as follows: • Starting by formalization of the stochastic Navier-Stokes equations with the random symmetry-breaking inputs and boundary conditions, we leveraged the spectral element method (SEM) as our solver along with the probabilistic collocation method (PCM) to construct a stochastic computational framework. • A global sensitivity analysis was performed in order to reducing the random space dimension. Therefore, the only dominant stochastic directions were kept. Using the selected stochastic directions, we computed the expected velocity fields. Temporal study of the velocity fluctuations revealed the volution of their probability distribution function (PDFs) and these statistical quantities reflect the instability dynamics and anomalous transport features. • Expected vorticity fields were computed in the next step to have a picture of the non-Gaussian statistical behavior at the onset of flow instability. These instabilities and vortices are showing a high potential to increase the memory effects in the hy- drodynamics which can be utilized as a mixing amplification factor in the engineering systems design. 5 Chapter 3: A novel dynamic nonlocal turbulence model was developed for the large and very large eddy simulations. Starting from the filtered Boltzmann equation, we could build a Germano identity in the divergence form of fractional Laplacian and propose a first of its kind dynamic nonlocal turbulence model. The model was tested in both a priori and a posteriori tests. • We developed three different data-driven approached for the determination of the frac- tional component. Therefore, the final dynamic model was totally automatic and free of need to the preknowledge for the tinning parameters. • The a priori results show that the dynamic model adds some important features to the base static model. These advantages are higher performance in the sense of correlations between the all stress components, forces and dissipation predictions. Moreover, it adds the back-scatter prediction to the base static model which is a intrinsic property of the turbulent flows. • The new dynamic model is capable to maintain the high-order structures in both the LES and VLES regions. This is very important achievement, since in the environmental applications of turbulence modeling, one can use significantly bigger filter sizes and reduce the total computational costs remarkably. Chapter 4: Observations of scalar turbulence on small scales have been found to exhibit anomalous diffusion characteristics and higher levels of intermittency than that of fluid turbu- lence. Consequently, scalar turbulence dynamics modeling at a subgrid level is a challenging problem. In this chapter we are introducing a new dynamic nonlocal model for the passive scalar turbulence. This model has been developed based on the Germano’s identity idea on the scalar flux in the Laplacian form. We tested the newly developed model in both a priori and a posteriori fashions and compared the obtained results versus the conventional static Prandtl-Smagorinsky and dynamic Prandtl-Smagorinsky and static nonlocal model. This 6 model is totally different from the one that we developed in the previous chapter, since it is based on the fluxes rather than the divergence of the fluxes. Also, in the physics of this problem we have a mean gradient in the second direction which challenges any model remark- ably. Despite strong anisotropies, our proposed subgrid-scale nonlocal turbulence model can be fit into coarse LES and VLES frameworks, and hence can be applied to environmental problems. • With our high-fidelity datasets pertaining to forced homogeneous isotropic turbulence, we examined the effect of fractional order and characteristic filter sizes in LES and VLES cases. • During enough large-eddy turnover times, we utilized ensemble-averaged quantities from ten separate three-dimensional snapshots to make final decisions. When the groundtruth force and the predicted SGS force were most correlated, the optimal frac- tional order was selected for each scenario. • we observed that the time-averaged LES solution obtained from utilizing the DNPS model is performing remarkably successful in maintaining a low-level error over the multitude of scales for spatial shift. Therefore, we realized that unlike the PSM and DPSM model, the DNPS model does a great job in prediction of the high-order and two-point (nonlocal) statistics of filtered scalar field especially in over the inertial- convective subrange. Chapter 5: A conclusion to the dissertation proposal describes the findings and possible future works. Moreover, we present the preliminary results of our newly developed model, dynamic tempered model. 7 CHAPTER 2 ANOMALOUS FEATURES IN INTERNAL CYLINDER FLOW 2.1 Background 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 processes that exhibit non-Markovian (long-range memory) effects, non-Fickian (nonlocal) interactions, non-ergodic statistics, and non-equilibrium dynamics [39]. 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 semi- conductors, 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 [40, 41]. 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. Flow within and around cylinders is a rich physical problem that involves complex geometry and nonlinear flow instabilities, with unsolved questions on flow/vortex structures and anomalous turbulent mixing [42]. 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 [43] 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 [44] and Griffin and Ramberg [45], respectively. These studies are categorized as external flows around cylinders and some significant contributions in this regard may be found in [46, 47, 48, 49, 50]. However, flow inside systems with fast rotation including cylinders, squares and annulus geometries are also of great importance. Turbo- 8 machinery, mixing process, gravity-based separators, geophysical flows, and journal bearing lubrication are all clear examples for these types of internal flows [51, 52]. 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 [53, 54, 55]. 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 [56, 57, 58] and numerically [59, 60, 61, 62]. In such problems, emergence of the adverse angular momentum is an important mechanisms, which initiates flow instability. More specifically, Lopez et al. [63] 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 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 [64], oscillating sidewalls [65] and, precessional forcing [66]. 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 simulation 9 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 stochasticity, (i.e., sources of disturbance). This urges for another level of modeling and investigation, 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 uncer- tainty 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 stochastically modeled once we obtain additional information about the system [67]. Uncertainty in modeling procedure and also inaccuracy of the measured data are two main factors in arising epistemic uncer- tainty. The uncertainty in modeling could be the result of a variety of possibilities including the effects of geometry [68, 69, 70, 71], constitutive laws [72, 73, 74, 75, 76], rheological models [77, 78, 79], low-fidelity and reduced-order modeling [80, 81, 82, 83, 84, 85, 86, 87], and random forcing sources in addition to the random field boundary/initial conditions [88, 89, 90, 91, 92, 93]. In the current work, we seek to fill a gap in the rich literature of investigating flow instabilities inside 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 computational 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 em- ploy spectral element method (SEM) along with the probabilistic collocation method (PCM) to formulate a stochastic computational framework. 10 • 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 structures”. These vortices increase the memory effects in the hydrodynamics and we characterize their impact as long-time “anomalous” time-scaling of enstrophy leading to effective enhancement in the mixing capacity of the system. The structure of the rest of this chapter is outlined as follows: In section 2.3, we formulate the stochastic version of the Navier-Stokes equations for incompressible flows and develop our stochastic modeling procedure. In section 2.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 2.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. 11 2.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 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. 2.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 12 equation, ∇ · V = 0, for Newtonian viscous fluids ∂V + V · ∇V = −∇p + ν ∇2 V , ∀(x, t; ω) ∈ Ω × (0, T ] × Ωs , (2.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 time-dependent wall angular velocity, θ̇(t; ω) = cos (α(ω)t) e−λ(ω)t , ∀(t; ω) ∈ (0, T ] × Ωs , (2.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 , (2.3) ∥x∥2 = 1, ∥rϵ(ω) ∥2 = ϵ(ω), 13 where ∥ · ∥2 denotes the L2 norm. Moreover, according to (2.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 , (2.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 [94], 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 rather work in the target space. Finally, the formulation of stochastic governing equations in (5.1) subject to the boundary/initial conditions in equation (2.3) can be posed as: Find V (x, t; ξ) : Ω × (0, T ] × Ψ → R such that ∂V + V · ∇V = −∇p + ν∇2 V , (2.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. 2.4 Stochastic Computational Fluid Dynamics Framework Discretization of Physical Domain and Time-Integration Spectral/hp element method [95] is a high-order numerical method to discretize the gov- erning equations (5.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. In 14 SNel e 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), (2.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 [95]. NEKTAR++ [96, 97], a parallel open-source numerical framework, provides a seamless plat- form offering efficient implementation of multiple SEM-based solvers in addition to the 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 [96]. 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 2.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 [98, 95] 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 [99], to construct the geometry and then the h-refined grid. The generated grid is illustrated in Figure 2.1a, which shows elemental nodes and 15 10−1 Error 10−2 10−3 0 200000 400000 600000 800000 1000000 1200000 DoF (a) (b) Figure 2.1: (a) Constructed structured grid with transitional h-refinement. (b) Grid conver- gence study based on the error in kinetic energy. h-refinement near the wall. For this h-refined grid, we employ a spatially-variable polyno- mial expansion [97] 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 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 2.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 A). 16 Stochastic Discretization Sampling from the parametric random space introduced in section 2.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) [100, 101, 102], 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. [103] 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 }J j=1 . In other words, we use the SEM setup described in section 2.4 to solve a set of deterministic problems in which the wall velocity field V ∂Ω (x, t; ξ) in equation (2.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 (ξ). (2.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ξ. (2.8) Ψ 17 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. [103] coincide these quadrature points with associated integration weights {wj }J j=1 , one can efficiently compute the approximation to the integral in equation (2.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 ). (2.9) j=1 In equation (2.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 ) = ( 12 )3 . Consequently, the approximate computation of the expectation integral (2.8) is simplified to J 1X E[V (x, t; ξ)] ≈ wj V (x, t; q j ). (2.10) 8 j=1 Similar to the MC approach and using (2.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; ξ)] . (2.11) 8 j=1 2.5 Variance-Based Sensitivity Analysis Grasping knowledge on the significance of sources of randomness in a stochastic modeling procedure could be very helpful in terms of reducing the computational cost and also deci- sion 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 [104, 105]. In practice, sensitivity of the QoI to each stochastic parameter is measured by the condi- tional variance in the QoI, which is caused by that specific parameter. In general, for a 18 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 [106], 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 , (2.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 24 (a) (b) Figure 2.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). 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 2.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. [111] studied the behavior time-averaged root- 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. [112] 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 25 part could be decomposed into V = V + V ′, (2.18) where V ′ represents the fluctuations of V and V denotes its ensemble average. Unlike the applied approach in [111, 112] 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 large enough number of realizations. Thus, what we obtain as the result of equation (2.10) is the representation of ensemble mean in a PCM setting [113]. The stochastic convergence analysis we performed in section 2.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; ξ)] . (2.19) According to the sensitivity analysis we performed in section 2.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 2.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 26 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 2.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 (2.18). The procedure of computing the fluctuations from SEM-based realizations is briefly explained in Appendix B. Figure 2.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 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 2.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 2.8a and 2.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 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 27 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 2.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. and asymmetries quickly leads to a highly non-Gaussian statistical behavior (see Figures 2.8b and 2.8d and compare with the standard Gaussian PDF). More specifically, Figure 2.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 2.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 28 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 2.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 2.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. 29 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 2.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 2.9a and 2.9b show that the non-zero skewness 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 2.8b. For u′θ , Figure 2.9a illustrates non-zero skewness factor values close to the wall at all the recorded times and Figure 2.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 2.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 2.9b and 2.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 30 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 2.9d). Finally, by comparing the PDFs of velocity fluctuations for 0.2 < t with the one associated with standard Gaussian (see Figure 2.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 [114] 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 failed to be addressed by Batchelor [115]. 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 [40]. 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 connection 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 2.11 provides this statistical information at t = 0.25, 0.5, and 0.75. Comparing the vorticity PDFs shown in Figure 2.11a to the standard Gaussian PDF makes it evident 31 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 2.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. that fingerprints 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 2.11b and 2.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 32 105 E(t) ∼ t−2 (I) ∼ t−1 ∼ t−1/2 Enstrophy 104 (III) 103 (II) 102 100 101 t (a) (b) Figure 2.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. 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 [116, 117, 118], subgrid-scale (SGS) motions and their nonlocal modeling for homogeneous and wall-bounded turbulent flows [119, 120], boundary layer flows [121, 122], and turbulent flows interacting with wavy-like moving/actuated surfaces (with application to reduction and control of flow separation) [123, 124]. Here, an interesting yet, practical question that could be raised is that if such “intensified” 33 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 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 [40]. 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 (2.16), we define the enstrophy, E(t), in our problem setting as Z 1 2 E(t) = ωz′ (x, t) dΩ. (2.20) µ(Ω) Ω By computing the record of enstrophy for relatively long times (obtained from the same flow 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 [88] 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 [100] 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. [102] 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 34 described in section 2.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 2.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 2.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 stage (III) in Figure 2.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. 35 CHAPTER 3 A DATA-DRIVEN DYNAMIC NONLOCAL SUBGRID-SCALE MODEL FOR TURBULENT FLOWS 3.1 Introduction The prohibitively high computational cost of direct numerical simulations (DNS) of re- alistic turbulent flows has motivated the community of research in turbulence to develop coarse-grained techniques including Reynolds-averaged Navier–Stokes (RANS) and large- eddy simulation (LES) methods to reduce the intractably large degrees of freedom in DNS studies[125, 126, 127, 38, 128, 129]. Using Reynolds averaging approach in the Navier–Stokes (N-S) equations provide temporally averaged quantities, however, in LES approaches one em- ploys a subgrid-scale (SGS) model, which represents the effects of the finer scales[130, 125]. Turbulence experimental and DNS features have confirmed that the turbulence is intrin- sically nonlocal, and its statistics are non-Gaussian, which means velocity increments have sharp peaks, heavy-skirts and also skewed[1, 2]. However, most of the turbulence models have been built based on the Boussinesq’s turbulent viscosity concept, in which one assumes turbulent stress tensor is proportional to the local mean velocity gradient at any point, and the proportionality coefficient is set to the turbulent viscosity. Prandtl in 1942 aimed to disregard this local constraint by introducing the extended mixing length concept for the first time. The new model was a migration from locality to nonlocality, however, the model and its implementation were not remarkably successful since the scale of nonlocality was comparable with the differential length scale. Afterward, Prandtl parametrized the model in a way that the mixing length was taken to be bigger than the differential length. This strategy was the same as adding a weak nonlocal concept to the model, hence called weak nonlocal. Bradshaw [3] in 1973 showed that Boussinesq’s hypothesis fails over curved surfaces and noted that form of the stress-strain relations is responsible for this failure. It should be mentioned that there were some important developments mostly based on polynomial series compared to the Boussinesq type modeling including the works done by Spencer and 36 Rivlin [4, 5], Lumley [6] and Pope[7]; however, they all lacked the accuracy that a “true” physical modeling should provide especially for the second-order and higher tensor series development. Using generalized-order derivatives is a relatively convenient approach to bring in nonlo- cality concept from mathematical point of view. The generalized-order operators represent the underlying heavy-tailed stochastic processes at the continuum level, which can be prop- erly utilized in incorporating the long-range interactions in various mathematical models including but are not limited to beam vibration analysis [131], anomalous rheology model- ing [132, 133, 134], damage modeling considering memory effects [8] and visco-elasto-plastic models [9]. Moreover, harnessing the generalized-order models capabilities can be obtained properly using the highly accurate numerical schemes for integer and fractional-order PDEs [10, 11, 12, 13, 14, 15, 16, 17, 135, 18, 19, 20], which is also an active research topic. In a pioneer work by Hinze et al. [21] in 1974, the authors described the memory effect in a turbulent boundary layer flow. They used the experimental data, produced downstream of a hemispherical cap attached to the lower wall of channel geometry and illustrated that, when one computes eddy-viscosity using Boussinesq’s theory in the lateral gradient of the mean flow, there is a significant non-uniform distribution that also exists in the outer region of the boundary layer. Interestingly, a nonlocal expression for the gradient of the trans- ported field was proposed in a novel approach by Kraichnan in the same year for the scalar quantity transport [22]. Afterward, fractional-order models based on the RANS approach were offered in [23, 24, 25, 26, 27]. One of the main contributions in the development of non- local RANS models is by Egolf and Hutter [25, 28]. They started from Lévy flight statistics and generalized the zero-equation local Reynolds shear stress expression to a nonlocal and fractional type. The method is based on Kraichnanian convolution-integral approach and utilizing different weighting functions. Using the mentioned weighting functions, one can make a bridge between the first-order gradient of the common eddy diffusivity models and the mean velocity difference term. Their proposed model is called the Difference-Quotient 37 Turbulence Model (DQTM) and is based on the four distinct steps that can be followed to change a local operator to a nonlocal one. In reality, the proposed model is a more general version of Prandtl’s zero-equation mixing length and shear-layer turbulence models. There is also emerging attention for the nonlocal LES models. Samiee et al. proposed a new model for the HIT flows based on the fractional Laplacian by employing Lévy stable and tempered Lévy stable distributions in kinetic level [29, 30]. They showed that the new non- local models can recover the non-Gaussian statistics of subgrid-scale stress motions. Laval et al. [31] analyzed the effects of the local and nonlocal interactions on the intermittency corrections in the scaling properties. They observed that nonlocal interactions are respon- sible for the creation of the intense vortices and on the other hand, local interactions are trying to dissipate them. Akhavan-Safaei et al. [32] proposed a fractional LES approach for the subgrid-scale modelings of the scalar turbulence. They utilized the two-point statistics for defining the optimal fractional-order of the new nonlocal model, and by using a priori assessment they showed that there is proper agreement between the probability distribution function (PDF) of the SGS dissipation and the one that comes from the filtered DNS data. Harmonious with this study, Akhavan-Safaei and Zayernouri developed a corresponding non- local spectral transfer model and a new scaling law for scalar turbulence in [136]. Their new analysis additionally reconciled the close similarities between this work and their earlier development in [32] when the filter scale approaches the dissipative scales of turbulent trans- port. There are also other related studies that one can consult with, including preliminary fractional modeling in wall-bounded turbulent flows [33], a priori survey of nonlocal eddy viscosity-based model for the isotropic and anisotropic (channel flow canonical test cases) turbulent flows [34], hybrid nonlocal model in the case of magnetically confined plasma [35], generalization of a deconvolution model with fractional regularization for the rotational N-S equations [36], and fractional Laplacian closure and its connection to Richardson pair dis- persion [37]. Going even beyond the scope of research in turbulence, a new comprehensive survey on the nonlocal models for several crucial applications, including anomalous subsur- 38 face dynamics, turbulence modeling, and extraordinary materials, was recently performed by Suzuki et al. [137]. Considering the nonlocal models in the literature, there are some important imperfections including the sensitivity to the model constant and fractional-order parameter, relatively low correlation coefficients ρ(τ ∆ , τ ∆,M odel ), and no back-scatter prediction of kinetic energy from small scales to large scales. However, the conventional and frequently utilized local LES turbulent models are being improved over time to be free from mentioned deficits. One of the methods in the improvement process is using the dynamic procedure for the determination of the model constant [38]. To fill the gap in the literature and provide an applicable and relatively easy to implement nonlocal LES model, we have developed a new dynamic nonlocal model that accounts for all the aforementioned downsides. In the new dynamic fractional subgrid-scale model (D-FSGS), both nonlocality and dynamic features have been leveraged together for the first time. This match between two important features provides a unique and higher performance than the dynamic local or static nonlocal models. Interestingly, the analysis showed that in the new model, we have remarkably less sensitivity to the fractional-order, which is needed to be specified in nonlocal models. This relative freedom is obtained thanks to the novel coupling between the dynamic procedure and the nonlocal nature of the base model. In the following, we derived and implemented the D-FSGS model in both a priori and a posteriori stringent tests and compared the results with the conventional local models including Smagorisnky (SMG) and dynamic Smagorinsky (D-SMG) models. This study is organized as follows: in section 5.2 we talk about the governing equations and development of the new model. Section 3.3 starts with introducing three main distinct approaches for the determination of the fractional-order, and then we do a comprehensive a priori assessment along with a comparative study with conventional models. In section 3.4, we study the numerical stability of the proposed model and do different tests including two- point diagnostics to have a complete overview of the model performances in a a posteriori 39 sense. 3.2 Model Development Implementation of a low-pass filter on the N-S equation forms a closure term on the right-hand side of the momentum equation, which needs to be modeled, representing the (unknown) SGS dynamics. The filtered incompressible N-S equations can be written as ∂ ū 1 + ū · ∇ū = − ∇p̄ + ν ∇2 ū − ∇ · τ, (3.1) ∂t ρ ∇ · ū = 0. In this equation, ū is the filtered velocity vector, ρ denotes the constant density, p̄ represents the filtered pressure, and ν shows the kinematic viscosity. The effects of the small scales arise in the so-called SGS stress tensor τij = ui uj − ūi u¯j , (3.2) forming the closure term to be modeled. Several models have been proposed during the last decades to close the filtered N-S equation in both functional and structural LES models [125, 127, 126]. In the classical Smagorinsky model (SMG)[138], the deviatoric part of the stress tensor is written as SM G = −2 ν S¯ , τij (3.3) t ij ∂ u¯ ∂ u¯ where νt denotes the eddy viscosity, and S¯ij = 12 ( ∂x i + ∂xj ) represents the filtered (resolved) j i strain rate tensor. In this model, νt is constructed based on the Prandtl’s Mixing Length hypothesis, νt = Cs L2δ |S̄|, (3.4) q where Lδ shows the effective grid scale, and |S̄| = 2S¯ij S¯ij exhibits the magnitude of the resolved scale strain rate tensor [127]. One of the main drawback of this model is its flow- dependent feature, which means that the model can not correctly predict final quantities with 40 a single universal constant in different scenarios such as shear flows, wall-bounded flows, or transitional flows. To overcome this challenge in the Smagorisnky model, Germano et al. [38] proposed a novel procedure for evaluation of the model coefficient, which is called the dynamic Smagorinsky model (D-SMG). The suggested procedure was a breakthrough in turbulence modeling, and several researchers utilized the same concept afterward [139, 140, 141, 142]. The dynamic procedure was designed based on the classic idea that one can extract useful information from the smallest resolved scales for modeling the subgird-scales; however, the way it was applied was indeed novel. It calculates the eddy-viscosity coefficient locally for each LES grid point as the calculation progresses, and there is no need for any predefined inputs in the model. This model was constructed based on the scale-invariance hypothesis and calculates the model constant using the information from the resolved section. Fractional SGS Model [29]: We have recently developed a nonlocal SGS model that we present here briefly for making this work self-contained. Starting from the Boltzmann kinetic level description of the flow, we employ the BGK model for the collision of the particles and have ∂f f − f eq + u · ∇f = − , (3.5) ∂t T where f = f (x, u, t) is called the single-particle probability distribution function and shows the particles density in the phase space (x, u) at time t. Moreover, f eq represents the local equilibrium distribution function and is defined based on the Maxwell distribution. ρ f eq (K) = F (K), (3.6) U3 |u−V |2 where F (K) = e−K/2 , K = 2 and U corresponds for the thermal agitation speed. U The left-hand side of Eq. 5.5 correspond to the streaming of the non-reacting particles and right-hand side shows the collision operator with a relaxation time T . One can solve Eq. 5.5 analytically by method of characteristic to find distribution in terms of equilibrium state 41 [143] Z ∞ f (t, x, u) = e−s f eq (x − uτ s, u, t − τ s) ds. (3.7) 0 Moments of the f would provide the macroscopic flow variables. Therefore, one can write R R ρ = f (t, x, u)du, and ρV (t, x) = uf (t, x, u)du to compute density and fluid velocity. Incorporating filtering procedure into the Eq. 5.5 would provide filtered Boltzmann equation as ∂ f¯ f¯ − f eq (∆) + u · ∇f¯ = − . (3.8) ∂t T However, the collision term is highly-nonlinear and filtering kernel cannot commute [144]. |u−V |2 Hence, a closure problem would be built by defining K = since U2 f eq (K) ̸= f eq (K). (3.9) One of the approaches for handling this problem is using a power-law distribution for model- ing f eq (K)−f eq (K). In this method, for closing the filtered collision operator, a heavy-tailed distribution function is used to account for the nonlinear effect of turbulence, i.e., f eq (K) − f eq (K) ≃ f M odel (K̄) = Dβ f β (K̄), (3.10) in which f β (K̄) = ρ3 F β (K), where F β (K) represents an isotropic Lévy β-stable distribution. U More details and discussion can be found in [29, 30]. Also, Dβ is a real-valued constant which is going to be addressed in the present dynamic model. In the fractional Laplacian model, the SGS forces are defined as (∇.τ )i = να (−∆)α ūi , α ∈ (0, 1], (3.11) where the filtered velocity ūi : Rd × (0, T ] → R, where d = 3 is the dimension of physical domain and T represents the simulation time, in addition, (−∆)α (·) denotes the space- fractional Laplacian of order 2α ∈ (1, 2], which can be defined as a singular integral operator 42 given by Z ūi (x) − ūi (y) (−∆)α ū i = cd,α dy, (3.12) |x − y|d+2α Rd 4α Γ(d/2+α) where cd,α = d/2 . It should be mentioned that under some regularity assumption π |Γ(−α)| for the velocity field (which all make sense in our problem), one can readily show that the fractional Laplacian can be written in divergence form via the Riesz transform. That can be employed in the Gauss’s divergence theorem to render the volume integral of this new representation as a surface integral, ensuring the conservation of momentum. Moreover employing the periodic boundary conditions, the corresponding Fourier transform of the fractional Laplacian in (5.12) is given by (see e.g., [145]) h i F (−∆)α ūi (x) = |k|2α F[ūi ](k), (3.13) in which F and k are the Fourier transform and Fourier numbers, respectively. Evidently, the integer-order Laplacian operator is recovered simply by putting α = 1. The Fourier transform 3.13 provides a rather convenient way of handling the fractional operators Fourier space and in our Fourier spectral method for simulating the problem. The corresponding “eddy-viscosity”-like model coefficient in (5.11) is then obtained as να = C F (α), (3.14) 2α where C represents an up-scaling model input being proportional to U T 2α−1 from the kinetics theory’s perspective, yet to be dynamically computed in the subsequent (continuum- level) simulations, moreover, F denotes a deterministic univariate function of the fractional- order α, explicitly given by 22α Γ( 2α+3 2 ) Γ(2α + 1) , F (α) = 3 (3.15) π 2 |Γ(−α)| rendering (5.11) as (∇.τ )i = C F (α) (−∆)α ūi , α ∈ (0, 1]. The term F (α) is introduced just for the sake of the convenience and the reader is referred to [29] for more details. In Appendix 43 A, we have provided the preliminary mathematical concepts and definitions. In what follows and in this generalized order context, we develop a new dynamic procedure to automatically compute C from data on-the-fly. Derivation of the Nonlocal Dynamic Model Procedure: We write (∇.τ )i = C1 F (α)(−∆)α ūi for the sake of simplicity. Implementing the second filtering process (test-level filter), gives the divergence of the SGS stresses at the test-level filter (subtest-scale stress) as (∇.T )i = C2 F (α)(−∆)α ub̄i . (3.16) In the above equation, (·) c indicates the test-level filtering, which is commonly chosen as twice the grid-level filtering. Now, the Germano identity [38], which relates the stresses at grid-level (τij ) and test-level (Tij ), is employed to make a bridge between the resolved scales and the subgrid-scales, Gij = ūd i ūj − ub̄i ūbj = Tij − τc ij , (3.17) which is a known quantity, and represents the resolved turbulent stress. In the divergence form, we have  ∇. ūd i ūj − ub̄i ūbj = (∇.G)i . (3.18) Therefore, we construct the Germano identity in the divergence form based on the previously introduced fractional Laplacian model. The Germano identity extrapolates and parametrize the model constant for the subgrid-scale part using the information in the smallest resolved scales,  V  (∇.T )i − (∇.b τ )i = C F (α)(−∆)α Vb̄i − F (α)(−∆)α V̄ i . (3.19) Here, we assume the model constant C being scale invariant around the grid and test filter size, (C = C1 = C2 ), and C is in general a space-time dependent tensor, yet, it has been assumed to be spatially uniform. One needs to solve Eq. 5.18 for C, which forms a tensorial 44 10 101 ( )u 1 100 10 2 10 1 ( )u 10 3 10 2 PDF 10 3 10 4 u 10 4 x 10 5 u 10 5 x 10 6 10 6 15 10 5 0 5 10 15 10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0 (a) (b) Figure 3.1: Difference between commuting of the low pass filter on the (a) integer order derivatives, and (b) fractional-order derivative (α = 0.6). equation. There are two common methodologies to get a scalar model constant. First, one may contract with S¯ij as mentioned originally in [38]. Second approach is proposed by Lilly [146], which is more preferable and commonly used. This approach is based on contracting with the tensor, which has been multiplied to the unknown coefficient in the right hand side of Eq. 5.18. One important point that should be mentioned in this section is that in the conventional local LES modeling, the main assumption is that the operation of filtering and integer-order differentiation commutes, then we get a set of filtered equations. Commuting means commutation between the filtered quantity and the spacial derivatives. However, the filtering procedure does not commute with the fractional-order operators in general. One simple example regarding this important matter is depicted in Figure 3.1 for α = 0.6, comparing the derivatives of the filtered velocity at test-level and filtered derivative for the integer and fractional-order cases. We have used Lilly’s approach for contracting purposes based on the Least Square Method that was described by Lilly in 1992 [146]. Therefore, one can write Eq. 5.18 as (∇.G)i − CNi )2 = e, (3.20) in which e is the squared error and the nonlocal term Ni is defined as V  Ni = F (α)(−∆)α Vb̄i − F (α)(−∆)α V̄i . (3.21) 45 Finding the unknown in minimized form is conveniently possible by putting the derivative 2 equals to zero, considering and testing ∂ e2 > 0 for the error minimization scope. Finally, ∂C the scalar model coefficient C can be computed dynamically as ⟨(∇.G)i Ni ⟩ C= . (3.22) ⟨Ni Ni ⟩ Numerical instability may be occurred due to the negative eddy-viscosity in prolonged pe- riods of time. As a remedy, one can perform an averaging over the directions of statistical homogeneity as suggested by Germano et al. [38] in the a posteriori assessments. Figure 3.2 illustrates the variations of model constants in D-SMG and D-FSGS models in the imag- inary center-line of a periodic domain in the first direction to have a comparison between the model constant variations in the context of an example. It is worth mentioning that in inhomogeneous flows we can not use the Eq. 3.13 since it is only correct for problems with periodic boundary conditions. However, there are other definitions of fractional Lapla- cian operator without this restrictions [145]. Moreover, the derivation of the nonlocal model from the Boltzmann-BGK model would be different and more complicated. For instance, in the boundary layer case, the analytical solution for the Boltzmann-BGK model would be different than the one that has been used in the derivation of the base model, FSGS. 3.3 A Priori Analysis We assess the performance of the D-FSGS and compare the results with the results of conventional LES models including SMG and D-SMG. Also, we discuss about the perfor- mance of the D-FSGS model using different filter sizes (covering both the LES and VLES regions). Before starting the a priori tests, we need to study the effects of the fractional-order parameter (α) and its proper selection for different scenarios. We utilize ten independent snapshots of forced HIT data distributed over enough eddy turnover times on a triply periodic domain, Ω = [0, 2π]3 , for each targeted Taylor Reynolds number. The covered Taylor Reynolds numbers are 170, 240, and 355, and the corresponding grid resolution for each case is 3203 , 5203 , 10243 , respectively. The turbulence dataset comes 46 1.5 1.0 0.5 Cs , C 0.0 0.5 1.0 1.5 D-SMG D-FSGS 0 20 40 60 80 100 120 x Figure 3.2: Comparing the model coefficients in dynamic Smagorinsky model (D-SMG) and dynamic FSGS model (D-FSGS) in the middle imaginary line in a periodic domain using Lδ = 4. from our open-access pseudo-spectral parallel code, which has been elaborated in [147]. The Inference of Optimum Laplacian-Order (αopt ) The D-FSGS model depends on the fractional-order, α. We provide three different data- driven approaches for the calculation of this parameter. First Approach: we perform a precise comparison between the obtained correlation coeffi- cients, considering the effects of all stress components in the SGS tensor and the ground-truth stresses that come from the DNS results. Subsequently, we find the αopt , where the maximum DN S , τ M odel ])⟩ is achieved, in which, ⟨·⟩ denotes the ensemble-averaged of ρij = ⟨Avg(ρ[τij ij on different snapshots of data, and for each snapshot of data, we consider the average of all corresponding tensor components. By sweeping the fractional-order values from zero to one and step size equals to 0.01, we determined the αopt for different characteristic filter sizes in three different Reynolds numbers. In Figure 3.3, we have shown the αopt values for different filter sizes for three different Taylor Reynolds numbers. Also, Lδ = ∆∗ ∆DN S , where ∆∗ is the ratio between the LES and DNS grids, and the DNS grid size is defined as ∆DN S = 2π N associated with N = 320, 520, 1024 for three different Reynolds numbers. These results are obtained based 47 on the ensemble-averaged quantities of ten snapshots of data for each Reynolds number to ensure meeting the necessary conditions regarding performing LES studies [148]. There is an interestingly noticeable trend between the filter size and the obtained αopt from the ensemble-averaged relation for each Reynolds number (see Figure 3.3a). Second Approach: we employ a logarithmic regression to come up with a correlation between the αopt and Lδ for each Reynolds number. The results are showing that αopt values obey a power-law form, which reads αopt = a Lδ b . (3.23) Increasing the filter size would incorporate more nonlocality, and that is the reason for the Table 3.1: Proper parameters for being used in second approach (Eq. 3.23) for determination of optimum Laplacian-order. Reλ a b 170 1.53 -0.62 240 1.39 -0.49 355 1.08 -0.30 increase in the error bars in the VLES section (see Figures 3.3b, 3.3c, and 3.3d ); however, all cases converge statistically to the realized power-law relations. Moreover, there is a direct relationship between the Taylor Reynolds number and the optimum fractional-order in each filter size. Using the above relation for finding the αopt would be the second method that can be applied conveniently and confidently since the coefficient of determination, R2 , for these equations are above 0.98, and indicates that the relation correlates very well with the measured data. The a, b values in Eq. 3.23 are being determined for different Taylor Reynolds numbers. Table 3.1 shows the proper parameters for being used in second approach of optimum Laplacian-order determination. The values in this table obtained based on the ensemble-averaged results of ten data snapshots for each Reλ . Moreover, interpolation methods can be conveniently utilized using this table to find the proper value of αopt for other Reynolds numbers. 48 0.9 Re =170 Re =240 0.8 Re =355 0.7 opt 0.6 0.5 0.4 0.3 4 8 10 16 20 26 32 (a) 1.0 1.0 25 Re =170 Re =240 0.7 Re =355 VLES Region 25 17.5 VLES Region 0.8 0.8 20 15.0 VLES Region 20 0.6 12.5 0.6 15 opt 0.6 15 0.5 10.0 0.4 10 0.4 10 7.5 LES Region LES Region 0.4 LES Region 5 5.0 0.2 5 0.2 0.3 2.5 0 2 3 4 6 8 10 16 2 4 8 10 20 26 4 6 8 16 20 30 (b) (c) (d) Figure 3.3: Effects of the filter size (Lδ ) and Taylor Reynolds number (Reλ ) on the optimum fractional-orders obtained by the first approach, (a) the big picture of ensemble-averaged results based on three datasets, (second row) amount of modeled turbulent kinetic energy corresponds to each filter size when (b) Reλ = 170, (c) Reλ = 240, and (d) Reλ = 355. Third Approach: we propose a third approach for obtaining the proper fractional-order for a quite wide range of (fine-to-coarse) LES region yet still in the a priori sense, shown in Figure 3.3. The right vertical axes (secondary axis) in Figures 3.3b, 3.3c, 3.3d show the volumetric averaged values of the modeled kinetic energy in percent (ϕ) at that specific filter size. To calculate this value, ϕ, we find the ratio between the resolved (Er ) and filtered (Ef ) kinetic energy as 1 Ef = u u, (3.24) 2 i i 1 Er = (ui ui − ui ui ). 2 Noting that almost all of the LES studies are designed to model about less than ten percent of the total kinetic energy, one can conclude based on the obtained plots for all Reλ numbers, 49 0.300 0.24 0.275 0.22 0.250 0.20 0.225 0.18 0.200 0.175 0.16 0.150 SMG 0.14 SMG D-SMG D-SMG 0.125 FSGS 0.12 FSGS D-FSGS D-FSGS 0.100 11 22 33 12 13 23 11 22 33 12 13 23 Stress Components Stress Components (a) (b) Figure 3.4: A priori assessment in the context of the ensemble-averaged correlation coefficient comparisons in different models for (a) Lδ = 4, and (b) Lδ = 8 . taking α = 0.5 could be a representative corresponding factor for the LES region. We will show that even with this rough estimation, the proposed model provides better correlations in all of the stress components than the other conventional models. Therefore, using the second or third approach, we will have a dynamic nonlocal model, which is totally free of any tuning parameter including fractional-order. Statistical Performance Assessment We use ten three-dimensional snapshots of velocity fields, ui (x), i = 1, 2, 3, corresponding to Reλ = 240 and perform a priori testing to show the model performance. The DNS datasets with Reλ = 240 and N = 520 will be utilized for the rest of the paper. We have depicted the ensemble-averaged correlation coefficients for all of the stress components in Figure 3.4 for the filter size equals to 4 and 8. We see a remarkable better performance in both filter sizes and all directions for the new proposed model, D-FSGS. Interestingly, the coupling between the new dynamic procedure and the nonlocal nature of the model helps significantly to compensate for the directions in which we have low correlations. We see that D-FSGS outperforms static nonlocal (FSGS) and dynamic local (D-SMG) models in both filter sizes and all directions. Its well-known that, dynamic local models do not necessarily increase the correlation coefficient remarkably. In contrast to the D-SMG model, 50 we see a significant raise in the correlation coefficient of the D-FSGS versus FSGS. Therefore, the new dynamic procedure not only adds the back-scatter prediction and being free from tuning constant capabilities, it also improve the performance as well. In another word, the dynamic procedure in the context of the nonlocal modeling seems to be more advantageous and fruitful than the local modeling. Q In Figure 3.5 we have plotted the = −⟨τij S¯ij ⟩, the SGS dissipation of kinetic energy Q for the same filter sizes mentioned previously. The negative values of mean we have back-scatter, where the kinetic energy is transferred from subgrid to resolved scales. As the SGS dissipation of the filtered DNS results (ground-truth) show, a significant amount of flow shows back-scatter phenomena. We see the new model, D-FSGS, and D-SMG models are the only ones that truly predict the back-scattering in energy flow. Moreover, comparing the PDFs of dissipation show that the D-SMG model over-predicts the values for the tails of the distributions, which can be a potential factor in facing unstable conditions in simulations with this model necessitating the averaging operation in a posteriori tests[127]. However, the new model provides closer results to the DNS ones. It should also be mentioned that a comparison of the percentage of grid points with back-scatter in filtered-DNS, D-SMG, and D-FSGS revealed that there is a better agreement for the D-FSGS model in the prediction of this quantity. As has been reported by Piomelli et al [149], this percentage is a function of the filter type that one is using. In this study, three-dimensional box filtering for both the grid-level and test-level was utilized. The percentage of grid points with back-scattering in DNS was 26%, in D-SMG was 20% and in D-FSGS was 24% for Lδ = 4. For bigger filter size, Lδ = 8, these numbers were 30%, 18% and 25% for DNS, D-SMG and D-FSGS, respectively, which quantitatively indicates the better prediction of back-scatter for the D- FSGS model. Figure 3.6 shows the PDF of the ensemble-averaged SGS forces in three directions using Lδ = 4 and Lδ = 8. As it is clear in all three directions, there is superior match between the results of the D-FSGS model and DNS ones in comparison to the results of the Smagorinsky-based, and FSGS model considering tail and peaks of the SGS forces. 51 101 101 DNS DNS SMG SMG 100 D-SMG 100 D-SMG FSGS FSGS 10 1 D-FSGS D-FSGS 10 1 10 2 PDF 10 3 10 2 10 4 10 3 10 5 10 4 10 6 4 2 0 2 4 6 8 2 1 0 1 2 3 4 5 6 7 Sij ij Sij ij (a) (b) Figure 3.5: PDF of ensemble-averaged SGS dissipation of kinetic energy using different models for (a) Lδ = 4, and (b) Lδ = 8. 101 101 101 DNS DNS DNS 100 SMG 100 SMG 100 SMG D-SMG D-SMG D-SMG 10 1 FSGS 10 1 FSGS 10 1 FSGS D-FSGS D-FSGS D-FSGS 10 2 10 2 10 2 PDF 10 3 10 3 10 3 10 4 10 4 10 4 10 5 10 5 10 5 10 6 10 6 10 6 6 4 2 0 2 4 6 6 4 2 0 2 4 6 6 4 2 0 2 4 6 ( . )1 ( . )2 ( . )3 (a) (b) (c) DNS DNS DNS 100 SMG 100 SMG 100 SMG D-SMG D-SMG D-SMG FSGS FSGS FSGS 10 1 D-FSGS 10 1 D-FSGS 10 1 D-FSGS PDF 10 2 10 2 10 2 10 3 10 3 10 3 10 4 10 4 10 4 4 3 2 1 0 1 2 3 4 4 3 2 1 0 1 2 3 4 4 3 2 1 0 1 2 3 4 ( . )1 ( . )2 ( . )3 (d) (e) (f) Figure 3.6: PDF of ensemble-averaged SGS forces in three directions using different models at Lδ =4 ((a), (b), (c)), and Lδ =8 ((d), (e), (f)). To employ the second and third approaches in the determination of the fractional-order, we do a test using α = αopt based on the precise determination method (first approach), then we apply the obtained power-law relations (second approach) to get αopt and finally α = 0.5 (third approach) to show its effect on the performance of the D-FSGS model. As has been depicted in Figure 3.7, even using the rough-estimation approaches for α, the D- 52 SMG 0.35 D-SMG FSGS D FSGS ( = opt) 0.30 D FSGS (2nd Method) D FSGS ( = 0.5) 0.25 0.20 0.15 2 3 4 5 6 7 8 9 10 Figure 3.7: Model performance comparison using ensemble-averaged correlation at different filter sizes using three proposed approaches for the determination of the fractional-order. FSGS model provides a higher correlation for all the filter sizes. Moreover, this plot again verifies the existence of power-law relation between the αopt and filter size. One of the main disadvantages of nonlocal models in general and nonlocal turbulence models, in particular, is their sensitivity to the fractional-order. Surprisingly, the dynamic nonlocal model exhibits a less sensitivity to this parameter around the value 0.5, since its effect is compensated by the automatic tuning parameter, thanks to the designed dynamic procedure. Moreover, the whole closure term in the LES framework indeed originates from filtering the first-order nonlinear convective term in the N-S equations. Here, the corresponding αopt = 0.5 in fact validates that the overall scaled up and nonlocal behavior of the closure term in our dynamic model emerges of total 2αopt = 1st order. 3.4 A Posteriori Analysis To assess the practical ability of the proposed dynamic nonlocal model to capture un- steady large-scale coherent structures and verify numerical stability, we carry out the cor- responding a posteriori study here. In this section, we test the performance of the new proposed model (D-FSGS) against the static and dynamic Smagorinsky models (SMG, D- SMG) as well as the ground-truth filtered DNS results. For this purpose, pseudo-spectral N-S solver, which has been discussed in [147], was utilized for a triply periodic domain 53 10 1 10 1 10 2 10 2 10 3 E(K) 10 3 10 4 DNS DNS SMG 10 5 SMG 10 4 D-SMG D-SMG D-FSGS 10 6 D-FSGS UNS UNS 10 5 K 5/3 scaling 10 7 K 5/3 scaling 100 101 102 100 101 102 K K (a) (b) Figure 3.8: Prediction of the kinetic energy spectra in different turbulence models at (a) t t τ ≃ 2, and (b) τ ≃ 4 . L L as Ω = [0, 2π]3 . The resolution for this dataset is 5203 and large-eddy turnover times is τL ≃ 2.7. First, we set up a decaying HIT case in which the initial condition is based on the statistically stationary data sets. Thereafter, DNS and filtered-DNS results were gathered for the sake of comparisons. Second, each LES model including D-FSGS, SMG, D-SMG, and also unresolved numerical simulation (UNS), which is basically a DNS solver using the LES grid, implemented using the same pseudo-spectral solver, initiated with the filtered DNS initial condition for that specific filter size. In the mean time, the time step of the problems is chosen to have the CFL less than unity to ensure a stable time-integration. We compared the kinetic energy spectrum obtained from each model with the ground- truth ones. In this section we are just showing the results of using Lδ = 4 since in bigger filter widths, there is not enough resolution especially in the small scale sections. Figure 3.8 illustrates energy spectrum in two relatively high integration times τt ≃ 2, 4 after initiation L of decay. Higher integration times were chosen to address the numerical stability of models, which can be a concern especially in the dynamic models. It is well noted that Smagorinsky model is one of the most dissipative models and this matter is evident in both plots. We see that the D-FSGS model predicts similar dissipative levels of D-SMG for the majority of wavenumbers; however, we see a good correction in the small scale section thanks to the 54 1.6 Filtered DNS 1.4 Filtered DNS D-SMG D-SMG 1.4 D-FSGS 1.2 D-FSGS 1.2 1.0 1.0 0.8 uiui 0.8 0.6 1 2 0.6 0.4 0.4 0.2 0.2 0.0 0.0 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 t/ t/ (a) (b) Figure 3.9: Decay of the resolved turbulent kinetic energy in different turbulence models using (a) Lδ = 4 as an LES case, and (b) Lδ = 20 as a VLES case. dynamic and nonlocal nature of the model. Moreover, the UNS case as one may expect should have an accumulation of noise due to the zero turbulent viscosity. Therefore, we can expect that the D-FSGS model should be ideally placed between highly dissipative SMG and UNS energy curves, which act as lower and higher bounds, respectively. We continued the implementation of the a posteriori analysis by comparing the perfor- mance of the models in decaying homogeneous isotropic turbulence with Reλ = 240 using LES filter width Lδ = 4, 20 as representative of LES and VLES regions, respectively. Since the DNS and LES solvers have not been initiated with the same initial conditions, we are comparing the results of filtered-DNS and LES models. Figure 3.9 shows the decay of re- solved kinetic energy with time. As the plots show, the D-FSGS model provides better agreement with the filtered-DNS results in both LES and VLES cases. Interestingly, this difference is more remarkable in using bigger filter sizes and can be due to the fact that enlargement of the filter size incorporates more details to the solver and local models are not successful enough to correctly handling them due to the nature of the model. On the other hand, the D-FSGS model shows better performance by taking into account the nonlocal features and fractional derivations. As the last and most stringent test in the context of the two-point diagnostics, higher-order structure functions for the velocity increments are 55 0.40 Filtered DNS SMG 0.000 0.35 D-SMG 0.005 0.30 D-FSGS 0.010 0.25 ( r uL )2 ( r uL )3 0.015 0.20 0.020 0.15 0.025 Filtered DNS 0.10 0.030 SMG 0.05 D-SMG 0.035 D-FSGS 0.00 101 102 101 102 r/ r/ (a) (b) Figure 3.10: Two-point diagnostics, velocity structure functions in different models and comparison to the filtered-DNS using Lδ = 4 at τt ≃ 2 for (a) n=2, and (b) n=3 in Eq. L 5.25. compared with the ground-truth filtered-DNS results. The second and third order structures functions are computed based on the following relation  n ⟨(δr ūL )n ⟩ = ⟨ ūL (x + reL ) − ūL (x) ⟩, n = 2, 3. (3.25) To get the structure functions, first, we compute the velocity fields at a certain time for each model then we shift them based on uL (x + reL ), and finally, we do the filtering operation on each obtained field. The spatial shift is related to the filter size as r = Lδ . Comparing the results for the second and third-order structure functions in Figure 3.10 for τt ≃ 2, reveals L that all models are almost preserving the second-order structures, however, in the third-order structures we have more deviation for the Smagorinsky-based models. The new proposed model is more successful in the prediction of the third-order structures and its trend in the LES region. 56 CHAPTER 4 DYNAMIC NONLOCAL CLOSURE MODELING FOR SCALAR TURBULENCE 4.1 Introduction Turbulence remembers and is fundamentally nonlocal. Such a longing portrait of turbu- lence originates from the delineation of coherent structures/motions, being spatially spotty, giving rise to interestingly anomalous spatio-temporal fluctuating signals [40].The statisti- cal anomalies in such stochastic fields emerge as: sharp peaks, heavy-skirts of power-law form, long-range correlations, and skewed distributions, which scientifically manifest the non-Markovian/non-Fickian nature of turbulence at small scales. Such physical-statistical evidence highlights that ‘nonlocal features’ and ‘global inertial interactions’ cannot be ruled out in turbulence physics. On a whole different (computational) level and in addition to the aforementioned picture, the very act of filtering the Navier-Stokes and the energy/scalar equations in the large eddy simulations (LES) would make the existing hidden nonlocality in the subgrid dynamics even more pronounced, to which it induces an immiscibly mixed physical-computational nonlocal character. This urges the development of new LES modeling paradigms in addition to novel statistical measures that can meticulously extract, pin-down, and highlight the nonlocal character of turbulence (even in the most canonical flows) and their absence in the common/classic turbulence modeling practice. One of the oldest and most conventional local SGS modeling is known as the Prandtl- Smagorinsky model model (PSM) and was initially conceptualized in [138]. Despite being a significant step forward in LES studies, this eddy viscosity-based model suffers from a low correlation ratio, lack of back-scatter prediction, and flow-dependent features. Backward transfer of kinetic energy (back-scattering) from small scales to large scales is an innate feature of the turbulent flows that does reflect in the DNS and experimental studies. Nev- ertheless, most of the proposed turbulence models only predict the cascade of the energy from large to small scales. Another import weakness for the static PSM model is that there 57 is not a single universal constant for the representation of different turbulent fields such as shear flows, rotating flows, or wall-bounded flows. As a remedy for the last two drawbacks, Germano et al. [38] proposed a new model, which is called the dynamic Prandtl-Smagorinsky model (DPSM). They designed a dynamic procedure for the computation of the model con- stants as the calculation progresses. This procedure is based on the local calculation of the eddy viscosity coefficient by sampling the smallest resolved scales and using the obtained in- formation in modeling the subgrid scales. Afterward, different dynamic models were designed and proposed based on the same concept [141, 140, 142]. Simulations of turbulent flows using DNS and experimental studies demonstrate that turbulence is intrinsically nonlocal [1, 2]. Nonlocality of turbulence emerges as the sharp peaks, heavy-skirts, and skewed probability density function (PDF) in statistics of the veloc- ity/scalar increments. Nevertheless, most of the turbulence models are based on Boussinesq’s turbulent viscosity concept, which assumes the turbulent stress tensor to be proportional to the local mean velocity gradient. Bradshaw [3] discussed that this assumption is not neces- sarily true everywhere and it does fail in some scenarios like curved surfaces. However, local models dominantly were utilized due to their easier implementations and absence of handy and feasible nonlocal models. Introducing the nonlocality concept to the mathematical models can be done in different ways for a variety of applications. The most applicable and convenient one is based on using generalized-order derivatives. In the recent years, there have been remarkable studies in utilizing generalized-order derivatives including anomalous rheology [132], damage modeling [8] and many more that can be found in [134, 137, 11, 13]. In turbulence modeling, there is also an emerging interest in recent years in the developments of nonlocal models. Egolf and Hutter [25, 28] introduced nonlocal models based on the Reynolds-averaged Navier–Stokes (RANS) coarse-grained technique. Some additional works in this area can be found in [27, 23, 24]. In the LES turbulence modeling, a fractional Laplacian-based model was developed for the homogeneous isotropic turbulent (HIT) flows in [29]. They utilized Lévy stable 58 distributions in kinetic level and finally derived a nonlocal model that addresses the non- Gaussian statistics of the turbulent flows. Later they extended their modeling approach through developing a tempered fractional model using a tempered Lévy stable distribution, which resulted in a promising performance in a priori and a posteriori tests [30]. A fractional eddy-viscosity-based model proposed in [34] and a priori tests were performed for the HIT and turbulent channel flow as canonical test cases. On the nonlocal turbulence modeling for the wall-bounded turbulent flows [33] proposed a class of turbulence model based on the fractional partial differential equations with stochastic loads. Additionally, one can consult with some preliminary related studies in [35, 36, 37]. Modeling of the residual flux for the LES of conserved passive scalars (such as temper- ature) transported in turbulent flow medium has also been an important direction in the computational turbulence research, natural, and engineering applications. Due to the advec- tive coupling with the turbulent velocity field, the fluctuations in passive scalar concentration field are known to be more intermittent and non-Gaussian compared to the velocity field [150, 151, 152]. This behavior results in considerably stronger deviations of passive scalar statistical temporal records (turbulent intensity and dissipation) from their mean values (in a stationary turbulent regime) when we compare them to their counterparts in the turbulent velocity fields [see e.g., 153, 154, 155, 156]. As a result, the residual scalar flux emerging in the governing equation for the LES of a turbulent passive scalar is naturally carrying a com- plex dynamics. Using a well-resolved filtered DNS dataset for the residual scalar flux, it has been shown that the subgird-scale dynamics has a “statistically nonlocal” nature that cannot be ruled out through conventional functional means of modeling such as those relying on the local eddy-diffusivity assumption [see 32, Sec. 3]. Therefore, using a detailed mathematical derivation starting Kinetic theory, [32] obtained a fractional-order SGS model for scalar flux that successfully reproduced the nonlocal and non-Gaussian behavior of the residual scalar flux. Similar to Smagorisky model for passive scalars, the nature of their model is static in terms of the model coefficient and for its identification they relied on an a priori regression 59 approach with respect to the FDNS data. In practice, this a priori regression step maybe found cumbersome and a dynamic procedure seems to be a proper modification to improve the generality of this model. In this work, we develop a new dynamic nonlocal passive scalar (DNPS) subgrid scale closure model. Both a priori and a posteriori assessments were performed on the model in order to test its performance. The results were compared with the conventional static and dynamic Prandtl-Smagorinsky models. As well as improving the performance of the static nonlocal passive scalar model (NPS) [32], the new model incorporates back-scatter prediction and does not require prior knowledge for the determination of the model constant. The remaining portions of this research are arranged as follows: in section 4.2, we provide the governing equations and derivation of the proposed model. In section 4.3 we elaborate on the importance of the fractional order and its identification method. Section 4.4 is dedicated to the a priori assessments and comparing the performance of the new model with the conventional PSM-based ones. Continuing the model performance tests, in section 4.5, we analyze the model performance and its numerical stability inside a solver. 4.2 Development of the new model (DNPS) Focusing on the incompressible flow regime and the transport of a conserved passive scalar (such as temperature field) in that flow medium, Navier-Stokes (NS) and Advection- Diffusion (AD) equations are the set of governing equations that constitute the dynamics [125]. In the Large-Eddy Simulation (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., 126] as ∂ ū 1 + ū · ∇ū = − ∇p̄ + ν ∆ū − ∇ · τ R ; ∇ · ū = 0, (4.1) ∂t ρ ∂ Φ̄ + ū · ∇Φ̄ = D ∆Φ̄ − ∇ · q R . (4.2) ∂t In these equations, u = (u1 , u2 , u3 ), p, and Φ are the velocity, pressure, and scalar con- centration fields, respectively. In (4.1), ρ denotes the fluid density, while ν represents the 60 viscosity of fluid, and in (4.2), D is the diffusivity of the passive scalar field. Moreover, fil- tering yields sources of closure in the LES governing equations as the divergence of residual stress, τ R = uu − ūū, and residual flux, q R = uΦ − ūΦ̄. 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 [127, 126]. The Reynolds decomposition for a general field such as scalar concentration, Φ = ⟨Φ⟩+ϕ, where ⟨Φ⟩ is the ensemble-averaged part of Φ, and ϕ denotes its fluctuating part [125]. 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) [157, 147]. As a result, the filtered AD equation (4.2) is rewritten as: ∂ ϕ̄ + ū · ∇ϕ̄ = −ū2 + D ∆ϕ̄ − ∇ · q R . (4.3) ∂t In (4.3), the residual scalar flux may be restated as: q R = uϕ − ūϕ̄. The goal of our study is to focus on developing “predictive” and “automated” approaches for modeling q R in a dynamic setting. Nonlocal modeling for the residual scalar flux An important element of a predictive modeling is the capability of the model to reproduce the main characteristics of the quantity that aimed to be predicted. In the LES of turbulent transport, nonlocality of the SGS dynamics requires a careful attention so that the prediction of the important statistical quantities such as resolved scalar variance, ⟨ϕ̄2 ⟩, would be realistic over the course of a long-term simulation. In a comprehensive study by [32], using a rich filtered DNS (FDNS) data set, it has been illustrated that Π = −⟨q R · ∇ϕ̄⟩, (which is the SGS contribution to the time evolution of resolved-scale scalar variance) has a strong nonlocal behavior towards the larger filter sizes in way that: (i) the normalized probability distribution function of Π exhibits heavier tails, (ii) the normalized two-point correlation function, −⟨q R (x) · ∇ϕ̄(x + r)⟩/Π, yields higher values at a fixed shift value, r, especially within the inertial-convective subrange. Moreover, they showed that classical eddy-diffusivity 61 modeling (EDM) for the SGS scalar flux fails to address this nonlocal behaviour regardless of the filter size. As a result, [32] developed a fractional order SGS model for the residual scalar flux that successfully reproduced the nonlocal behavior they observed in the filtered DNS data set. Their mathematical modeling was originated from investigating the source of LES closure at the kinetic level from the filtered Boltzmann transport equation (FBTE), ∂ḡ ḡ − g eq (L) + v · ∇ ḡ = − . (4.4) ∂t τg FBTE (4.4) governs the time-evolution of distribution function for a single passive scalar particle , g = g(t, x, v) at time t with particle’s spatial location x and velocity v. In (4.4), the well-know Bhatnagar–Gross–Krook (BGK) kinetic model for the collision of two particles is utilized and it is characterized with a single parameter called relaxation time (τg ). Moreover the BGK model assumes a local equilibrium distribution function g eq (L) (known as the Maxwell distribution) for the two-particle collision, which is a normally distributed function of L that is parameterized by the locally conserved quantities. In the FBTE, the source of closure stems from the fact that the filtering operator does not commute with the ¯ Given the fact the g eq (L) is a filtered collision operator, which results in g eq (L) ̸= g eq (L). normal distribution and has a multi-exponential nature, [32] modeled its behavior with an α-stable Lévy distribution as a function of L̄, which closes the FBTE. Using the ensemble- averaging over the 3-D space, they derived the continuum-level filtered AD equation with the modeled residual scalar flux as: q R = −Dα R(−∆)α−1/2 ϕ̄, α ∈ (0, 1]. (4.5) In (4.5), R(−∆)α−1/2 (·) represents the fractional order gradient operator through the Riesz transform (see the Appendix A), and Dα is a positive real-valued model coefficient. For more details on the derivation, the interested readers are referred to [32, Sec. 4 and Appendix B]. In particular, obtaining Dα requires a priori model identifications such as regression utilizing the true values of q R from filtered DNS data. This procedure, is inherently making 62 the modeling procedure impractical and more complicated. In order to address this issue and elevate the modeling framework to an automated level, we develop a dynamic procedure based on the fractional order SGS scalar flux given in (4.5). Dynamic nonlocal modeling for the residual scalar flux Following the concept introduced by [38], the fluxes at test-level filter (subtest-scale fluxes) can be written as QR b̄ α−1/2 ϕ. i = −Dα Ri (−∆) (4.6) c assigned for showing filtering at the test-level. The ratio between the test-level where, (·) and grid-level filter sizes are usually chosen equals to two. Resolved scales and subgrid-scales are being related through the Germano identity. Gi = QR cR i − qi , (4.7) replacing the previously achieved relation in (4.7) gives V Gi = b̄ + D R (−∆)α−1/2 ϕ̄. −Dα Ri (−∆)α−1/2 ϕ (4.8) α i By assuming scale-invariance condition for the model constant and considering the known quantities for the definition of this identity, we can simplify the Germano identity as V  V  b̄ ϕ̄ ūi − ϕ ub̄i = −Dα Ri (−∆) α−1/2 b̄ ϕ − Ri (−∆) α−1/2 ϕ̄ . (4.9) It is a common assumption that the filtering procedure commutes with the integer-order derivatives and this assumption is basically one of the main steps of deriving the filtered equations in LES methods. However, we show that this is not a valid premise when one deals with the fractional operators. It can be seen in (4.9) that the filtering process does not commute with the fractional operators. The difference between the fractional derivative of filtered scalar field and the filtered fractional derivative at the test-level filter is not a negligible amount. Actually, we are using this difference to estimate the behavior in the modeled section. Several methods have been suggested in this section for the calculation 63 of the model constants. Since the scope is having a scalar model constant at the end, we need to first contract the tensorial quantities. One approach that was utilized originally by Germano et al. [38] for the dynamic Smagorisnky method is based on contracting using the filtered strain rate tensor. The other usual method which is more robust and fruitful is suggested by [146]. The second approach is based on the least square method (LSM) and assumes a squared error as V 2 b̄ u ϕ̄ ūi − ϕ b̄i + Dα Mi = e. (4.10) in which Mi is defined as V Mi = Ri (−∆)α−1/2 ϕ b̄ − Ri (−∆)α−1/2 ϕ̄. (4.11) Now, we can put the derivative of squared error equals zero to find the model constant. Also, 2 we double-check for the correctness of ∂ e2 > 0 in order to get the minimized values for the ∂C error. Therefore, the scalar model constant would be calculated as V ⟨(ϕ̄ ūi − ϕb̄ u b̄i ) Mi ⟩ Dα = − . (4.12) ⟨Mi Mi ⟩ We also added averaging operators to avoid numerical instability due to the negative eddy diffusivity. This approach of averaging over the directions of statistical homogeneity is pro- posed by [38]. This paper focuses on the development of a dynamic, functional nonlocal model. The same concept can also be applied to build a hybrid model (functional + struc- tural), but that is not the focus of this paper. In Appendix A, we have provided the necessary mathematical concepts and definitions that being used in the model development procedure. 4.3 Data-driven Identification of Optimum Fractional Order, αopt We introduce a data-driven approach to determine the optimum fractional order in the derived DNPS closure model. The new nonlocal proposed model, like every other nonlocal one, has a tuning knob which is called fractional order (α). According to (4.5), this value can take values between zero and one, and setting different fractional orders would result in 64 different fractional derivatives. Therefore, an optimization algorithm is needed to determine this value before deployment of the model in a priori and the a posteriori tests. Several approaches and criteria could be imagined for this optimization procedure such as (1) se- lection of α when gaining maximum correlations for the average of the scalar fluxes and ground-truth DNS results in three directions or (2) when there is a maximum correlation for the divergence of fluxes from the model with the ground-truth results, and (3) when we have a minimum mean squared error comparing the divergence of the fluxes (closure terms). Preliminary analysis showed that the second and third approaches are providing very close results and also we have a better prediction of back-scatter phenomenon in case we use the second (or third) approach rather than the first one. Therefore, the second approach is selected and utilized in this research. We construct our database based on 10 sample snapshots of DNS simulations using the pseudo-spectral parallel code elaborated in [147]. Using the mentioned framework we generate a stationary HIT flow with 5203 resolution. The computational domain is a cube as Ω = [0, 2π]3 , and the Taylor-scale Reynolds (Reλ ) and Schmidt (Sc = ν/D) numbers are 240 and 1, respectively. The 10 snapshots are sampled out from the DNS simulation after being sure of reaching statistically stationary condition over 10 large-eddy turnover times [147, 32]. Thereafter, the discussed optimization procedure is implemented on each snapshot and finally, the ensemble-averaged quantities are reported and utilized for the discussions and analysis. As it has been reported in [158, 30, 29, 32], changing the LES filter size would change the optimum fractional order. Moreover, Lδ = ∆∗ ∆DN S , where ∆∗ is the ratio between the LES and DNS grids, and the DNS grid size is defined as ∆DN S = 2π N , N = 520 is the DNS resolution. In this study, we utilize three different filter sizes of Lδ = 4, 10, 20. The first filter size (Lδ = 4) would represents a conventional LES filter size, the second filter size (Lδ = 10) is called the coarse LES, and finally the third filter size (Lδ = 20) is trying to test the model performance in the very large-eddy simulation (VLES) studies. Our criteria for this classification are based on the amount of filtered scalar variance that is being modeled 65 0.7 Lδ = 2 Lδ = 10 y = 1.0002Lδ −0.131 (a) 0.90 (b) Lδ = 4 Lδ = 20 0.6 0.85 0.5 0.80 ρ αopt 0.4 0.75 0.3 0.70 0.65 0.2 0.2 0.4 0.6 0.8 1.0 2 4 10 20 26 α Lδ Figure 4.1: Obtaining αopt at different filter sizes in one of the snapshots of data and (b) Power-law relation between the αopt and Lδ based on the ensemble-averaged results. in the turbulence modeling procedure. The amount of the modeled filtered scalar variance for these three filter sizes are 4%, 9% and 18%, respectively. Figure 4.1(a) shows the process of finding optimum fractional order based on the highest correlation for the forces in different filter sizes in one of the sample snapshots of data. As it is clearly seen, the maximum values of the plots are moving toward zero by increasing the filter size. This reverse relation has also been reported in [158, 29] and is due to the fact that by increasing the filter size, more nonlocality is incorporated (see section 4.2). Increasing the filter size, usually results in lower correlation coefficients (ρ); however, one can see that by going from Lδ = 10 to Lδ = 20 there is not a significant reduction in the performance. This interesting feature will be elaborated more in the upcoming sections. Finding the exact values of the optimum fractional orders based on the ensemble averaging over the 10 snapshots of filtered DNS data, demonstrates that there is a power-law relation between the filter size Lδ and αopt . As it has depicted in Figure 4.1(b), the trend-line is clearly showing a power-law behavior in which the R2 values, the coefficient of determination, for these regression procedures are above 0.996. 4.4 A Priori Tests We perform a comprehensive a priori test on the newly developed DNPS model and compare the results with the results of the conventional eddy-diffusivity-based model includ- 66 F DN S DP SM 100 F DN S DP SM (a) P SM DN P S (b) P SM DN P S 10−1 10−2 PDF 10−3 PDF 10−5 10−4 −40 −20 0 20 40 −10 0 10 h∇.qi h∇.qi F DN S DP SM 100 P SM DN P S (c) 10−1 PDF 10−2 10−3 −5 0 5 h∇.qi Figure 4.2: Normalized PDF of ensemble-averaged SGS forces using different models and different filter sizes, (a) Lδ = 4 for LES, (b) Lδ = 10 for coarse LES, and (c) Lδ = 20 for VLES. ing static Prandtl-Smagorinsky (PSM) and dynamic Prandtl-Smagorinsky (DPSM) models. The ground-truth results are achieved using the filtered DNS data (FDNS). We have depicted the ensemble-averaged forces related to the closure term in Figure 4.2 using previously mentioned filter sizes to cover all scenarios regarding the characteristic LES filter size. The first graph (a) which belongs to the conventional LES scenario, Lδ = 4, shows that the new DNPS model provides better prediction in capturing the peak and both right and left wings of the PDF. Also, the one-point correlation coefficient between the DNPS and FDNS results is 0.48 which is remarkably higher than PSM and DPSM models with 0.19, 0.40 correlation coefficient, respectively. Increasing the filter size to Lδ = 10 in second plot 67 F DN S DP SM F DN S DP SM 100 (a) P SM DN P S 100 (b) P SM DN P S 10−2 PDF PDF 10−2 10−4 10−4 10−6 −50 0 50 100 −20 0 20 −hq.Ḡi −hq.Ḡi F DN S DP SM (c) P SM DN P S 100 10−1 PDF 10−2 10−3 10−4 −10 −5 0 5 10 −hq.Ḡi Figure 4.3: Normalized PDF of ensemble-averaged SGS dissipation using different models and different filter sizes, (a) Lδ = 4 for LES, (b) Lδ = 10 for coarse LES, and (c) Lδ = 20 for VLES. (b) clearly shows that conventional PSM-based models are deviating from the ground-truth values. In the third plot (c), we have the results for the Lδ = 20 which is categorized as a VLES case. In this level and previous filter size, the DPSM and PSM models perform closely. Therefore, using dynamic model does not necessarily increase the performance. However, we observe the DNPS model still maintains the performance fairly high, according to the predictions of the peaks and tails of the distributions and correlation coefficients. During the next test of the a priori section, we assess the models’ ability to predict the back-scattering phenomenon. Computing the SGS dissipation of the scalar variance provides a clear insight into the capability of the SGS models in reproducing the backward scattering 68 in the turbulent cascade. This quantity is defined as Π = −⟨q R · Ḡ⟩, (4.13) in which, Ḡ represents the gradient of the filtered scalar field. Static models fail to capture the negative values for dissipation. However, the left wings do exist in the PDF of the ground-truth filtered DNS data, so it is part of the true physics. Furthermore, one of the major purposes of dynamic procedures is adding this important feature to the static models. Results in Figure 4.3 illustrate the predicted dissipation using different models and filter sizes according to the above definition. As the three plots are showing, the PSM model fails to predict the left wings, back-scatter phenomena. DPSM model does predict the flow of scalar intensity from small scales to large scales, but the DNPS model has better agreement with the FDNS results. Moreover, in large filter sizes (Lδ = 10, 20) we see that the DPSM model fails in predicting the peak of the PDF. Obtained results illustrate that dynamic nonlocal modeling not only adds the capability of back-scatter prediction to the static model, but also significantly increases its performance. 4.5 A Posteriori Tests Evaluating the performance of any SGS model is ultimately targeted in an LES setting, where instead of utilizing the filtered DNS variables to construct the modeled closure terms, one can use the LES-resolved flow variables, and apply the modeled closure for solving the LES equations through time. This method of assessment is called a posteriori analysis which is coined by [159] highlighting that the turbulence model is examined after being implemented in a numerical solver. Similar to the a priori testing the reference values for comparisons are obtained from filtering the DNS-resolved flow variables. As a common practice in a posteriori analysis of LES models, the time records of turbulent intensities are compared with their counterparts obtained from filtering the DNS results. For instance, in assessment of a model for the SGS scalar flux, resolved-scale scalar variance, 12 ⟨ϕ̄2 ⟩, is the target turbulent intensity. Moreover, in more robust assessments, complex statistical quantities such as high-order 69 structure functions of resolved-scale turbulent fields are compared with the ones obtained from filtered DNS [see e.g., 155, 30, 158]. This type of examination, provides a sophisticated information on the two-point (indicating nonlocal correlations) and high-order statistical performance of the SGS model in an LES. In order to perform the large-eddy simulations on the problem setting introduced in section 4.2, we employ the open-source pseudo-spectral solver developed in [147]. We slightly modify this DNS framework to account for multiple SGS models of interest for q R , including the DNPS model we developed. This DNS framework has already been successfully utilized for the LES of decaying HIT flows with a variety of implemented models for the SGS stress tensor [30, 158]. In an LES setting, both SGS stresses and fluxes in the filtered NS and AD equations are modeled, and a specific choice of SGS model for τ R could have a dominant or mixed effects on the resolved scalar concentration field [155]. In our study, we choose to freeze these potential effects; therefore, we would be able to fully concentrate on the modeling aspects and performance for q R in (4.3). As a result, in our numerical setup, we resolve the NS equations on DNS resolution, and throughout explicit filtering after each time-step we provide the velocity field on the LES resolution which is fed to the equation (4.3). This procedure has been employed in earlier studies such as [155, 160]. According to the turbulent regime with Reλ = 240 and Sc = 1 (as utilized in the calibration of αopt ), we choose a fully-developed turbulent state for the velocity and scalar concentration from a well-resolved DNS. In order to reach to this turbulent state, the NS equations are resolved for approximately 15 τLE while an artificial forcing mechanism is enforced to the low wavenumbers (energy containing range) to maintain the turbulent kinetic energy [147]. Afterwards, a we start resolving the AD equation from an initial fluctuating concentration of ϕ0 (x) = 0 while imposing a uniform mean-gradient as described in section 4.2. By resolving the NS and AD equations for approximately 15 large-eddy turnover times, the skewness and flatness records of the fluctuating passive scalar gradient reaches to a statistically stationary state, ensuring the fully-developed turbulent scalar fluctuations are 70 (hφ̄2imodel − hφ̄2iref.)/hφ̄2iref. (hφ̄2imodel − hφ̄2iref.)/hφ̄2iref. P SM NP S (a) (b) 0.2 DP SM DN P S 0.6 0.4 0.1 0.2 0.0 0.0 −0.2 −0.1 0 1 2 3 4 5 0 1 2 3 4 5 t/τLE t/τLE Figure 4.4: Relative error (with respect to the FDNS) in the records of scalar variance predicted in the LES for (a): Lδ = 10 and (b): Lδ = 20. achieved. This procedure was successfully exercised in [147, 32, 156]. By the explicit filtering of the velocity and scalar concentration fields, we take the initial condition for our LES tests. In our tests, we implement two static SGS models (PSM and NPS) in addition to two dynamic ones (DPSM and DNPS), and employ them in time-integration of filtered AD equation (4.3) for 5 large-eddy turnover times. As mentioned in section 4.2, we utilized a stabilization process to avoid numerical instability. In this process, an averaging operator over the directions of statistical homogeneity is performed [38]. Since our solver provides spectral accuracy in space and is not dissipative enough to handle the negative eddy-viscosity in long-time integration, we need to have this step in the implementation of the dynamic models in the a posteriori assessments. Records of scalar variance Evolution of the scalar variance is an important indicator in reliable prediction of the turbulent intensity. Figure 4.4, shows the temporal records of relative error in the resolved- scale scalar variance using different SGS scalar flux models with respect to the obtained time record from filtering the DNS solution as the reference temporal record. The scalar variance errors are reported for the cases of coarse LES (Lδ = 10) and VLES (Lδ = 20) filter scale resolutions. The LES test cases as well as the reference DNS are conducted for 5 large-eddy turnover times. 71 For the large-eddy simulations on Lδ = 10 resolution (Figure 4.4a), we observe that at initial stage of the simulation t/τLE < 2, the DPSM model exhibits the lowest error com- pared to the other models in a way that on average over this time-span, the absolute value of the errors in the NPS, PSM, and DNPS models are 2.9, 2.5, and 1.4 times higher than the DPSM, respectively. However, for 2 < t/τLE , the records of relative error indicate that nonlocal models (NPS and DNPS) start to perform better compared to the PSM and DPSM models. For instance, for the time-interval of 2.5 < t/τLE < 5, a comparison between the time-averaged errors from each model shows that the absolute value of error in the PSM, DPSM, and NPS models are 25, 18, and 7 times higher than what the DNPS model yields, respectively. In the LES tests with Lδ = 20 (Figure 4.4b), we observe approximately 9-10 times larger errors for the PSM and DPSM models compared to the DNPS model, when the records of error are averaged over the first 2 large-eddy turnover times of LES. This clearly shows that unlike the PSM and DPSM models, the DNPS model is reliably capa- ble of keeping the modeling error at an acceptably low level while the LES is adjusting to the initial condition. For 2 < t/τLE < 5, the DPSM model exhibits lower errors com- pared to its non-dynamic version (approximately 25% less error); however, DPSM model’s error is approximately 3 times larger than the error recorded in the DNPS model. These numerical/statistical observations indicate that the dynamic procedure effectively improves the prediction of turbulent intensity in the LES at the long-time integration region. More- over, they prove that the nonlocal models exhibit a remarkably better performance in the long-term prediction of resolved-scale scalar variance for different filter-scale resolutions. Scalar structure functions and statistical nonlocality Investigating the nonlocal behavior of turbulent regime is a vital and ultimate task in testing testing the performance of an SGS model in evolution of the turbulent field in LES [148]. 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 72 2.0 P SM DP SM 2.0 P SM DP SM (a) DN P S DN P S (b) 1.5 1.5 Ern1.0 Ern 1.0 0.5 0.5 0.0 102 102 r/η r/η P SM P SM 6 DP SM (d) DP SM (c) DN P S 5 DN P S 5 4 4 Ern Ern 3 3 2 2 1 1 0 102 102 r/η r/η Figure 4.5: Time-averaged relative errors in the computed ⟨δr ϕ̄n L ⟩ from LES with different SGS models with respect to the FDNS as reference solution using Lδ = 20. The time- averaging is done over 4 < t/τLE < 5 for (a): n = 2, (b): n = 3, (c): n = 4, and (d): n = 5. direction where 2 ≤ n [151]. These structure functions of order n are defined as  n ⟨δr ϕ̄n L⟩ = ϕ̄L (x + reL ) − ϕ̄L (x) ; n = 2, 3, . . . , (4.14) where r is the size of spatial increment, L represents the longitudinal direction (the direction along the imposed uniform mean-gradient) [161, 155, 156], and eL specifies the unit vector along the longitudinal direction. Considering the relative error between the ⟨δr ϕ̄n L ⟩ obtained from the LES solutions using an SGS model and the ground-truth FDNS solution for ϕ̄(x) at a specific time, this error function is defined as: ⟨δr ϕ̄n ⟩FDNS − ⟨δr ϕ̄nL ⟩LES . Ern = L (4.15) ⟨δr ϕ̄n L ⟩FDNS 73 Focusing on the temporal region 4 < t/τLE ≤ 5 that the LES solution has undergone long time-integration, we select uniformly distributed samples of full-domain ϕ̄(x) in time to compute Ern up to n = 5. Since we are dealing with a statistically-stationary problem, we are allowed to take the temporal-average of these error function obtained from the sampled nth-order structure 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 field in a long-time integrated LES. Figure 4.5 illustrates this time-averaged Ern against the normalized spatial shift, r/η, for the PSM, DPSM, and the DNPS models. Fig- ures 4.5a and 4.5c are showing the relative errors is the even-order structure functions, where the DNPS model has a considerably better performance. In particular, for the second-order structure functions (Figure 4.5a) the error for the DPSM model within r/η < 200 region is considerably higher than what is observed for the DNPS model. In fact, for r/η < 100 DPSM yields errors above 10 times greater than the reasonably small and steady errors we compute for the DNPS model, and for 100 < r/η < 200 the DPSM errors are approximately 5 times larger. This important observation indicates significance of the dynamic nonlocal model in successful prediction of nonlocal statistics of ϕ̄ (within the inertial-convective sub- range) compared to its conventional counterpart (DPSM). Moreover, for 200 < r/η we can still observe that the DNPS model’s predictions for the time-averaged ⟨δr ϕ̄2L ⟩ is still 2-3 time more accurate compared to the DPSM model. Interestingly, for the time-averaged Er4 (Figure 4.15c) the similar behavior as described for the second-order structure function is observed supporting our argument on the effectiveness of the dynamic nonlocal model in pre- diction of high-order statistical nonlocality. On the other hand, the DNPS model exhibits better performance in predicting the odd-order ⟨δr ϕ̄2L ⟩ compared to the DPSM model (see 3 Figures 4.15b and 4.15d). In Figure 4.15b, it is observed that Er<200η for the DNPS model is maintained at a reasonably low and steady level; however, the DPSM model returns up to 4 times larger errors. This remarkable performance of the DNPS model is even more 5 highlighted when we look at the Er<200η in Figure 4.15d, the DPSM model returns up to 74 3 9 times larger errors compared to its nonlocal counterpart. For Er>200η (Figure 4.15b), the 5 performance of both models are the same; however, by looking at Er>200η (Figure 4.15d) the DNPS model always return lower errors. This comparison indicates that in prediction of the odd-order ⟨δr ϕ̄n L ⟩, employing the DNPS model effectively improves the accuracy especially for the high-order (3 < n) structure functions. 75 CHAPTER 5 A DYNAMIC TEMPERED NONLOCAL TURBULENCE MODEL FOR LARGE EDDY SIMULATION OF CHAOTIC FLOWS 5.1 Introduction For centuries, scientists have studied the phenomenon of turbulence, which occurs fre- quently in nature. In spite of developments in the statistical theory of turbulence, there are no useful analytical solutions to turbulent flows in geometries of engineering interest. The most straightforward approach to solving turbulent flows is a direct numerical simu- lation (DNS). Numerical solutions are used in DNS to solve discretized equations. Using a fine mesh that can resolve even the smallest motion scale, and implementing a scheme that minimizes numerical dispersion and dissipation errors, it is possible to obtain an ac- curate three-dimensional, time-dependent solution of the equations without any modeling assumptions. In this case, the only error which remains is the numerical approximation errors[125, 126, 127]. There are some important limitations to perform the DNS simula- tions. Resolving all the motion scales is possible when we have enough grid points which 9 is proportional to the 9/4 power of the Reynolds number (Re 4 ) and costs of the compu- tations scales with Re3 . More details in this regard can be found in [162]. In order to decrease the computational costs and make the simulations tractable, one approach is to use Reynolds-averaged Navier–Stokes (RANS). In this approach, we solve the equations for only the averaged quantities and all the effects of the instantaneous scales are modeled by a turbulence model. Due to the low computational costs, the RANS approach has been the backbone of industrial CFD applications in the previous decades. It is however necessary in many practical cases to understand the transient behavior of the flow and therefore the RANS approach will not be sufficient. The Large Eddy Simulation (LES) approach does not adopt the time or ensemble-averaging techniques of the RANS and instead tries to com- pute the large-scale motions directly. In the LES approach, only the small-scale motions (subgrid-scales(SGS)) are being modeled and this modeling decreases the costs associated 76 with computation significantly. Between the DNS and RANS, the LES sits somewhere in the middle in terms of computational cost and provides a higher degree of accuracy than the RANS-based approaches. Another important benefit of LES approaches over RANS is the fewer needed adjustments. The reason for this issue relies on the fact that the small scales tend to be more homogeneous and universal and hence they would be less affected by the boundary conditions than the large scales. Most of the developed LES models are based on the local closure modeling point of view and they differ in the determination of the eddy viscosity or proportionality factor[138, 38, 163, 164]. From a local point of view, the turbulent stress tensor at any point is pro- portional to the local mean velocity gradient. However, experimental and DNS studies have shown that turbulence has intrinsically nonlocal properties, as well as non-Gaussian statistics, meaning velocity increments have sharp peaks, heavy skirts, and are also skewed [34, 165, 158, 2, 136, 166, 29]. From a mathematical perspective, generalized-order derivatives are an easy way to introduce nonlocality. At the continuum level, generalized-order opera- tors lay out the fundamental heavy-tailed stochastic processes that can be properly applied to incorporating long-range interactions in a variety of contexts like damage and rheology modeling [132, 133, 134, 167], visco-elasto-plastic [9] and beam vibration analysis [131]. Fur- thermore, numerical schemes for fractional-order PDEs and integer-order PDEs can be used to harness the capabilities of generalized-order models[10, 11, 13, 168, 15, 17, 18, 19]. One of the first RANS-based nonlocal turbulence models was proposed by Egolf and Hutter [25, 28]. Their fractional model is based on the Lévy flight statistics and a generalized form of the zero-equation local Reynolds shear stress. On the nonlocal LES models, a new model is proposed based on the fractional Laplacian for the homogeneous isotropic turbulence (HIT) flow by applying the Lévy stable distributions at the kinetic level [29]. The obtained results in this research showed the potential of the fractional LES model in preserving the non- Gaussian statistics of the subgrid-scale stress motions. Akhavan-Safaei et al. [32] utilized the same concept to build a nonlocal LES model for scalar turbulence. Two-point statistics 77 were utilized in the determination of the optimum fractional parameter and comparing the probability distribution functions (PDF) of the different quantities showed the proper agree- ment for the new model’s result with the filtered DNS ones. In accordance with this study, a new nonlocal spectral transfer model, as well as a new scaling law for scalar turbulence were developed [136]. Inspired by the cascade of energy in turbulent flows and its exponen- tial decay in the dissipation rage, Samiee et al. [166] developed a tempered nonlocal model by applying tempered heavy-tailed distribution within the previously developed fractional framework. By using fractional and tempering parameters, one can characterize nonlocal structures in the turbulent inertial and dissipation ranges. An eddy viscosity model is pro- posed by Clark Di Leoni et al. [34] and tested for both the HIT and channel flow turbulent test cases using the a priory analysis. In addition to the mentioned studies, there are also other related researches including [31, 33, 36, 37]. The first dynamic fractional subgrid-scale model (DFSGS) model was developed by Seyedi and Zayernuri [158]. The new model is based on the fractional Laplacian operator, however, the model constant is computed autonomously by a novel dynamic procedure. Moreover, three data-driven approaches were proposed for the optimal determination of the fractional parameters at different filter sizes and Reynolds numbers for both the large and very large eddy simulation (VLES) of the HIT test cases. By coupling the nonlocal modeling and dynamic procedure, it is shown that the performance is increased, the model coefficient is computed automatically, and the kinetic energy backscatter is predicted correctly. The development of the dynamic nonlocal model, DFSGS, was a big step forward and highlighted the great potential of these kinds of models. In this paper, we develop a dynamic tempered nonlocal LES model in order to mimic the finite variance feature of the common patterns in nature and also the exponential decay of the kinetic energy spectra in the dissipation range. Moreover, we test the performance of the newly developed model versus the dynamic Smagorinsky (DSMG) and the DFSGS models in different test cases including HIT and multi-layer temporal shear layer in the context of the a priori and a posteriori analysis. 78 In section 5.2 we talk about the governing equations and the steps we took to develop the DTF model. In section 5.3 various a priori tests applied for both the HIT as the isotropic turbulence test case and multi-layer temporal shear layer case as the anisotropic test case. Also, the optimization algorithms and results for the determination of the fractional and temporal components are provided. Testing the real performance of the models including the numerical stability and preserving the true statistics in long integration times within a spectral numerical solver is accomplished in section 5.4. 5.2 Model Development In large eddy simulations of turbulent flows, large-scale quantities are obtained by filtering the velocity and pressure fields. The filtered governing equations for these large scales in incompressible conditions are formulated as ∂ ū 1 + ū · ∇ū = − ∇p̄ + ν ∇2 ū − ∇ · τ, (5.1) ∂t ρ ∇ · ū = 0. In which, ū is the vector of the filtered velocity , ρ denotes the constant density, p̄ represents the filtered pressure, and ν shows the kinematic viscosity. The so-called SGS stress tensor reflects the effects of small scales and shown by τij = ui uj − ūi u¯j . (5.2) This unknown term requires closure modeling. Over the last decades, several functional and structural LES models have been proposed [125, 127, 126, 169]. One of the most famous and well-accepted LES models among the community are the static and dynamic Smagorinsky models[138, 38]. The classical static Smagorinsky model (SMG) is based on an eddy viscosity- based approach and follows the Boussinesq hypothesis. SM G = −2 ν S¯ , τij (5.3) t ij 79 ∂ u¯ ∂ u¯ in which νt shows the eddy viscosity, and S¯ij = 12 ( ∂x i + ∂xj ) denotes the resolved (filtered) j i strain rate tensor. Prandtl’s Mixing Length Hypothesis is used here to construct νt . νt = Cs L2 |S̄|. (5.4) q In the above relation L reflects the effective grid scale, and |S̄| = 2S¯ij S¯ij shows the magnitude of the resolved scale strain rate tensor [127]. In different scenario, such as shear flows, wall-bounded flows, or transitional flows, this model cannot correctly predict final quantities with a single universal constant because it is flow-dependent. As a solution to this challenge, Germano et al. [38] proposed a novel procedure for estimating the model coefficient, known as the dynamic Smagorinsky model (DSMG). Researchers have utilized the same concept since the proposed procedure became a breakthrough in turbulence modeling [139, 140, 141, 142]. The model calculates the eddy-viscosity coefficient for each LES grid point locally as the calculation progresses, so no predefined inputs are needed. The scale- invariance hypothesis was used to construct this model, and the model constant was derived from the resolved section to determine the constant. Tempered Fractional Subgrid-Scale Model Our recent work included the development of a non-local tempered SGS model and we present it here briefly in order to make this work self-contained [166]. To describe the collision of particles, we use the Boltzmann kinetic model based on Boltzmann’s kinetic description of the flow. ∂f f − f eq + u · ∇f = − , (5.5) ∂t T here f = f (x, u, t) indicates the single-particle probability distribution function and shows the particles density in the phase space (x, u) at time t. Moreover, f eq represents the local equilibrium distribution function and can be defined based on the Maxwell distribution as follows ρ f eq (K) = F (K), (5.6) U3 80 |u−V |2 in which F (K) = e−K/2 , K = 2 and U represents the speed of thermal agitation. On U the left-hand side of Eq. 5.5, one can see the non-reacting particle streaming, while on the right-hand side, one can see the collision operator with a relaxation time T when it comes to the non-reacting particles. Eq. 5.5 can be solved analytically using a characteristic method and the distribution can be computed based on the equilibrium state [143] Z ∞ f (t, x, u) = e−s f eq (x − uτ s, u, t − τ s) ds. (5.7) 0 Macroscopic flow variables like density and fluid velocity can be obtained using moments of R R the f . Then, one can write ρ = f (t, x, u)du, and ρV (t, x) = uf (t, x, u)du and get the required fields. The filtered Boltzmann equation would be obtained by incorporating the filtering proce- dure into Eq. 5.5. ∂ f¯ f¯ − f eq (∆) + u · ∇f¯ = − . (5.8) ∂t T Due to the highly nonlinear nature of the collision term, the filtering kernel is not able to commute with it[144]. Therefore, there is a closure problem and it can be addressed as |u−V |2 K= since U2 f eq (K) ̸= f eq (K). (5.9) Recently, a new approach based on the tempered Lévy stable distributions is suggested in [166] to model the f eq (K) − f eq (K), i.e., f eq (K) − f eq (K) ≃ f M odel (K̄) = Dβ,λ f β,λ (K̄), (5.10) where f βλ (K̄) = ρ3 F βλ (K), where F βλ (K) represents a tempered Lévy β-stable distribu- U tion. Moreover, Dβ is called the model constant and will be calculated automatically in the new dynamic tempered fractional model (DTF). The subgrid-scale forces in this new model are obtained as [166]: (∇.τ )i = να,λ (−∆α, λ )ūi , , α ∈ (0, 1/2) ∪ (1/2, 1), λ ≥ 0. (5.11) 81 In Eq. 5.11 the resolved velocity ūi : Rd × (0, T ] → R, where d = 3 is the dimension of physical domain and T represents the simulation time, in addition, (−∆)α,λ (·) shows the tempered fractional Laplacian which can be defined as a singular integral operator given by Z ūi (x) − ūi (y) (∆ + λ)α ū i = cd P.V. dy, (5.12) |x − y|d+2α exp(λ|x − y|) Rd −Γ(1/2) 1 in which cd = d/2 . In Appendix C, we have provided the preliminary mathe- 2π Γ(−2α) cos(πα) matical concepts and definitions in this regard. It should be mentioned that λ = 0 is the case that we don’t have a tempering term and would result in the same model that we developed recently [158]. As the detailed derivations in [166] demonstrates, the SGS force term in this model is achieved using the summation of a fractional Laplacian and another term which includes the tempering parameter, λ, implemented on the filtered velocity and multiplied to a model constant. (−∆α, λ ) = ϕ10 (α)(−(−∆)α ) + ϕ11 (α, λ)(∆ + λ)α , α ∈ (0, 1/2) ∪ (1/2, 1), λ ≥ 0. (5.13) In the above relation ϕ10 (α) and the ϕ11 (α, λ) are defined based on the fractional and tem- pering coefficients and Gamma function as follows 1 ϕ10 (α) = (Γ(2α + 1) − Γ(2α)) (5.14) 2α + 3 2α + λ ϕ11 (α, λ) = (Γ(2α − 1)). 2α + 3 Calculation of the SGS stress components is easily possible using the α−Riesz potential as described in [158, 29, 166, 32, 165]. Dynamic Tempered Fractional Model (DTF) As we saw in the previous section, the final SGS force term in the new model can be simply written as (∇.τ )i = να,λ (−∆α, λ )ūi . Now, we are developing a unique dynamic procedure to automatically calculate the model coefficient. To this end, we apply the second filtering procedure which is also called test-level filter and we show it by (·). c Doing so provides the 82 divergence of the stresses in the test-level (subtest-scale) as (∇.T )i = να,λ (−∆α, λ )ub̄i . (5.15) This paper uses a common ratio and doubles the size of the test-level filter compared to the grid-level filter. Inspired by the identity introduced by Germano et al. [38] which connects the scales of the grid-level to the test-level, we show that Gij = ūd i ūj − ub̄i ūbj = Tij − τc ij , (5.16) in the divergence form the above relation can be formulated as  ∇. ūdi ūj − ub̄i ūbj = (∇.G)i . (5.17) By extrapolating and parametrizing the model constant at the smallest resolved scales, the Germano identity determines the model constant for the subgrid-scale part. To construct the Germano identity in the divergence form using the introduced tempered fractional Laplacian relations in the previous section, we obtain  V  (∇.T )i − (∇.bτ )i = να,λ (−∆)α,λ Vb̄i − (−∆)α,λ V̄i . (5.18) In the derivation of the above relation, we used the same assumption that has been used in several dynamic-based LES approaches[38, 158, 165] which is the invariant of the model constant. Therefore the filtering procedures would commute and we can simplify the relation. The next step is finding the να,λ value in Eq. 5.18. However, the mentioned equation is a tensorial relation and the model constant is not a scalar. To reach a constant scalar value, one common approach in this section is Lilly’s method [146] which contracts the equation by the right-hand side itself. We use this approach in this project based on the least-squares method (LSM) and using the previous relation for the Germano identity in divergence form, we attain (∇.G)i − να,λ Mi )2 = e. (5.19) 83 In Eq. 5.19 the e indicates the squared error and Mi shows the nonlocal term which reads  V  α,λ b̄ Mi = (−∆) Vi − (−∆) V̄i . α,λ (5.20) One important point worth mentioning is that in integer-order and common dynamic models, we assume that the filtering procedures are commutable with the integer-order derivatives. However, as we illustrated previously [158] (see Fig 1) the filtering operators are not com- muting with the fractional(general) order derivatives. This is the reason that the right-hand side of the Eq. 5.20 is not zero. In the dynamic Smagorinsky (DSMG) model which is based on the integer-order derivatives, the right-hand side of the corresponding equation would be zero if there was no strain rate tensor multiplied by the first-order derivative term. In the tempered fractional Model (DTF) we don’t have the strain rate tensor but the right-hand side is still holding due to the mentioned important difference in the integer and fractional-order derivatives. Now one can calculate the unknown model constant by putting the derivative of the error equal to zero in Eq. 5.19. Also, the scope of error minimization can be determined by considering and testing ∂2e > 0 for the mentioned relation. Lastly, one can define the ∂να,λ 2 scalar model constant as follows ⟨(∇.G)i Mi ⟩ να,λ = . (5.21) ⟨Mi Mi ⟩ In the above relation, we also added the averaging operator over the directions of statistical homogeneity as suggested by Germano et al [38] to prevent the numerical stability which may be cased in long time interactions due to the negative eddy viscosity values. 5.3 A Priori Tests In this section, we compare a priori performance of the DTM, DFSGS, and DSMG with respect to the filtered DNS data as the ground-truth results. To this end, first, we need to gather a high-resolution and reliable DNS dataset and then do the statistical analysis. We used our in-house pseudo-spectral parallel code elaborated in [147] to generate a stationary forced homogeneous isotropic turbulence (HIT) flow in a cube with 5203 resolution for a 84 periodic computation domain as Ω = [0, 2π]3 . The time integration is performed using a fourth-order Runge-Kutta (RK4) scheme while the CF L < 1. The targeted Taylor Reynolds number was 240 ( Reλ = 240). We stored 10 snapshots of the DNS data at a fully turbulent and statistically stationary state in which kmax η ≈ 1.5 over the span of 10 large-eddy turnover times. Additionally, Lδ = ∆∗ /∆DN S , where ∆∗ is the LES grid size, and the DNS grid size is defined as ∆DN S = 2π N , N = 520 is the DNS resolution. In this paper, we have tried to target the larger filter sizes and do the analysis in the coarse LES and very large eddy simulation (VLES) rather than the conventional large eddy simulation LES grid sizes. Therefore both filter sizes of Lδ = 10, 20 are tested which represents the coarse and very large eddy simulation cases based on the amount of the modeled turbulent kinetic energy corresponds to each filter size which has reported in [158]. In this section, we test the performance of the new model and compare it with different LES models in terms of the prediction of the probability density function (PDF) of the SGS forces and the dissipation of turbulent kinetic energy in two separate scenarios. The first test case is the HIT and the second one is the multi-layer temporal shear layer case. In the former test which is based on a canonical flow condition, we test the performance of the models for the fully isotropic turbulent flow situation using the ensemble-averaged quantities. In the latter, we analyze the performance of the model when we have significant anisotropy in the turbulent flow. Homogeneous Isotropic Turbulence (HIT) Flow As the first test case, we use the gathered 3D high-resolution DNS datasets to do the a priori analysis. However, first, we need to define the optimum values and sensitivity analysis for the fractional and tempered input parameters for the new model. Since we have two input parameters, we use a grid search approach for the α, λ and find the best pairs of values for the model inputs when the average correlations for the three SGS forces have the highest value with the corresponding ground-truth ones. In another word, we select our two input fractional and tempering parameters when we have the highest correlation coefficients 85 0.120 0.155 0.115 0.150 0.110 ρ 0.145 ρ 0.105 0.100 0.140 0.095 0.135 0.090 0.130 0.085 1.0 1.0 0.8 0.8 0.6 0.6 0.0 0.0 0.2 0.4 0.2 0.4 0.4 0.2 λ 0.4 0.2 λ 0.6 0.6 α 0.8 0.0 α 0.8 0.0 1.0 1.0 (a) (b) Figure 5.1: Optimization process for the α, λ in the HIT test case based on the grid search approach for Lδ =10 (a) , and Lδ =20 (b). for the average of the SGS forces and the ground-truth ones. The reason for choosing this optimization process lies in the fact that we want to cover all three directions of the SGS forces in the grid search process and not to be biased to a certain direction. Moreover, the initial results showed better performance using this approach rather than the processes introduced in [158, 165]. Samiee et al. [166] used high-order structure functions in one of their optimization processes, however, in the new model, the model constant is calculated dynamically and there is no need to implement another level of optimization. In Fig. 5.1 we see the results of the grid search based on the average of the three SGS force correlations for the filter size of 10 and 20 as the coarse LES and VLES cases. As can be seen, there is not necessarily a single pair of optimum results for this process. The optimum results can be gained using different values for the α and λ. In fact, the effect of each parameter can be compensated by the other one and an imaginary line (or surface) in the 3D results would provide the best correlations. The optimization process provides a proper intuition regarding the selection of the frac- tional and temporal values. Now we use the obtained values in the previous section and do 86 100 F DN S DSM G 100 F DN S DSM G 100 F DN S DSM G DF SGS DT F DF SGS DT F DF SGS DT F 10−2 10−2 10−2 PDF PDF PDF 10−4 10−4 10−4 −20 −10 0 10 20 −20 −10 0 10 20 −20 −10 0 10 20 h(∇.τ )1i h(∇.τ )2i h(∇.τ )3i (a) (b) (c) F DN S DSM G F DN S DSM G F DN S DSM G 100 DF SGS DT F 100 DF SGS DT F 100 DF SGS DT F 10−1 10−1 10−1 PDF 10−2 PDF 10−2 PDF 10−2 10−3 10−3 10−3 10−4 10−4 10−4 −5 0 5 −5 0 5 −5 0 5 h(∇.τ )1i h(∇.τ )2i h(∇.τ )3i (d) (e) (f) Figure 5.2: PDF of ensemble-averaged SGS forces for the HIT test case in three directions using different models at Lδ =10 ((a), (b), (c)), and Lδ =20 ((d), (e), (f)). the a priori tests. In Fig 5.2 we see the PDF prediction of ensemble-averaged SGS forces in three directions for different models and their comparisons with the ground-truth filtered DNS (FDNS) ones. The first row illustrates the PDFs using Lδ =10 which is a coarse LES grid size and the second row shows the VLES predictions using Lδ =20. As can be seen, the new model has improved the performance of the previous non-tempered predictions in almost all directions and filter sizes. This improvement is more noticeable in capturing the long tails distributions. In a comparison of the local models and the nonlocal ones, we see that the DSMG model as the representative of the local LES models, doesn’t capture the heavy tails distributions which are initiated from the nonlocal nature of the turbulence. Fig. 5.3 depicts the prediction of the ensemble-averaged SGS dissipation of kinetic energy in different models at two filter sizes of Lδ = 10, 20. Since all three compared models are dynamic, we have the prediction of the back-scattering phenomena by producing a negative dissipation. As discussed previously in the optimization process for the determination of the fractional and tempered parameters, we optimized the values based on achieving the highest SGS forces. However, we see a relatively proper prediction for the dissipation distributions 87 F DN S DSM G 101 F DN S DSM G 100 DF SGS DT F DF SGS DT F 100 10−1 10−2 PDF PDF 10−2 10−4 10−3 −5 0 5 10 −2 0 2 4 6 8 −hS¯ij τij i −hS¯ij τij i (a) (b) Figure 5.3: PDF of ensemble-averaged SGS dissipation of kinetic energy using different models for (a) Lδ = 10, and (b) Lδ = 20. and a slightly higher improvement than the non-tempered DFSGS model for the new model. Multi-layer Temporal Shear Layer As the next test case in the context of the a priori tests, we want to test the performance of the developed model along with the other dynamic models when we have significant anisotropy in the flow. To this end, we create a multi-layer temporal shear layer (or temporal mixing layer) case based on the setup discussed in detail in [170]. We initiate a turbulence flow from the linear stability theory [171]. Therefore, we introduce the nondimensional initial streamwise velocity profile u0 as ∂ψ u0 (x, y) = u0¯(y) + cnoise , (5.22) ∂y in which −y 2 L2 δ02   ψ(x, y) = e cos(8πx) + cos(20πx) . (5.23) In the above relation, δ0 shows the initial vortex thickness and L is the size of the domain which is 2π. We used the hyperbolic tangent velocity profile definition that is suggested in [172] for the mean initial profile which reads 2yL  u0¯(y) = tanh . (5.24) δ0 88 By choosing the noise factor c0 = 0.001, we make sure to maintain a small percentage of mean velocity perturbations and not to dissatisfy the Taylor hypothesis. In order to have periodicity in the normal direction as well as the other directions, an even number of shear layers was applied in this project, namely 10. Having 10 temporal shear layer accelerates the mixing and decrease the needed computational time to reach chaotic turbulent flow condition significantly. The framework here makes use of periodic boundary conditions in all three directions: streamwise, spanwise, and normal. Thus, direct numerical simulations become more feasible due to the smaller computational domain. In the next step, we add the ensemble-averaged snapshots of the HIT case over the span of 10 eddy turnover times in order to add the randomness along the dominant shear layers. This average 3D field is added by a 1 percent ratio and after multiplying an exponentially decaying function that almost kills the magnitude everywhere but in the vicinity of the temporal layers. Therefore, we see the growth of the Kelvin–Helmholtz eddies in different magnitudes and sizes along the shear layers (See Fig. 5.4). The same in-house DNS solver which is elaborated previously initialized with the prepared initial condition in the decaying fashion using a fixed time step size of dt = 0.001 and checking the CFL number to be under 1 for the numerical stability purposes. The same procedure of sampling out the DNS dataset at different time steps was carried out, and we gathered the DSN snapshots of the temporal shear layer case. Thereafter, we implemented the previously discussed optimization and grid search algorithm to find the optimum values of the fractional and tempering parameters in this anisotropy test case. Fig. 5.5 represents the sensitivity analysis graphs for the decaying temporal shear layer case using characteristic filter size of Lδ =20 initiated with the mentioned initial condition at different time steps. The procedure depicted in Fig. 5.5 provides a proper notion of the input parameters of the model. Using the information provided, we compared the PDF predictions of the SGS forces at a different phase of the mixing in different LES models. To this end we utilized the largest filter size, Lδ =20, to show the new model’s capability in the VLES cases. It should be mentioned that using the smaller filter size for this case,Lδ =10, 89 (a) (b) (c) (d) Figure 5.4: Slice of the vorticity fields of the DNS temporal shear layer case at t=0.01 (a), t=0.02 (b), t=0.05 (c) and t=0.1 (d). we obtained almost the same performance for the new model and the DSMG one. However, using the bigger filter size, we basically incorporate more nonlocal interaction information 90 0.105 0.125 0.100 0.120 ρ ρ 0.095 0.115 0.090 0.110 0.105 0.085 1.0 1.0 0.8 0.8 0.6 0.6 0.0 0.0 0.2 0.4 0.2 0.4 0.4 0.2 λ 0.4 0.2 λ 0.6 0.6 α 0.8 0.0 α 0.8 0.0 1.0 1.0 (a) (b) 0.1225 0.120 0.1200 0.1175 0.115 0.1150 ρ ρ 0.110 0.1125 0.1100 0.105 0.1075 0.100 0.1050 1.0 1.0 0.8 0.8 0.6 0.6 0.0 0.0 0.2 0.4 0.2 0.4 0.4 0.2 λ 0.4 0.2 λ 0.6 0.6 α 0.8 0.0 α 0.8 0.0 1.0 1.0 (c) (d) Figure 5.5: Optimization process for the α, λ in the temporal shear layer test case based on the grid search approach for Lδ =20 at t=0.125 (a), t=0.15 (b), t=0.175 (c) and t=0.2 (d). into the models and the nonlocal models can provide better agreement in comparison to the local ones. Comparing the new tempered and the previous non-tempered model, DFSGS, indicates that tempering is helping the model to be improved in capturing the peak and tails of the distribution. In the next a priori analysis for the temporal shear layer test case, we examine the model performance in predicting the dissipation distribution and backscattering phenomena. The DFSGS and DSMG models do predict the flow of energy from small scales 91 101 F DN S DF SGS DSM G DT F 101 F DN S DF SGS DSM G DT F 101 F DN S DF SGS DSM G DT F 100 100 100 PDF 10−1 PDF 10−1 PDF 10−1 10−2 10−2 10−2 10−3 10−3 10−3 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 (∇.τ )1 (∇.τ )2 (∇.τ )3 (a) (b) (c) F DN S DSM G F DN S DSM G F DN S DSM G DF SGS DT F DF SGS DT F DF SGS DT F 101 101 10 1 100 100 100 PDF PDF PDF −1 10−1 10 10−1 10−2 10−2 10−2 −0.4 −0.2 0.0 0.2 0.4 −0.4 −0.2 0.0 0.2 0.4 −0.4 −0.2 0.0 0.2 0.4 (∇.τ )1 (∇.τ )2 (∇.τ )3 (d) (e) (f) Figure 5.6: PDF of SGS forces for the temporal shear layer test case in three directions using different models at Lδ =20, snapshot=100 ((a), (b), (c)), snapshot=150 ((d), (e), (f)). F DN S DSM G 101 F DN S DSM G 100 DF SGS DT F DF SGS DT F 100 10−1 10−2 PDF PDF 10−2 10−4 10−3 −5 0 5 10 −2 0 2 4 6 8 −hS¯ij τij i −hS¯ij τij i (a) (b) Figure 5.7: PDF of SGS dissipation of kinetic energy using different models for Lδ = 20 (a) snapshot=100, and (b) snapshot=150. to large scales thanks to their dynamic natures, but the DTF model has slightly better agreement with the FDNS results. It should be mentioned that if the optimization process has been defined based on the dissipation capturing criteria, we would have had better results in this section. However, the SGS forces are the ones that we targeted and still, we have acceptable dissipation predictions at both two shown time steps. 92 5.4 A Posteriori Tests In this section, we do the ultimate test for the newly developed model and test its performance and numerical stability in solving the LES equations through time. Piomelli et al. [159] introduced a posteriori term for implementing the model within the numerical solver, and the use of numerical solvers for implementing the model became a common practice afterward. In order to have the reference and ground-truth results for the sake of comparison, we carried out a DNS case with high resolution and the same conditions. The filtered DNS results (FDNS) are obtained each time that is targeted for the LES model comparisons. The parallel spectral solver that discussed in previous sections is utilized for the LES and DNS simulations for a triply periodic domain Ω = [0, 2π]3 . For the new model, DTF, we can claim that the model is numerically stable since the utilized solver is spectral- based and is not dissipative like the finite difference of finite volume approaches. Therefore, if any model works within this solver and tolerates the large integration times, its stability is guaranteed within the other common solvers. In this set of experiments, as a common practice, we first test the performance and numerical stability of the models in decaying homogeneous isotropic turbulence with Reλ = 240 using the LES (VLES) filter width Lδ = 20. Fig 5.8 shows the time records of the turbulent intensities in terms of the decay of the resolved kinetic energy for different models and compares them with their counterparts obtained from filtering the DNS results. The graph shows that using this large filter size which corresponds to the VLES region, the new model provides a better agreement than the non-tempered DFSGS model and the DSMG models in short and long integration times. The unresolved numerical simulation (UNS) results belong to the case in which we run a DNS solver using an LES grid. Therefore the turbulent viscosity would be equal to zero in this case and one can expect a high-level prediction of the kinetic energy. It should be mentioned that using Lδ = 10 there was not a significant improvement over the other LES models and in this test all the models could capture the true first-order statistics. In the next a posteriori test, we analyze the prediction of the kinetic energy spectrum 93 F DN S DSM G UNS DT F DF SGS 1.0 1 2 hūi ūi i 0.5 0.0 0 2 4 6 8 t (a) Figure 5.8: Decay of the resolved turbulent kinetic energy in different turbulence models using Lδ = 20. at two different times. To this end, a comparison was made between the kinetic energy spectrum of each LES model with the ground-truth filtered DNS ones as shown in Fig. 5.9. We compared the results after one and three eddy turn-over times, which is almost 2.7 for this case. It is important to choose large integration times to address the numerical stability of models, which is an issue especially when a model is dynamic. Our plot in this section only shows the results of using Lδ = 4, since bigger filter widths do not provide enough resolution, especially at small scales. As it can be noticed from the graph, all LES models are obeying the − 53 law in the inertial subrange section and the new model provides slightly higher dissipative levels in this filter size than the other models. However, the prediction of the small-scale sections in the new model is more accurate than the other ones thanks to the tempering nature of the model and its exponential decaying rates. Testing the performance of an SGS model in evolutionary simulations of the turbulent field in LES requires investigating the nonlocal behavior of turbulent regime [148]. As the last test for the model in the context of the two-point diagnostics, we compare the higher-order structure 94 F DN S DSM G 10 −1 UNS DF SGS DT F K −5/3 scaling 10−1 10−3 E(K) E(K) 10−3 10−5 10−5 10−7 F DN S UNS DSM G DT F DF SGS K −5/3 scaling 100 101 102 100 101 102 K K (a) (b) Figure 5.9: Prediction of the kinetic energy spectra after one (a) and three (b) eddy turn-over times. functions with their counterpart filtered DNS ones. This is a relatively stringent test and most of the models usually fail to preserve the true statistics in a large filter size and after a long time of integration. The equation below shows the mathematical relation to obtain this quantity based on the filtered fields and their shifted values in space. 1.6 FDNS DSMG 0.00 DFSGS DTF 1.4 0.02 1.2 0.04 1.0 0.06 ( ruL)2 0.8 0.6 ( ruL)2 0.08 0.10 0.4 0.12 0.2 0.14 FDNS DSMG 0.0 DFSGS DTF 0.16 101 102 101 102 r/ r/ (a) (b) Figure 5.10: Velocity structure functions using different models and comparison to the ground-truth results using Lδ = 20 for n=2 (a), and n=3 (b).  n ⟨(δr ūL )n ⟩ = ⟨ ūL (x + reL ) − ūL (x) ⟩, n = 2, 3, (5.25) where r is the size of spatial increment, L represents the direction, and eL specifies the unit vector along that direction. The Fig 5.10 shows that all LES models are performing almost 95 the same in the prediction of the second-order structure functions. Third-order structure functions are being better captured by the new model, DTF. One can see that the DFSGS model results are slightly better than the DSMG model and the new tempered model adds another level of improvement in the prediction of the correct third-order structure functions. 96 CHAPTER 6 SUMMARY AND FUTURE WORKS In this study, our main goal was to devise a framework that provides a data-driven approach for nonlocal modeling of anomalous transport phenomena. We started with the stochastic study of the turbulent flows in the context of the internal cylinder flows to see the anomalous and non-Gaussian behaviour of the flow. Thereafter, we developed a novel dynamic nonlocal turbulence model that accounts for the non-Gaussian statistics. In the next chapter, we construct a new dynamic model for the scalar turbulence which has the capacity to works with very large LES filter sizes. Finally, we proposed the dynamic tempered fractional model for the large and very large eddy simulation of turbulent flows and we tested them for the HIT and temporal shear layer test cases. 6.1 Concluding Remarks • In chapter 2, our study leverages the outcome of stochastic modeling and simulations 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 math- ematical representation of the stochastic incompressible Navier-Stokes equations was presented. Further, a high-fidelity stochastic CFD framework was introduced, which employs spectral/hp element method in the forward solver and later on the stochas- tic 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 characterize the most effective stochastic factor on the total variance of kinetic energy, consequently, the “ec- centric rotation” was learned to be the dominant source of stochasticity. Later on, the 97 expected solution from a very fine uni-variate PCM discretization on the dominant ran- dom 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 spot- ted 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 ve- locity fluctuations and therefore fluctuating vorticity field as the flow evolves in time. Moreover, the statistics of flow fields were quantitatively 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 vortic- ity 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 emphasize 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 sys- tems 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 98 non-Gaussianity in flow fluctuations, we sought to study the effects of 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 ef- fects, 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 generated 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. • In chapter 3, We introduced a novel dynamic nonlocal turbulence model for the isotropic turbulent flows, which can be applied for both LES and VLES purposes. The model has been developed based on the fractional Laplacian derivative of re- solved velocity field, and a unique dynamic procedure defines the model constant in the divergence form. The effect of the fractional-order, Re number, and characteris- tic filter size in LES and VLES cases are scrutinized using multiple high-fidelity and well-resolved DNS datasets belongs to forced homogeneous isotropic turbulence. Fi- nal decisions were made based on the ensemble-averaged quantities gathered from ten separate three-dimensional snapshots distributed over enough turnover times. The optimum fractional-order for each scenario is chosen when the maximum correlation is obtained for the ground-truth and predicted stresses of the model. The obtained results indicate that the relation between the filter size and fractional-order is closely obeying a power-law form for all the assessed Re numbers. Also, there is a direct rela- tionship between the Re number and the fractional-order. Considering all the attained statistical results, two other data-driven straightforward methods for the determina- tion of the fractional-order suggested and tested successfully. Our analysis included both a priori and a posteriori assessments. In the first one, we showed that there is a higher correlation between the results of the proposed method and the ground-truth 99 DNS results in comparison to other conventional LES models. Also, the capability of the model in the prediction of the back-scatter is discussed in different filter sizes. In the a posteriori assessments we tested the model performance in long-time integration in the context of a real N-S solvers, and its ability in capturing the large-scale coher- ent structures examined. Analysis were performed in a decaying isotropic turbulence scenario, and the new D-FSGS, static Smagorinsky and dynamic Smagorinsky models implemented separately. The results show that the new model is more successful in the prediction of the resolved turbulent kinetic energy in both LES and VLES stud- ies. To test the models’ performance in the long-time integration and its numerical stability, kinetic energy spectra are compared with the filtered-DNS results. Finally, two-point diagnostics were accomplished to compare different models’ performance in the context of preserving second and third-order structures. The results demonstrate that the new model behaves better in preserving the high-order structures than the conventional ones. This study shows the promising potential of using bigger filter sizes rather than the conventional LES filter sizes. This important characteristic can be uti- lized to compensate the fractional model’s high-cost demand. Moreover, accelerated evaluation methods for the fractional operators including learning-based approaches and fast solvers can be leveraged to pave the road for the real and more practical applications of nonlocal models. • In chapter 4, we developed a novel dynamic nonlocal closure model for the subgrid- scale scalar field in the context of the large and very large eddy simulation (LES, VLES). With our high-fidelity datasets pertaining to forced homogeneous isotropic turbulence, we examined the effect of fractional order and characteristic filter sizes in LES and VLES cases. During enough large-eddy turnover times, we utilized ensemble- averaged quantities from ten separate three-dimensional snapshots to make final deci- sions. When the ground-truth force and the predicted SGS force were most correlated, the optimal fractional order was selected for each scenario. We initially tested the 100 proposed model in the context of a priori assessment. Based on the results, we showed that the new DNPS model is more accurate in predicting SGS dissipation and force terms. Furthermore, the DNPS model retained a fair performance in the VLES cases, unlike other conventional models that did not return accurate prediction of LES closure at very large filter sizes. In an LES setting, we managed to examine the performance of SGS models in an a posteriori sense. We resolved the filtered AD equation for 5 large-eddy turnover times while ū contribution in the advective coupling obtained from the explicit filtering of the DNS solution. We looked at the relative error (with respect to the FDNS) in the records of the ⟨ϕ̄2 ⟩ from the LES solution using differ- ent SGS models. Tracking these records showed that the combination of nonlocal modeling and dynamic procedure (DNPS modeling) is an effective way for accurate prediction of the resolved-scale turbulent intensity (scalar variance) especially when the goal is to study the long-term behavior. Moreover, we examined the prediction of the longitudinal resolved-scale scalar structure functions, ⟨δr ϕ̄nL ⟩, for 2 ≤ n ≤ 5. Com- pared to their time-averaged FDNS solution over 4 ≤ t/τLE ≤ 5, we observed that the time-averaged LES solution obtained from utilizing the DNPS model is perform- ing remarkably successful in maintaining a low-level error over the multitude of scales for spatial shift. Therefore, we realized that unlike the PSM and DPSM model, the DNPS model does a great job in prediction of the high-order and two-point (nonlocal) statistics of ϕ̄ especially in over the inertial-convective subrange. In conclusion, we showed that employing larger filter sizes instead of conventional LES filter sizes in the proposed model is a promising direction. In this way, the high computational costs associated with fractional modeling can be compensated, and one can achieve stable, reasonably accurate, and fast simulations using the new closure model. • In chapter 5, an improved version of the dynamic fractional and tempered fractional SGS modeling [158, 166] is derived based on the following steps: (i) applying tempered Lévy stable distributions at the kinetic level (ii) obtaining the SGS forces based on the 101 tempered fractional Laplacian operator (iii) implementation of the new dynamic proce- dure for the computation of the model constant. Each scenario’s optimum parameters are determined when the predicted and ground truth SGS forces are correlated the most. The newly developed model is tested based on both a priori and a posteriori analysis for the LES and VLES cases. In the a priori tests, we analyzed the perfor- mance of the models in the prediction of the SGS forces and the SGS dissipation for the isotropic and anisotropic test cases. The non-isotropic test case is set up based on a multi-layer temporal shear layer turbulence case in which there are dominant anisotropic directions. The results of this section indicate that the new model does have some improvement in the prediction of the different force quantities for the HIT test case in comparison to the dynamic non-tempered and dynamic Smagorinsky mod- els for both the LES and VLES cases. The results of the multi-layer temporal shear layer test case indicate that the new model has almost the same performance as the non-tempered model using Lδ = 10 for the LES region. However, by enlarging the characteristic filter size to Lδ = 20, the new model performs better in the prediction of the statistics, especially the SGS forces. In the a posteriori analysis of the new model, we utilized a highly accurate spectral-based solver to test the numerical stability and performance of the model in preserving the correct statistics after long integration times. The results of this section indicate that the new dynamic tempered model has better performance in the prediction of the resolved turbulent kinetic energy for a de- caying HIT test case. Thanks to the tempered nature of the model, there is also good agreement for the prediction of the energy spectrum for the new model, especially in the dissipation range. The last test represents the capability of the models in pre- serving the high-order structure functions. Based on the results of this test, the new model shows significant potential, particularly when larger filters are used. It appears that the next generation of turbulence models might be based on combining AI-based models with nonlocal models that provide higher correlations. It is also suggested to 102 use the learning-based approaches for the determination of the fractional and temporal parameters in different modeling scenarios. 6.2 Future Works In the future, there will be some open topics to address in this field. The following are some of them: • Constructing a nonlocal turbulence model for temporal mixing layer and bound- ary layer problems: One can generalize and expand the current derivation process of ho- mogeneous isotropic flows to the turbulent flows near a solid walls. This section will allow to test the performance of the nonlocal models near the solid walls and cases with significant anisotropy. • Nonlocal Deep Neural Network (NDNN) surrogate modeling: We will train a deep neural network model based on the dense and convolution architectures to account for the nonlocal effects of the turbulent transport. Almost all the deep learning based turbulence models are being trained using the local approaches and they don’t address the true nonlocal physics especially in the bigger filter sizes. 103 BIBLIOGRAPHY [1] Michael Wilczek, Dimitar G Vlaykov, and Cristian C Lalescu. Emergence of non- gaussianity in turbulence. In Progress in Turbulence VII, pages 3–9. Springer, 2017. [2] 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. [3] P Bradshaw. Agardograph, no. 169. Nato Science and Technology Organisation, USA, 1973. [4] Anthony James Merrill Spencer and Ronald S Rivlin. The theory of matrix polyno- mials and its application to the mechanics of isotropic continua. Archive for rational mechanics and analysis, 2(1):309–336, 1958. [5] Anthony James Merrill Spencer and Ronald S Rivlin. Further results in the theory of matrix polynomials. Archive for rational mechanics and analysis, 4(1):214–230, 1959. [6] JL Lumley. Toward a turbulent constitutive relation. J. Fluid Mech, 41(2):413–434, 1970. [7] SBj Pope. A more general effective-viscosity hypothesis. Journal of Fluid Mechanics, 72(2):331–340, 1975. [8] 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. [9] 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. [10] Marta D’Elia, Qiang Du, Christian Glusa, Max Gunzburger, Xiaochuan Tian, and Zhi Zhou. Numerical methods for nonlocal and fractional models. arXiv preprint arXiv:2002.01401, 2020. [11] Ehsan Kharazmi and Mohsen Zayernouri. Fractional sensitivity equation method: Ap- plication to fractional model construction. Journal of Scientific Computing, 80(1):110– 140, 2019. [12] 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. [13] Yongtao Zhou, Jorge L Suzuki, Chengjian Zhang, and Mohsen Zayernouri. Implicit- explicit time integration of nonlinear fractional differential equations. Applied Numer- ical Mathematics, 156:555–583, 2020. 104 [14] 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. [15] Anna Lischke, Mohsen Zayernouri, and Zhongqiang Zhang. Spectral and spectral ele- ment methods for fractional advection–diffusion–reaction equations. Numerical Meth- ods, 157, 2019. [16] Seyedhadi Seyedi. Multiresolution solution of burgers equation with b-spline wavelet basis. arXiv preprint arXiv:1812.10117, 2018. [17] Mohsen Zayernouri and George Em Karniadakis. Fractional spectral collocation method. SIAM Journal on Scientific Computing, 36(1):A40–A62, 2014. [18] Ehsan Kharazmi, Mohsen Zayernouri, and George Em Karniadakis. Petrov–Galerkin and spectral collocation methods for distributed order differential equations. SIAM Journal on Scientific Computing, 39(3):A1003–A1037, 2017. [19] Mohsen Zayernouri and George Em Karniadakis. Fractional Sturm–Liouville eigen- problems: theory and numerical approximation. Journal of Computational Physics, 252:495–517, 2013. [20] Ehsan Kharazmi and Mohsen Zayernouri. Fractional pseudo-spectral methods for distributed-order fractional pdes. International Journal of Computer Mathematics, 95(6-7):1340–1361, 2018. [21] Julius Oscar Hinze, RE Sonnenberg, and Peter Johann Henry Builtjes. Memory effect in a turbulent boundary-layer flow due to a relatively strong axial variation of the mean-velocity gradient. Applied Scientific Research, 29(1):1–13, 1974. [22] Robert H Kraichnan. Direct-interaction approximation for shear and thermally driven turbulence. The Physics of Fluids, 7(7):1048–1062, 1964. [23] Wen Chen. A speculative study of 2/ 3-order fractional Laplacian modeling of tur- bulence: Some thoughts and conjectures. Chaos: An Interdisciplinary Journal of Nonlinear Science, 16(2):023126, 2006. [24] Brenden P Epps and Benoit Cushman-Roisin. Turbulence modeling via the fractional Laplacian. arXiv preprint arXiv:1803.05286, 2018. [25] Peter W Egolf and Kolumban Hutter. Fractional turbulence models. In Progress in Turbulence VII, pages 123–131. Springer, 2017. [26] Fujihiro Hamba. An analysis of nonlocal scalar transport in the convective boundary layer using the green’s function. Journal of Atmospheric Sciences, 52(8):1084–1095, 1995. [27] Fujihiro Hamba. Nonlocal analysis of the reynolds stress in turbulent shear flow. Physics of Fluids, 17(11):115102, 2005. 105 [28] Peter William Egolf and K Kutter. Nonlinear, nonlocal and fractional turbulence. Graduate Studies in Mathematics. Springer,, 2020. [29] 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. [30] Mehdi Samiee, Ali Akhavan-Safaei, and Mohsen Zayernouri. Tempered fractional LES modeling. Journal of Fluid Mechanics, in press, (arXiv preprint arXiv:2103.01481), 2021. [31] JP Laval, B Dubrulle, and S Nazarenko. Nonlocality and intermittency in three- dimensional turbulence. Physics of Fluids, 13(7):1995–2012, 2001. [32] 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. [33] Brendan Keith, Ustim Khristenko, and Barbara Wohlmuth. A fractional pde model for turbulent velocity fields near solid walls. Journal of Fluid Mechanics, 916, 2021. [34] 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. [35] Alexander V Milovanov and Jens Juul Rasmussen. A mixed soc-turbulence model for nonlocal transport and lévy-fractional fokker–planck equation. Physics Letters A, 378(21):1492–1500, 2014. [36] Hani Ali. Theory for the rotational deconvolution model of turbulence with fractional regularization. Applicable Analysis, 93(2):339–355, 2014. [37] Max Gunzburger, Nan Jiang, and Feifei Xu. Analysis and approximation of a fractional Laplacian-based closure model for turbulent flows and its connection to richardson pair dispersion. Computers & Mathematics with Applications, 75(6):1973–2001, 2018. [38] Massimo Germano, Ugo Piomelli, Parviz Moin, and William H Cabot. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A: Fluid Dynamics, 3(7):1760– 1765, 1991. [39] R. Klages, G. Radons, and I. M Sokolov. Anomalous transport: foundations and applications. John Wiley & Sons, 2008. [40] Peter Alan Davidson. Turbulence: an introduction for scientists and engineers. Oxford university press, 2015. [41] Peter William Egolf and Kolumban Hutter. Nonlinear, Nonlocal and Fractional Tur- bulence: Alternative Recipes for the Modeling of Turbulence. Springer Nature, 2019. 106 [42] 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. [43] JH Gerrard. The mechanics of the formation region of vortices behind bluff bodies. Journal of fluid mechanics, 25(2):401–413, 1966. [44] GH Koopmann. The vortex wakes of vibrating cylinders at low reynolds numbers. Journal of Fluid Mechanics, 28(3):501–512, 1967. [45] 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. [46] 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. [47] 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. [48] 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. [49] 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. [50] 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. [51] 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. [52] Marcello Lappa. Rotating thermal flows in natural and industrial processes. John Wiley & Sons, 2012. [53] 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. [54] 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. [55] 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. [56] J. P. Gollub and Harry L. Swinney. Onset of turbulence in a rotating fluid. Phys. Rev. Lett., 35:927–930, Oct 1975. [57] 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. 107 [58] Yufeng Lin, Jerome Noir, and Andrew Jackson. Experimental study of fluid flows in a precessing cylindrical annulus. Physics of Fluids, 26(4):046604, 2014. [59] 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. [60] 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. [61] 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. [62] 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. [63] 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. [64] J. M. Lopez and F. Marques. Instabilities and inertial waves generated in a librating cylinder. Journal of Fluid Mechanics, 687:171–193, 2011. [65] 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. [66] Juan M. Lopez and Francisco Marques. Rapidly rotating precessing cylinder flows: forced triadic resonances. Journal of Fluid Mechanics, 839:239–270, 2018. [67] Armen Der Kiureghian and Ove Ditlevsen. Aleatory or epistemic? does it matter? Structural Safety, 31(2):105 – 112, 2009. Risk Acceptance and Risk Communication. [68] 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. [69] 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. [70] 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. 108 [71] Chunsong Kwon and Daniel M. Tartakovsky. Modified immersed boundary method for flows over randomly rough surfaces. Journal of Computational Physics, 406:109195, 2020. [72] 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. [73] 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. [74] 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. [75] 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. [76] 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. [77] 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. [78] 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. [79] 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. [80] 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. [81] 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. [82] 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. 109 [83] Mirmehdi Seyyedi and Bita Ayati. Design and optimization of a sequential and hybrid advanced oxidation process system using response surface methodology. Journal of Applied Water Engineering and Research, pages 1–13, 2022. [84] 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. [85] 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. [86] 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. [87] 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. [88] Dongbin Xiu and George Em Karniadakis. Modeling uncertainty in flow simulations via generalized polynomial chaos. Journal of computational physics, 187(1):137–167, 2003. [89] Hessam Babaee, Sumanta Acharya, and Xiaoliang Wan. Optimization of Forcing Pa- rameters of Film Cooling Effectiveness. Journal of Turbomachinery, 136(6), 11 2013. [90] 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. [91] Hessam Babaee, Minseok Choi, Themistoklis P. Sapsis, and George Em Karniadakis. A robust bi-orthogonal/dynamically-orthogonal method using the covariance pseudo- inverse with application to stochastic flow problems. Journal of Computational Physics, 344:303 – 319, 2017. [92] 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. [93] Fatemeh Habibi, Mirmehdi Seyyedi, and Bita Ayati. Synthesis and application of reusable and magnetic rgo/fe3o4 nanocomposites in br46 removal from an aqueous solution; future prospects of an efficient adsorption platform. J Mater Environ Sci, 13(08):900–913, 2022. [94] Bernt Oksendal. Stochastic differential equations: an introduction with applications. Springer Science & Business Media, 2013. 110 [95] George Karniadakis and Spencer Sherwin. Spectral/hp element methods for computa- tional fluid dynamics. Oxford University Press, 2013. [96] 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. [97] 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, and Spencer J. Sherwin. Nektar++: Enhancing the capability and application of high- fidelity spectral/hp element methods. Computer Physics Communications, 249:107110, 2020. [98] 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. [99] 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. [100] 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. [101] 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. [102] 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. [103] 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. [104] 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. [105] 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. 111 [106] I.M Sobol’. On sensitivity estimation for nonlinear mathematical models. Matem. Mod., 2(1):112 – 118, 1990. [107] Olivier Le Maı̂tre and Omar M Knio. Spectral methods for uncertainty quantification: with applications to computational fluid dynamics, 2010. [108] 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 (submitted). [109] 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. [110] 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. [111] 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. [112] Siegfried Grossmann, Detlef Lohse, and Chao Sun. Velocity profiles in strongly turbu- lent Taylor-Couette flow. Physics of Fluids, 26(2):025114, 2014. [113] Dongbin Xiu and Daniel M Tartakovsky. Numerical methods for differential equations in random domains. SIAM Journal on Scientific Computing, 28(3):1167–1185, 2006. [114] George K Batchelor. Computation of the energy spectrum in homogeneous two- dimensional turbulence. The Physics of Fluids, 12(12):II–233, 1969. [115] 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. [116] M. Zayernouri and M. Metzger. Coherent features in the sensitivity field of a planar mixing layer. Physics of Fluids, 23(2):025105, 2011. [117] 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. [118] 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. 112 [119] 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. [120] P Clark Di Leoni, Tamer A Zaki, George Karniadakis, and Charles Meneveau. Two- point stress-strain rate correlation structure and non-local eddy viscosity in turbulent flows. arXiv preprint arXiv:2006.02280, 2020. [121] Xiang IA Yang, Sergio Pirozzoli, and Mahdi Abkar. Scaling of velocity fluctuations in statistically unstable boundary-layer flows. Journal of Fluid Mechanics, 886, 2020. [122] 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. [123] 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. [124] AM Akbarzadeh and I Borazjani. Controlling flow separation on a thick airfoil using backward traveling waves. AIAA Journal, pages 1–8, 2020. [125] Stephen B Pope. Turbulent flows, 2001. [126] Pierre Sagaut. Large eddy simulation for incompressible flows: an introduction. Springer Science & Business Media, 2006. [127] Charles Meneveau and Joseph Katz. Scale-invariance and turbulence models for large- eddy simulation. Annual Review of Fluid Mechanics, 32(1):1–32, 2000. [128] Parviz Moin, Kyle Squires, W Cabot, and Sangsan Lee. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A: Fluid Dynamics, 3(11):2746–2757, 1991. [129] Oleg V Vasilyev, Thomas S Lund, and Parviz Moin. A general class of commutative filters for les in complex geometries. Journal of computational physics, 146(1):82–104, 1998. [130] Stephen B Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New journal of Physics, 6(1):35, 2004. [131] 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. [132] 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. 113 [133] 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. [134] 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. [135] Mirmehdi Seyyedi and Bita Ayati. Treatment of petroleum wastewater using a se- quential hybrid system of electro-fenton and nzvi slurry reactors, future prospects for an emerging wastewater treatment technology. International Journal of Environment and Waste Management, 28(3):328–348, 2021. [136] Ali Akhavan-Safaei and Mohsen Zayernouri. A nonlocal spectral transfer model and new scaling law for scalar turbulence. arXiv preprint arXiv:2111.06540, 2021. [137] Jorge 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. arXiv preprint arXiv:2110.11531, 2021. [138] Joseph Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiment. Monthly weather review, 91(3):99–164, 1963. [139] Thomas S Lund and Hans-Jakob Kaltenbach. Experiments with explicit filtering for les using a finite-difference method. 1995. [140] Fady M Najjar and Danesh K Tafti. Study of discrete test filters and finite differ- ence approximations for the dynamic subgrid-scale stress model. Physics of Fluids, 8(4):1076–1088, 1996. [141] Ugo Piomelli and Junhui Liu. Large-eddy simulation of rotating channel flows using a localized dynamic model. Physics of fluids, 7(4):839–848, 1995. [142] Sandip Ghosal and Michael M Rogers. A numerical study of self-similarity in a turbu- lent plane wake using large-eddy simulation. Physics of Fluids, 9(6):1729–1739, 1997. [143] 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–309, 2010. [144] Sharath S Girimaji. Boltzmann kinetic equation for filtered fluid turbulence. Physical review letters, 99(3):034501, 2007. [145] 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. [146] Douglas K Lilly. A proposed modification of the germano subgrid-scale closure method. Physics of Fluids A: Fluid Dynamics, 4(3):633–635, 1992. 114 [147] Ali Akhavan-Safaei and Mohsen Zayernouri. A parallel integrated computational- statistical platform for turbulent transport phenomena. arXiv preprint arXiv:2012.04838, 2020. [148] Charles Meneveau. Statistics of turbulence subgrid-scale stresses: Necessary conditions and experimental tests. Physics of Fluids, 6(2):815–833, 1994. [149] Ugo Piomelli, William H Cabot, Parviz Moin, and Sangsan Lee. Subgrid-scale backscatter in turbulent and transitional flows. Physics of Fluids A: Fluid Dynam- ics, 3(7):1766–1771, 1991. [150] Boris I Shraiman and Eric D Siggia. Scalar turbulence. Nature, 405(6787):639–646, 2000. [151] Zellman Warhaft. Passive scalars in turbulent flows. Annual Review of Fluid Mechan- ics, 32(1):203–240, 2000. [152] Katepalli R Sreenivasan. Turbulent mixing: A perspective. Proceedings of the National Academy of Sciences, 116(37):18175–18183, 2019. [153] DA Donzis, KR Sreenivasan, and P Kc Yeung. Scalar dissipation rate and dissipative anomaly in isotropic turbulence. Journal of Fluid Mechanics, 532:199–216, 2005. [154] 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. [155] Gavin D Portwood, Balasubramanya T Nadiga, Juan A Saenz, and Daniel Livescu. Interpreting neural network models of residual scalar flux. Journal of Fluid Mechanics, 907, 2021. [156] Ali Akhavan-Safaei and Mohsen Zayernouri. A Nonlocal Spectral Transfer Model and New Scaling Law for Scalar Turbulence. arXiv preprint arXiv:2111.06540, 2021. [157] 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. [158] S Hadi Seyedi and Mohsen Zayernouri. A data-driven dynamic nonlocal subgrid-scale model for turbulent flows. Physics of Fluids, 34(3):035104, 2022. [159] Ugo Piomelli, Parviz Moin, and Joel H Ferziger. Model consistency in large eddy simulation of turbulent channel flows. The Physics of fluids, 31(7):1884–1891, 1988. [160] Antoine Vollant, Guillaume Balarac, and C Corre. A dynamic regularized gradient model of the subgrid-scale stress tensor for large-eddy simulation. Physics of Fluids, 28(2):025114, 2016. 115 [161] 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. [162] Ugo Piomelli. Large-eddy simulation: achievements and challenges. Progress in aerospace sciences, 35(4):335–362, 1999. [163] Bert Vreman, Bernard Geurts, and Hans Kuerten. Large-eddy simulation of the tempo- ral mixing layer using the clark model. Theoretical and Computational Fluid Dynamics, 8(4):309–324, 1996. [164] Franck Nicoud and Frédéric Ducros. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, turbulence and Combustion, 62(3):183–200, 1999. [165] S Hadi Seyedi, Ali Akhavan-Safaei, and Mohsen Zayernouri. Dynamic nonlocal passive scalar subgrid-scale turbulence modeling. arXiv e-prints, pages arXiv–2203, 2022. [166] Mehdi Samiee, Ali Akhavan-Safaei, and Mohsen Zayernouri. Tempered fractional les modeling. Journal of Fluid Mechanics, 932, 2022. [167] 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. [168] Guangping Li, Jalil Manafian, Mirmehdi Seyyedi, R Sivaraman, and Subhiya M Zey- nalli. Periodic, cross-kink, and interaction between stripe and periodic wave solutions for generalized hietarinta equation: Prospects for applications in environmental engi- neering. Advances in Mathematical Physics, 2022, 2022. [169] Rahul Agrawal, Michael P Whitmore, Kevin P Griffin, Sanjeeb T Bose, and Parviz Moin. Non-boussinesq subgrid-scale model with dynamic tensorial coefficients. arXiv preprint arXiv:2202.05502, 2022. [170] M Zayernouri and M Metzger. Coherent features in the sensitivity field of a planar mixing layer. Physics of Fluids, 23(2):025105, 2011. [171] Saad A Ragab and JL Wu. Linear instabilities in two-dimensional compressible mixing layers. Physics of Fluids A: Fluid Dynamics, 1(6):957–966, 1989. [172] Alfons Michalke. On the inviscid instability of the hyperbolictangent velocity profile. Journal of Fluid Mechanics, 19(4):543–556, 1964. [173] GP Neitzel and Stephen H Davis. Energy stability theory of decelerating swirl flows. The Physics of Fluids, 23(3):432–437, 1980. [174] 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. 116 [175] 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. [176] 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? arXiv preprint arXiv:1801.09767, 2018. [177] Elias M Stein. Singular integrals and differentiability properties of functions (PMS-30), volume 30. Princeton university press, 2016. 117 APPENDIX A 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 ρ − =− , (A.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 (A.1) can be solved through the Laplace transform on the variable t [173, 53]. 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), (A.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 [173] and [53] for derivations). Using equation (A.2) and implementing the same initial and boundary conditions in the numerical setup for a low Re number, a 118 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 A.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. comparison in different times was made (see Figure A.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. 119 APPENDIX B 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 [174, 175]. 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. 120 APPENDIX C FRACTIONAL-ORDER DIFFERENTIAL OPERATORS According to Lischke et al. [176], 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 (ξ) , (C.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 (ξ), (C.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, (C.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, (C.4) Rd |x − s| 22α Γ(α+d/2) where Cd,α = d/2 for α ∈ (0, 1] and Γ(·) represents Gamma function [176]. The π Γ(−α) α-Riesz potential is also formulated [177] as n  o Iα u(x) = (−∆)−α u(x) = F −1 |ξ| −2α F u (ξ) . (C.5) Considering (C.5), the Riesz transform is then given by n iξj  o Rj u(x) = ∇j I1 u(x) = F −1 − F u (ξ) , (C.6) |ξ| which is utilized in formulating the SGS scalar flux and D-FSGS models. 121