Discovering new scaling laws in turbulent boundary layers via multi-expression programming
Abstract
Flow turbulence modeling is an expensive computational operation that is often run with simulations using physical simplifications to reduce the cost. Large-eddy simulations (LES) of turbulent flow often make use of wall models in order to lower the computational cost of the simulation in the regions near solid walls, where typically the flow activity contains the smallest structures. Instead of resolving the boundary layer spatiotemporal scales, an algebraic expression for the fluid velocity field yields the wall shear stress, which is then used as a boundary condition to the outer flow, also known as wall modeling. Recently, novel wall model formulations are being developed using data-driven methods which exploit datasets generated with accurate wall-resolved large eddy simulation. However, it is still unclear what dimensional groups should be used to normalise such datasets when the wall shear stress is the target output, rather than a known quantity. The first goal of this dissertation is to explore such relevant groups across the dataset using machine learning applications as multi-expression genetic programming, a tool that allows to derive optimal expressions based on a population of initial candidates and a fitness function. The second goal is to be able to use it to find a correcting expression that yields the wall shear stress given other data and the defined groups.
Full text
BACHELOR FINAL THESIS Discovering new scaling laws in turbulent boundary layers via multi-expression programming Author: Irene Simó Director /Co-director: Arnau Miró Bernat Font Degree: Bachelor in Aerospace Engineering Examination session: Spring 2022 Document: Report
Page ii of 53 Abstract Flow turbulence modeling is an expensive computational operation that is often run with simulations using physical simplifications to reduce the cost. Large-eddy simulations (LES) of turbulent flow often make use of wall models in order to lower the computational cost of the simulation in the regions near solid walls, where typically the flow activity contains the smallest structures. Instead of resolving the boundary layer spatiotemporal scales, an algebraic expression for the fluid velocity field yields the wall shear stress, which is then used as a boundary condition to the outer flow, also known as wall modeling. Recently, novel wall model formulations are being developed using data-driven methods which exploit datasets generated with accurate wall-resolved large eddy simulation. However, it is still unclear what dimensional groups should be used to normalise such datasets when the wall shear stress is the target output, rather than a known quantity. The first goal of this dissertation is to explore such relevant groups across the dataset using machine learning applications as multi-expression genetic programming, a tool that allows to derive optimal expressions based on a population of initial candidates and a fitness function. The second goal is to be able to use it to find a correcting expression that yields the wall shear stress given other data and the defined groups. Keywords: wall model, LES, shear stress, dimensional group, machine learning
Page iii of 53 Acknowledgments My most sincere gratitude to my supervisors Dr. Arnau Mir´ o and Dr. Bernat Font, for all the guidance, support, criticism, and insight throughout these months. All the work has been possible thanks to them. A special acknowledgment to the CASE - Large-scale Computational Fluid Dynamics (Barcelona Supercomputing Center) research group for the occasional help that the rest of their researchers have given me. A hearty thank you to all the people I have met during this journey of obtaining the aerospace engineering diploma for being the core of this experience. To my family, for the support.
Page iv of 53 Declaration of honour This work is the final product of my undergraduate’s dissertation for the aerospace engineering program in the ESEIAAT Faculty, Universitat Polit` ecnica de Catalunya. The subject of study as well as the methodology were proposed to me by my tutors, who very gently have introduced me into and guided me through the complex physics of turbulence along the winter and spring months during which this thesis has been carried out. The dissertation is developed within the framework of the CASE - Large-scale Computational Fluid Dynamics (Barcelona Supercomputing Center) research group, focused in computational fluid mechanics and high performance computational mechanics. I declare that the work in this Bachelor Thesis is completely my own work, no part of this Bachelor Thesis is taken from other people’s work without giving them credit, all references have been clearly cited, I’m authorised to make use of the research group CASE - Large-scale Computational Fluid Dynamics (Barcelona Supercomputing Center) related information I’m providing in this document. I understand that an infringement of this declaration leaves me subject to the foreseen disciplinary actions by Universitat Polit` encica de Catalunya. Signature Irene Sim´ o Mu˜ noz [email protected] Barcelona, June 22nd, 2022
CONTENTS Page v of 53 Contents List of Tables vi List of Figures vii Nomenclature viii 1 Introduction and motivation 1 1.1 Problemdescription ........................................... 1 1.2 Goalsofthiswork ............................................ 2 1.3 Methodology ............................................... 2 1.4 Scopeandrequirements......................................... 3 2 Background and fundamentals 5 2.1 Fundamentals .............................................. 5 2.2 Literature ................................................. 6 2.3 Thelawofthewall ............................................ 8 2.4 Wall-boundedflowsinLES ....................................... 12 3 Methodology 15 3.1 Computationalapproach......................................... 15 3.2 Case study: channel flow at Reτ= 180 ................................. 20 3.3 Datasampling .............................................. 22 3.4 Profileanalysis.............................................. 24 3.5 Numericalimplementation........................................ 26 3.6 Datasetcreation ............................................. 28 3.7 MultiExpressionPorgramming ..................................... 30 4 Results and discussion 32 4.1 Profileanalysis.............................................. 32 4.2 Datasetanalysis ............................................. 36 4.3 Multi Expression Programming analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5 Conclusions and further work 42 5.1 Summary ................................................. 42 5.2 Conclusions................................................ 42 5.3 Furtherwork ............................................... 43 References 45
LIST OF TABLES Page vi of 53 List of Tables 1 Tablewithsublayers ........................................... 9 2 Shear velocities characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Simulationparameters.......................................... 21 4 Realtive errors between shear velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
LIST OF FIGURES Page vii of 53 List of Figures 1 Turbulent boundary layer with some eddies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Lawofthewall .............................................. 9 3 WallResolvedandWallModelmesh.................................. 12 4 Schematic of wall models functioning principles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 5 Tree diagram of the modeling techinques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6 Wall Model algebraic expression to solve for uτ............................ 16 7 Snapshot of different shear velocities for half a channel (δ)...................... 20 8 Simulationmeshshots.......................................... 21 9 Channelmesh .............................................. 21 10 Spatial Pearson’s coefficient plots for velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 11 Temporal Pearson’s coefficient plots for velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 12 Local Reynolds and shear velocity behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 14 Structureofthedataset ......................................... 28 13 FlowchartoftheUML .......................................... 29 15 Velocitygradientresults ......................................... 33 16 Error associated to the wall model shear velocity with wall distance . . . . . . . . . . . . . . . . . 34 17 Local Reynolds for different wall model shear velocities . . . . . . . . . . . . . . . . . . . . . . . 35 18 Data obtained for three wall distances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 19 Histogram of the MEP dataset variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 20 Scatter of the final dataset information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 21 Scatter of the correcting expression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
LIST OF FIGURES Page viii of 53 Nomenclature •wMagnitude at the wall δBoundary layer thickness κvon K´ arm´ an constant, see equation (2.10), page 12 uVelocity field BLog law constant, see equation (2.10), page 12 Re Reynolds number RelLocal Reynolds number ReτViscous Reynolds number ubBulk velocity ¯ •Filtered quantity h•i Averaged quantity µDynamic viscosity νKinematic viscosity ρDensity u+Non-dimensional streamwise velocity uτShear velocity y+Non-dimensional wall distance DNS Direct Numerical Simulation LES Large-Eddy Simulation RANS Reynolds-Averaged Navier Stokes equations WM Wall Modeled WR Wall Resolved
Introduction and motivation Page 1 of 53 1 Introduction and motivation 1.1 Problem description Turbulence is composed of structures of diverse scales that behave in a rather chaotic manner. The energy, velocity, and gradients in turbulent flows can take a wide range of values. Consequently, turbulent structures can differ a lot from each other. This makes them difficult to study because turbulent structures do not appear in any particular order, except that for wall-bounded flows, smaller structures tend to appear in the regions nearest to the solid wall. Another obstacle is that turbulent structures’ range of magnitudes is too wide to be covered entirely. Both large and small turbulent structures, as well as the entirety of the flow domain, can be studied with the Navier Stokes equations. As the flow fields contain an incredibly vast quantity of information, these equations are time-consuming to solve when used for turbulent flows, since they describe every detail of every turbulent structure; both large and small. This makes high-Reynolds-number flows (very turbulent flows) unbearable to model with the Navier Stokes equations. Advances on this topic have been limited by this. One of the lines of research uses a filtering operation to filter down to structures of a certain energetic scale only. This technique reduces the computational cost because only the larger scales are directly solved, whilst the smaller ones get filtered and are typically modeled. This approach is given the name of Large-Eddy Simulations or LES. Filtering by energetic scale has been a widespread technique to study the range of these structures. Sometimes, however, it is not an optimal approach to study the innermost region of the boundary layer, where, again, turbulent structures are smaller, since more detailed results require finer meshes. The finer the mesh, the better the turbulent information obtained. LES meshes that increase the density of nodes for the lower regions of the boundary layer - making them fines in this domain - are called wall resolved LES. These can provide exact solutions for the field in said region, but also require a higher computational cost than coarser meshes. Coarser meshes rely on the use of models to bear with the inner layer thus prioritizing a lower computational cost; the wall modeled LES. Wall models are mathematical and algorithmic entities capable of modeling the flow behavior in these innermost regions. These models can be used to complement coarser meshes that do not cover the boundary layer with such detail. With this technique, only the fluid structures down to the outer sublayers are directly resolved. The main line of research for the last decades has been to increase the accuracy of simulations at a lower cost; to improve the wall modeled simulations. In recent years many improvements have been made in the field of near-the-wall turbulence modeling. The various possible approaches differ notably as their formulation depends not only on the type of grid and domain used near the wall but also on the mathematical nature of the models, among other characteristics. In this sense, many articles have stridden away from classical approaches; see per instance wall modeling based on control theory (Nicoud et al., 2001, [9]) or on predictive neural networks (Yang et. al, 2019, [15]). All these wall models use the information gathered from the mesh to provide results for the inner layer region. However, because the meshes for these simulations are not as fine, the information they retrieve is not as accurate. Therefore, wall models still carry an error in their results. Wall resolved simulations, on the other hand, are much richer in local information in the inner sublayer. This makes them, as previously stated, more accurate than wall-modeled ones. Again, for wall-modeled simulations, the information provided by the model is less precise. Therefore, a need for a correction of the wall-modeled results arises. This correction would complement the function of the current wall models. With a correction like this, the wall model results would better match the wall resolved ones; an equally accurate simulation would be feasible at a lower computational cost.
2.3 The law of the wall Page 8 of 53 wall models, hence why LES requires less model tuning than other approaches such as RANS. These small structures are located in the dissipation range’s scales, energetically speaking. LES has both the reduced computational cost of not directly simulating every single structure while still explicitly representing the larger ones, which happen to be the least flow-dependent ones as they highly depend on the problem geometry. The filtering operation for LES begins with the separation of the velocity field into the filtered components ¯ u(x, t)(larger turbulent scales) and the residual - sub-grid - component u0(x, t)[11] u(x, t) = ¯ u(x, t) + u0(x, t) The governing equations are still the Navier Stokes, but spatial-filtered. Note that the filtering operation has some properties such as φψ 6=φ¯ ψ φ +ψ=¯ φ+¯ ψ φ06=0 Which have effect on the filtering of the Navier-Stokes equations (2.1). Once space-filtered, these can be rewriten [3]: ∇ · (¯ u) =0 ∂¯ u ∂t +∇ · (¯ u¯ u) = − ∇¯p+∇ · τ(¯ u) + ∇(¯ u¯ u−uu)(2.4) =− ∇¯p+∇ · τ(¯ u)− ∇τSGS The additional term ¯ u¯ u−uu is the stress tensor deriving from the residual or sub-grid velocity. This tensor is the entity to be modeled, as it is an unsolvable entity in the equations - note that the associated component u0(x, t)has been filtered out i.e. unresolved -. This tensor is τSGS =−(¯ u¯ u−uu)(2.5) =−νSGS (∇¯ u+∇¯ uT) + 1 3tr(τSGS)δ Where νsgs is the subgrid viscosity and δis the unitary tensor. Equation (2.4) can be rewriten: ∂v ∂t +∇ · (vv) = −∇˜p+∇ · (ν+νSGS )(∇v+∇vT)(2.6) Where vis an approximation of the filtered velocity ¯ uderived from (2.5). The pressure ˜pincludes the hydrostatic part of the Boussinesq’s relation. [1] Overall, equation (2.6) is still the original linear momentum equation from the Navier Stokes equations (2.2), but with a viscosity division. The solved viscosity ν, and the subgrid viscosity νSGS, that accounts for the smallest structures that remain unresolved. The equations are solved numerically for the filtered velocities v. For the unresolved entities, a LES closure model is used to evaluate the unresolved stress tensor, also called SGS tensor τSGS . For wall-bounded flows, modeling the subgrid stress tensor for the inner layer (region closest to the wall) might be insufficient for accurate result obtention, if the mesh is not fine enough. In these cases, wall models are used to model this region, on top of turbulence models. 2.3 The law of the wall The boundary layer is composed of different sublayers that present various behaviors. Within the thickness of the boundary layer δ, three regions can be identified; the inner layer (y/δ < 0.1), the buffer layer, and the defect layer or log-law region (y/δ > 0.1for both, the latter the outermost). For a laminar boundary layer, the viscous forces are predominant over the inertial ones. For a turbulent boundary layer, this only happens in the viscous sublayer; for the rest of the boudary layer (y/δ > 0.1) the inertial forces are more relevant thant the viscous ones.
2.3 The law of the wall Page 9 of 53 10−1100101102103104 0 10 20 30 y+ u+ Law of the wall u+=y+ u+=1 κln(y+) + B Figure 2: Law of the wall. Own elaboration with Spalding’s approximation [13] Table 1: Table with sublayers. Extracted from [11, Chap. 7] One of the quantities varying the most throughout the boundary layer in the wall-normal direction - and consequently in its sublayers, too - is the mean velocity profile in the streamwise direction. In kinematic and viscous equilibrium, this variation is represented throughout the widely accepted law of the wall, depicted in Figure 2. Equilibrium conditions (therefore, the law of the wall) apply for wall-parallel turbulent flows in the absence of adverse pressure gradients and geometrical discontinuities, such as flat-plate turbulent boundary layers, channels, and pipe flows. [1] The parameter uτor shear velocity, which requires the description of other values such as the bulk velocity or the shear stress at the wall to be defined conceptually. The law of the wall uses non-dimensional wall distance y+and streamwise velocity u+to describe this variation in order to provide a universal model. These groups are characterized by the shear velocity. The entire domain of a fluid can be considered to flow at an upstream velocity u∞. As there is one unique velocity capable of representing said flow, one can also define a Reynolds number for it. For a flow in a channel - this is, two walls parallel at h= 2δdistance -, some new descriptors can be defined. τwShear stress at the wall The shear stress is proportional to the velocity gradient at the wall. This proportionality is dependent on the fluid dynamic viscosity µ: τw=µ∂u ∂y y=0 (2.7) uτFriction or shear velocity From these wall shear stress τwand viscosity νone can define a viscous scale that is the appropriate
2.3 The law of the wall Page 10 of 53 velocity scales in the near-wall region. [11, Chap. 7] It can be writen as: uτ=rτw ρ(2.8) Sometimes with a change in notation, written as u∗ This magnitude is dependent on the gradient of velocity at the wall, which for turbulent flows can drastically change both instantly and spatially. ReτReynolds τ Based on the shear velocity, the Reynolds number can be writen: Reτ=uτδ ν(2.9) And Reichardt’s expression [12] states: Reτ= exp 1 0.88 ln(Re/0.09)(2.10) ubBulk velocity Once a body is met by the fluid, viscosity of said fluid originates the boundary layer. In particular, at the point of contact between the fluid and the solid body the fluid adheres completely to the solid. This is called the no-slip condition. Then the previously defined upstream velocity u∞and Reynolds number Re, which were global for the whole domain, must be redefined. Bulk velocity ubresults from the integration along the height of the channel 2δof all the streamwise velocities; ub=1 δZδ u(y)dy Re Reynolds number based on the bulk velocity The Reynolds number can be rewriten [11, Chp. 7]: Re =uL ν=ub2δ ν Both the bulk velocity and the Reynolds number based on the bulk velocity, give generic information about the flow, but fail to inform about the local behaviour. These groups are helpful to reduce the number of variables when describing the flow behavior inside the boundary layer, and also are useful to define non-dimensional groups. The great advantage of them is that the velocity profile scales with them and therefore they allow the obtaining of a universal law for the previously defined conditions. In particular, the law of the wall is described with non-dimensional wall distance and streamwise velocity; hu+iNon-dimensional average streamwise velocity The non-dimensional form of the streamwise velocity used in the law of the wall uses the averaged velocity and the shear velocity as the reference value, such that: hu+i=hui/uτ(2.11) y+Non-dimensional distance from the wall The non-dimensional form of distance from the wall used in the law of the wall uses the distance from the wall, the shear velocity and the kinematic viscotisity νas the reference values, such that: y+=y/δv=yuτ ν(2.12)
2.3 The law of the wall Page 11 of 53 Law of the wall The law of the wall is derived by assuming that the turbulence near a solid wall is a function only of the flow conditions pertaining at that wall and is independent of the flow conditions further away. Its derivation is simply a dimensional analysis that begins with Prandtl’s proposal: u uτ =fuτy ν In terms of averaged velocity, and rewriting the previous expression and introducing the viscous Reynolds eq. (2.9), one can adress the mean velocity changes as: hui=uτF0y δ, Reτ Where F0is a universal non-dimensional group to be determined. But it is the velocity gradient dhui/dy the most relevant quantity, dynamically speaking. Not all velocity variations are as relevant: ∂hui ∂y =∂hui ∂y ,∂hvi ∂y ,∂hwi ∂y ∂hui ∂y ∂hvi ∂y ,∂hwi ∂y And therefore it is better to write: ∂hui ∂y =uτ yΦy δv ,y δ(2.13) Where Φis a universal non-dimensional function. Note that the two parameters δvand δcontain the same information, but their scales are more appropriate for the viscous and outside the viscous sublayers respectively. One of the current challenges in fluid mechanics is to parametrize the law of the wall for the entire domain of the boundary layer. See in Figure that for both y+<5(linear region) and y+>30 (log-law region) there are accurate mathematical expressions representing the curves; u+=y+for the former and the log-law eq (2.14) for the latter; whilst for the buffer sublayer there is still no clear model. 2.3.1 Log-law In regards to the outer layer, where viscosity has the least impact, the log-law - see equation (2.14) - is an expression used to universally represent said profile. For a turbulent flow, it can be obtained from the law of the wall general expression eq. (2.13) with the assumptions: y δ1y+1 The previous conditions imply that the previously proposed non-dimensional function Φwill be constant, such that; Φy δv ,y δ→Φy+=1 κ And back to the original expression (2.13); dhui dy =uτ y 1 κ Where, substituting the expressions of y+(2.12) and u+(2.11); dhu+i dy+=1 y+κ Integration for y+results in the log-law: hu+i=1 κln(y+) + B(2.14) Where κ= 0.41 and B= 5.2are constants [11, Chp. 7] It is a time-averaged expression that provides an agreement between all near-the-wall flows starting from approximately y+>30 as seen in Figure 2, so that its mean velocities behavior can be compared in different scenarios.
2.4 Wall-bounded flows in LES Page 12 of 53 h Figure 3: Wall Resolved (left) and Wall Model (right) mesh. The nodes for a Wall Resolved simulation need to capture accurate information in the inner layer, so their node density increases in this region. Wall Modeled meshes rely on wall models and do only accurately capture information down until the exchange location h. The information below this distance contains errors. Own elaboration. 2.4 Wall-bounded flows in LES Because LES doesn’t explicitly simulate every structure, it has a lower computational cost. However, it still explicitly represents the larger structures, which are also the ones that depend most on the problem geometry and are hence the least flow-dependent. For wall-bounded flows, the challenge is to correctly caputre the velocity gradient in the viscous sublayer. In these scenarios, two main apporaches are used in LES; Wall Resolved and Wall Modeled flows. 2.4.1 Wall Resolved and Wall Modeled Wall Resolved simulations use a fine mesh for the entirety of the simulation domain. This mesh covers in more detail the wall distances down to the very wall. This approach successfully captures the velocity gradients at all wall distances and accurately models the smaller structures with the same precision as the far-field flow is resolved. This is because the inner layer is resolved explicitly, and for this, a key requirement is the first off-wall node to be at y+<1. A typical mesh for wall resolved simulations is depicted in Figure 3, left. Wall Modeled LES relies on a model for resolving the inner layer of fluid. The meshes used are not as fine as the wall resolved ones, and this causes the velocity gradient to be poorly retrieved in the viscous sublayer, as seen in Figure 3 (right). Although it also accurately resolves the larges structures and greater wall distances, it does not directly approach the closest to the wall region of fluid. Instead, the wall shear stress is modeled via a wall model (WMLES). However, the value of the wall shear stress is not correct because the wall modeled grid does not accurately capture the velocity gradient. Part of the present work requires the creation of a wall model. To check the accuracy of the results given by said wall model in the inner layer, data from a wall resolved simulation will be used. 2.4.2 Wall Modeling Wall Models in Large-Eddy Simulations (WMLES onwards), in their turn, can be of various approaches. A possible classification divides them into Hybrid LES/RANS and wall shear stress models, depending on the original simulation’s nature [6, pp. 5–8]. Hybrid LES/RANS combine RANS and LES domains to further reduce the computational cost, whilst wall shear stress models are based on the LES approach only. The latter is used in this work. For wall shear stress models, the grid used in the LES solver covers the full domain. For all simulations (free flows and wall-bounded) LES’s goal is to model the SGS stress tensor. The goal of the wall shear stress model, however, is to evaluate the wall shear stress only. In other words, the wall shear stress models’ goal is to model the velocity gradient, which is precisely one of the disadvantages of WMLES with respect to WRLES. This model takes the LES resolved velocity as the input and produces the wall shear stress as the output. This is then fed back to the LES solver as a boundary condition needed to numerically solve the filtered Navier Stokes equations, as depicted in Figure 4.
2.4 Wall-bounded flows in LES Page 13 of 53 yyyyyyy Wall Model h δu(y) ui, yi τw Figure 4: Schematic of wall models functioning principles. At the exhange location (h), the model retrieves data such as the streamwise velocity and wall distance and operates with it, returning the shear stress or an equivalent magnitude as the shear velocity uτ. Own elaboration adapted from [5,6] The main goal is to supply the accurate solving of the velocity gradient of WRLES with a wall model, which is less computationally demanding. These wall shear stress models are, in their majority, built on RANS equations. These are discretized separately from the main simulation, meaning the wall model is fairly independent and should be treated as modularly as possible. Depending on to what extent said models solve the RANS equations and if they do have wall-parallel grid connectivity, they can be divided into PDE, ODE, and algebraic [6, pp. 5–8]. Partial and Ordinal Differential Equations Partial Differential Equations solve the Reynolds Averaged Navier Stokes equations (or variants, such as the turbulent boundary layer equations (TBLE) or just the diffusive term) entirely and have grid connectivity in all coordinate directions. This is because they need a separate near-wall grid. Grid connectivity has similar requirements to the RANS grids’. Ordinal Differential Equations use some simplifications for RANS solving and do not have grid connectivity. They require numerical solving of the ODEs in the wall-normal direction. Data-driven Data-driven wall models have gained popularity in recent years. They typically use datasets from already known models and infer them to obtain new functional forms. These can then be injected into simulations to obtain more accurate predictions. [2] Algebraic They use analytical or known velocity profiles to obtain expressions that establish relationships between velocity and velocity gradient with shear stress, thus allowing direct solving. One of these is typically is the law of the wall. Exchange location hThe exchange location is the distance from the wall at which the LES and wall stress model interchange values. This distance, which contains the viscous and buffer sublayers, is modeled by the wall model. The value at which this happens is relevant because it must match the value at which the LES simulation begins to fail at providing accurate results. The wall modeled results at the exchange location are perfectly accurate but contain errors for smaller wall distances. This is because WMLES cannot capture the velocity gradient correctly, as previously stated. See Figure 3. If the exchange location is too close to the wall, the log-layer mismatch can appear, although it is still not very clear how this phenomenon originates [5]. In summary, the different approaches presented by the literature and the possibilities for LES are presented in Figure 5
2.4 Wall-bounded flows in LES Page 14 of 53 DNS LES Wall-bounded flows Wall Resolved (WRLES) Wall Models (WMLES) Wall shear stress PDE Data-driven Algebraic (Wall functions) Log-law Hybrid LES/RANS Free flows RANS Figure 5: Tree diagram of the modeling techinques. Own elaboration adapted from [1,6,8]
Methodology Page 15 of 53 3 Methodology This section covers the main workflow and work packages obtained in the development of this dissertation. Firstly, subsection 3.1 narrates the computational approach and expected challenges for this work. Secondly, the configuration, geometry and properties of the data used are presented in subsection 3.2. Then the statistical studies to ensure valid results are detailed in subsection 3.3, and next, the general overview of the study can be found in subsection 3.4. This latter subsection requires a great understanding of the different concepts previously covered in section 2. Later, in order to clarify the general computational implementation of the analysis performed as well to present the Python scripts developed to the reader, subsection 3.5 provides an overview of the whole script set. Finally, discussion on the proposed groups and dataset used for MEP is detailed in subsection 3.6. The main goal of this work is to create a dataset containing the optimal information so that the Multi Expression Programming algorithm can find a fitting expression for the shear velocity uτcorrection. 3.1 Computational approach The computational approach to log-law solving is filled with challenges and possibilities. This section covers the nature of the inputs given to the log-law, the solving mechanism, and some particularities to bear in mind - initial conditions, approximations, and such -, and sources of error. The issues covered are, when indicated, particular for the studied case; this is a channel and an algebraic wall model. To cover some of the aforementioned aspects of the computational implementation, the reader shall picture the given data as a wall resolved simulation’s results; velocity field (u), velocity gradient (∇u)and position of each of the mesh’ nodes (x, y, z)and fluid properties such as density and all viscosities (ρ, µ, ν). 3.1.1 Solving the log-law for the shear velocity Recalling the log-law expression eq. (2.14) and the non-dimensional groups u+(2.11) and y+(2.12), the log-law can be rewriten: hui uτ =1 κln yuτ ν+B(3.1) Where the solution for the shear velocity uτis not trivial as it cannot be yielded directly. This requires a solver, later discussed in this subsection. Note also that the given inputs, wall distance, and streamwise velocity, yand urespectively, can be diverse; particularly u, that can be instant, local, or averaged for both. Averaging possibilities are discussed in the following paragraphs of this section and more depth in section 3 as it is the main focus of this work. As stated, the developed Wall Model uses the forementioned log-law from equation (3.1) to obtain the shear velocity uτgiven the distance from the wall yand streamwise velocity uof a node. In order to do so, the equation is solved for uτusing the Newton-Raphson method via the scipy optimize newton [4] 2python function for computational cost optimization, tough preliminary work for this project included a Newton Raphson Python function as well. Some particularities to discuss when solving using Newton Raphson, and in general, any iterative solver, include convergence, the initial conditions, and the tolerance. Convergence and initial conditions are discussed in the following paragraph 3.1.2. For tolerance, a value that provides results in a reasonable time is set; for this project, 10e-8 was used. Tolerance is not one of the main sources of errors for log-law solving and is therefore not worth deeper analysis. In the case of the log-law equation (3.1), the natural logarithm adds one extra restriction for the wall, where y= 0. Because of the mathematical properties of the natural logarithm, an input of y= 0, therefore ln(0) is not possible. For the wall, however, the velocity is also zero u= 0 because of the no-slip condition, which makes uτstill unsolvable even if yielded from the linear region expression u+=y+. This restricts the applicability of both the log-law and the law of the wall to the first off-wall information, as the shear velocity for the wall will never be solvable. However, shear information at the wall is already given by the shear stress τwin eq. (2.7). 2See the documentation here
3.1 Computational approach Page 16 of 53 uτ f(x)Eq. (3.2) f0(x)Eq. (3.3) •Solution Figure 6: Wall Model algebraic expression to solve for uτ, eq. (3.2) and its derivative eq. (3.3) Of course, the linear expression in the viscous sublayer (i.e. where it holds), would be much more accurate. In the case of the log law, as seen in Figure 2, it only holds for y+>30, meaning that the inaccuracy for inputs closest to the wall will be high. The error varies along with the distance from the wall, being minimal in the exchange location and the log-law region (this is between y+≈30 up to y+≈500). In the exchange location, the error compared to the wall modeled results is minimal because the model provides the actual value used in that distance. 3.1.2 Initial guess for converging and correct iterative solving of the log-law One key aspect of iterative solving is wisely choosing a suitable initial guess; both for converging ensurement and for the desired solution-finding. As the Newton Raphson method is based on the derivative’s slope - whether this is negative or positive -, both convergence and correct solving are dependent on this parameter. Converging The initial guess’ position with respect to the solution is determined by the mathematical properties of the derivative in the sense that for a proposed solution x: •f0(x)<0or f0(x)>0Solution is at a different value. The value for which the derivative at that point is zero will be proposed as a solution for the following iteration. •f0(x) = 0Local minimum or maximum found; solution found. For a function with logarithmic properties as the one to be solved, a badly chosen initial guess can result in no convergence. To avoid this, one can take the original function eq. (3.1) and manipulate it to obtain the expression that the Newton Raphson solver will face. Substituting the adimensional terms and constants, we can write the function that the Wall Model solves for uτ: f(x) = 0 = 1 0.41 ln yuτ ν+ 5.2−hui uτ (3.2) Knowing the solver will also operate with the derivative; f0(x) = 0 = 1 0.41uτ +hui u2 τ (3.3) Next, the convergence analysis requires an evaluation of the tendencies these functions show. Figure 6 plots both the original and derivative functions for fixed values of uand y, irrelevant for the general analysis. Note that the derivative is always positive for the natural domain uτ>0, which is the only region where the log law (3.2) holds, as the logarithm only exists for said region. This means that for a given initial guess, the solver will always escalate to higher solutions. If the initial guess is greater than the desired solution, there will be no convergence.
3.1 Computational approach Page 17 of 53 Note too that the logarithm can never be negative. This means that an initial guess smaller than zero, or zero, is impossible to solve. For convergence, one can state: 0< x0≤xs Where x0is the initial guess and xs, the solution. Desired solution yielding The expected values of uτwill differ depending on the geometry of the problem as well as the fluidic conditions such as the bulk velocity. It is of great relevance to obtain a generic value of the shear stress uτbefore the Newton Raphson solving to be able to critically evaluate the obtained result. This can be done with the most general expression of uτusing expression (2.8) with the bulk velocity and Reynolds number. Note that this will yield the most general value of uτfor the studied scenario, therefore serving as a great value to use for contrast and testing. As seen in Figure 6, however, there is only one local minimum, therefore no other solutions will be found. The convergence analysis also highlights that for a greater than the desired solution initial guess, there will be no convergence. It is still a good praxis to compare the obtained uτvalue with the general one to ensure no errors are left undetected. As previously commented, given the mathematical properties of the function to solve and the way of working of the Newton Raphson method, the initial guess must always be smaller than the resulting root. Because this must be computed for various points in the profile at first, and the Newton result for all of them can vary notably (carrying different errors respectively, of course), the safest option is to take a fixed initial guess that can be used for all wall distances and velocities. Trying to find an initial guess smaller than the solution but close enough is inefficient and problematic for repetitive solving of the function. There is no limit when it comes to setting a smaller initial guess except that this must be positive, for two reasons: • From a computational and mathematical point of view, an initial guess x0<0would imply no convergence, as previously stated. Not only that but uτdefinition (2.8) states that the shear velocity is defined as the square root of the division of the shear stress and the density; a square root will never be negative. • From a physical point of view, again the definition of uτcancels any possibility of the value being negative. Even ignoring the mathematical impossibilities of uτbeing negative, an uτ<0would imply physical nonsenses as a negative distance to the wall, negative Reynolds number, and such. 3.1.3 Linearity and other approximations When using a finite element methodology for the simulation, some techniques used for computational cost reduction might have an impact on the results. In addition, the results given by the simulation are never continuous - the mesh can never be, by nature -, and while grid size is set to ensure the results’ quality, there still might be scenarios where approximations have to be made. Velocity gradient in the viscous sublayer One of the basic conditions of the gradient is its continuity. For the simpler resolution of the finite element method, a lumped mass matrix is sometimes used, as the system to solve is drastically simplified by using a diagonal matrix. This, however, is problematic for nodes with boundary conditions such as nodes at the wall. Consequently, in the case of the gradient, the values at the wall are not accurate. A simple solution for this is to use the first-order numerical approximation to compute the gradient for these nodes. Using the first-order finite-difference - as the fluid domain would not allow any other order finite difference, since this is being calculated at the very boundary -. f(x+ ∆x)−f(x) ∆x
3.4 Profile analysis Page 24 of 53 3.4 Profile analysis This subsection contains the discussion on what dimensional groups to select for the MEP algorithm and how to operate with the given data to obtain them. 3.4.1 Goal definition and formal description The basic step to properly approach the problem to solve is to establish a clear definition of the expected output; which is an expression capable of correcting the shear velocity obtained with wall-modeled data and instantaneous. For formalization and nomenclature simplification purposes, one can name the shear velocity to be corrected u∗ τ, and profile-local wall-resolved shear velocity: uτ. Then, the desired output is the function fsuch that: uτ=f(u∗ τ, ...)→uτ=f(u∗ τ,{γ}) The problem is not only to define such function but also to propose other variables relevant to such function, this is, propose dimensional a set of groups γto complement the u+ τor u∗ τ. To obtain a satisfactory function, both the mentioned function fand the proposed groups γmust meet a set of requirements: •Groups To make the solution universal, such groups must be non-dimensional. They must also be dependent on variables known before the solving of the wall model, this is, variables given by the LES solver at the exchange location or parameters of the simulation. This narrows down the space to only a set of variables that will serve as building blocks for the groups γ. γ=f(ρ, µ, ν, ub, Re, y, u)(3.6) Note that the possibilities do not only include non-dimensional combinations of these parameters but also averaging options must be taken into account. Setting limitations on this is relevant as using spatiotemporal, spatial, or temporal averaging would change the nature of the solution yielded. As the intention is to find an expression to locally correct the obtained friction velocity u+ τ, the streamwise velocity must be, at least, not spatially averaged. For the study, both instant uτand u∗ τare used too, but temporal averaging could be discussed in further work. Therefore, all values must be local and instant, if their nature allows it. It is relevant to note that there are some parameters known in the case study, a channel, such as the boundary layer thickness δ-2δ=hfor a channel -, that are not so easily calculated for other geometries. To retrieve a fairly universal result, such information must be considered unknown when constructing the function. •Function The function must be valid for, at least, the case study. Other than properties and requirements natural to mathematical functions and machine learning regressions, no specific conditions are set. Group proposal Focusing on extracting a simple expression as a basis, the first approach only includes only one group γ. The proposed group is the local Reynolds, defined as: Rel=yu ν=f(y, u)(3.7) It is only dependent on the wall distance and streamwise velocity as the study is performed with data retrieved from the same simulation; same fluid; therefore the kinematic viscosity νis constant. To study the impact of the local Reynolds only, the wall distance ymust be fixed and constant, so that the only variable is the velocity u. This means that the obtained expression will a priori be valid only for a certain wall distance. Of course, this is rather inconvenient as the wall distance is a highly particular parameter that can differ enormously for different geometries; therefore it cannot be selected directly. Rather, a non-dimensional wall distance y+should be selected as constant to ensure the general character of the yielded solution. For this, it is important to note that the selected non-dimensional wall distance, or reference non-dimensional wall distance y+ r, must belong in the outer layer, particularly where the log-law holds, this is: 30 < y+ r<30.000
3.4 Profile analysis Page 25 of 53 Furthermore, recalling the expression of the non-dimensional wall distance eq. (2.12), one can see a shear velocity is needed. Given that the wall resolved shear velocity is the magnitude to be known, and that the wall modeled one is the one to be corrected, the logical choice is to use the theoretical one. In addition, this is also the most generic one of them and the most easily calculated, given that it can be set even before running the simulation. In summary, to perform the first run firstly a reference wall distance y+ rmust be set. This will help yield, for each profile of study, a wall distance yusing the theoretical shear velocity uτ-, that will set both itself and its associated streamwise velocity local and instantaneous uas inputs for the wall model. The model, with said yand uwill retrieve a faulted u∗ τ, that will be corrected with an expression coming from a machine learning process to match the wall resolved uτ. This is depicted more schematically in the section 3.5, where a flow diagram of the whole dataset creation algorithm is provided. 3.4.2 Discussion on local Reynolds The only group proposed for a first approach is the local Reynolds, but this simple-looking number is worth a previous study as well. The local Reynolds is simply a non-dimensional group identical to the Reynolds number but resulting from local values. It can give information on the turbulent status of the flow locally, hence the reason why it is interesting to use as a characterizing parameter. For the studied region, where the log-law holds, one can write u uτ =1 κln yuτ ν+B With the definition of Rel=yuτ/ν Rel=yuτ ν1 κln yuτ ν+B Rel=y+1 κln y++B But for the sake of understanding how Relconditions the shear velocity - this is, how local values affect uτ-, note that the dependency between the local Reynolds Reland the shear velocity uτis heavily influenced by the coefficient of the wall distance and the kinematic viscosity, such that: Rel=fuτ,y ν This, given that is an argument to a natural logarithm, is of high relevance. Note that the proportionality between the local Reynolds and the shear velocity is conditioned by the value of this ratio, as for arguments smaller than 1 the logarithm will be negative, but for arguments greater than 1 it will be positive. This is, in addition, relevant, because a Reynolds number, per definition, cannot be negative as it lacks physical sense. This also indicates that for a single wall distance, the local Reynolds number is almost proportional to the shear velocity. For this, as the ratio of wall distance vs kinetic viscosity is now constant, one can name it C=y/ν. Note in Figure how the yielded equation form the log-law and fits perfectly the data, but also, how it is almost linear. This is due to the mathematical properties of the function: Rel=uτC1 κln(uτC) + B Rel=uτC κ[ln(uτ) + ln(C)] + uτBC Rel=uτln(uτ)A1+uτA2(3.8) Where: A1=y νκ A2=By ν+ ln y νy νκ This expression requires nodal values of the shear velocity, as the resulting Reynolds number is also locally nodal.
3.5 Numerical implementation Page 26 of 53 uτ Rel Figure 12: Local Reynolds and shear velocity behavior. Expression (3.8) for two values of the constants, showing how it changes with different wall distances. Figure 12 plots expression (3.8) for some values of the constants. As it shows, for a flow with a certain kinematic viscosity ν, different wall distances ydefine different curves. In summary, there is an existing and predictable relationship between the local Reynolds number and the wall modeled shear velocity for constant wall distance and kinematic viscosity, for the log-law region. This will likely have an impact on the obtained results and posterior correcting expression. 3.5 Numerical implementation For the given simulation data, work has been carried out before preparing the dataset to work with. For each snapshot, several profiles have been used for data gathering. To do so, the previous subsection described the problem to solve and the needed data to obtain to solve it. As discussed, the profile analysis will require the extraction of the values of the wall model shear velocity u∗ τ and local Reynolds number Relof each of the nodes studied. Given the information provided by the simulation results, several operations must be performed. Because the snapshots are separated in different files but all share the same mesh and properties, before allocating space for the dataset and filling it, some previous calculations and variable declarations shall be made. As each snapshot will be used to retrieve various profiles, it is useful to have a dictionary with arrays of node coordinates to easily select the same profiles for each snapshot. To do so, a snapshot is loaded and their node coordinates are copied into arrays. Filtering and sorting their components to ensure uniqueness and ascending order, these arrays are then stored. For the evenly distributed axis; xand z, which are the ones that determine the position of the profile selected, the spatial sted between nodes dx and dz respectively is also stored. This will serve as an auxiliary value to narrow down the selecting cube for each profile, where the cube containing two profiles is defined as: X=x−dx 2, x +dx 2Y= (0 −dy0, h +dy0)Z=x−dz 2, z +dz 2 Where X, Y, Z are the lengths of the cube in each direction, and dy0is any value that follows dy0>0. Any other calculations that do not depend on only local nor temporal values can be performed now too, such as the theoretical shear velocity uτt, only dependent on the channel Reynolds number Re, the kinematic viscosity νand the boundary layer thickness δ=h/2for a channel. Finally, the allocation for the dataset is also done before the iteration. As each node studied will have its data, the total nodes are calculated as n=nt st·nx sx·nz sz·ny Where each ncorresponds to the total available samples - number of snapshots t, number of nodes in the x direction, and number of nodes in the zdirection -, and each sis the temporal or spatial step. For the spatial
3.5 Numerical implementation Page 27 of 53 steps, as the number of nodes per direction is the same and the spatial correlations shared a very close-to-zero value for steps of 5 nodes, the criteria chosen is sx=sz. These fractions may carry computational difficulties if they are not exact, some adaptations - including extra space, for instance -, might be done, as all three quotients must be integrers. This is the condition that triggers the ceiling restriction on the calculation of the quotient. Note that nyis 2 since each column of nodes contains two profiles; one per each half of the channel. Once these are performed, temporal and spatial iterations can commence. For each instant, a new snapshot has to be loaded. And for each instant, various columns are studied; each named with the superscript iin the flowchart in Figure 13. For this, for each file various cubes are delimited by their coordinates and dimensions, as previously discussed, to select only one profile per iteration. Then, for each profile the following are calculated: τwIs necessary for the yielding of the wall resolved shear velocity, but is at the same dependent on the velocity gradient. y+Needed to compute for the same y+for each profile studied. As previously discussed in the Computational Approach subsection 3.1, from the array of non-dimensional wall distances computed for each profile, the nearest element to the fixated value selected is retrieved. The real wall distance associated with this element is used for the following computations. The theoretical shear velocity uτ, obtained in the previous calculations of the script, is used. ReτAlthough it is not needed for any posterior calculation, it provides a more intuitive understanding of the obtained shear stress. Note that for both the non-dimensional wall distance y+and the wall shear stress τw, one is still missing data. For the former, the calculations must be done for each half of the channel, and given that the coordinates in the wall distance direction yrange from 0 to h, it is impractical to obtain the wall distance for the upper half nodes. For the latter, the second-order finite difference approximation for the streamwise velocity - this is, the gradient - is missing; and the wall distances for the upper half present the same problem. To solve this, two functions are coded: •midflip Flips the first half of the input array and allocates it in the second half of the same array. When used with the wall distances - the coordinates vector one previously stored -, it returns an array with 0 to half channel wall distances, and then half channel to 0 wall distances. In essence, this array contains the distance to the nearest wall for all the nodes in the ydirection •gradv Computes the second-order finite difference approximation for the streamwise velocity and wall distance, thus returning the desired component of the velocity gradient. Both of them are run previously to obtain the τw,y+and Reτvalues. After this, the analysis must be done two times for each column, one for each half of the channel. As all the arrays at disposal at this point have the same length and are sorted so that each element corresponds to the same yposition, in ascending order from 0 to h, working with indexes is the simplest approach. Because of this, the first thing to do is find the two indexes associated with the fixated y+ r. For this, another function is coded •find nearest Which returns the index of the element in the given array which is closest to the indicated value. This function is run firstly for the lower half of the ycoordinates, then for the upper one, thus obtaining the indexes corresponding to the two nodes that will be studied. For each of the nodes, the following data is retrieved: yNeeded as an input for the wall model uNeeded as an input for the wall model uτW R Or the correct shear velocity. This is obtained from the previously calculated shear stress array, just by getting the values associated with either the lower or upper wall.
3.6 Dataset creation Page 28 of 53 The next step is to compute the shear stress calculated by the wall model uτ W M , of course, for which the Newton Raphson algorithm is used. The final step is to also compute the local Reynolds, easily calculated using the wall distance yand streamwise velocity uof the node, and the kinematic viscosity. The data is saved in the dataset, noted in Figure 13 with superscript d, and a new iteration kcan begin. A flow chart presenting the Unified Modeling Language (UML) diagram of the whole process is presented in Figure 13. 3.6 Dataset creation The main purpose of the dataset creation is to obtain an element that contains data structured in a manner such that it can be used for training and testing machine learning algorithms - a MEP algorithm, in particular -. However, the following methodology and description are focused solely on the allocation of space, saving of information, and structuration and order of the introduced data. This data is stored in a two-dimensional array that can well be used for applications other than machine learning, or that can be convenient for other purposes. The creation of a suitable dataset that not only is tractable by the MEP algorithm but also contains relevant groups for the expression desired is key to the success of all the work done until this point. The used MEP tool allows the input of information to be quite straightforward, therefore the requirement of the dataset to be tractable by the algorithm is broadened to be computationally speaking easy to operate with. For this reason, the selected data is stored in a two dimensional array and then exported to both text and .hdf5 files - in forms of a simple text structured in columns and a file with fields, respectively - when all columns of nodes and instants have been analysed. The retrieved data is stored in the allocated space such that: uτ0u∗ τ0|y+=100 Rel0|y+=100 uτh u∗ τh|y+=100 Relh|y+=100 uτ0u∗ τ0|y+=100 Rel0|y+=100 uτh u∗ τh|y+=100 Relh|y+=100 . . .. . .. . . uτ0u∗ τ0|y+=100 Rel0|y+=100 uτh u∗ τh|y+=100 Relh|y+=100 Test Train 3 nsamples Column (x, z)1 Column (x, z)2 Lower half profile (0< y < h/2) Upper half profile (h/2< y < h) Figure 14: Structure of the dataset The scheme in Figure 14 shows the order in which the samples are stored in the dataset, too. Note that the profile analysis operates by taking one snapshot, then analyzing all the selected profiles in it, then loading the following snapshot. This means that for every snapshot, it firstly selects one column, analyses both the lower and upper profiles, and moves to the next column. Because of this disposition, the lower and upper halves’ profiles get stored in an interleaved manner, each pair corresponding to one particular column. This allows the algorithm to minimize the loaded information at once, as it can replace the space of old profiles’ data as soon as it is not relevant and also establishes the pattern of assignment of odd and even indexes of the array to each half of the channel.
3.6 Dataset creation Page 29 of 53 uτt =Reν h/2 coord = x= [...] y= [...] z= [...] Import new snapshot yw i=(|yi−y0|yi≤h/2 |yi−yh|yi> h/2 Select new column of nodes xi;zi ui;yi uy|w=ui 1−ui w y1−yw uy=ui (y+h)−ui (y−h) 2h τw=µuy u(k) τ=µτw y+ i=yw i·u(k) τ/ν j∀yw(j) i≈y+ r y≤h/2y > h/2 yNR;uNR yNR =yw(j) uNR =ui(j) u∗(k) τ=u∗ τ0 f(u∗ τ) = 1/κ ln(yNRu∗ τ/ν) + B−uNR/u∗ τ ud(k) τ;u∗d(k) τ;Reld(k) ud(k) τ=u(k) τ u∗d(k) τ=u∗(k) τ Reld(k)=yi(k)·ui(k)/ν xi=xk+1 zi=zi xi=xi zi=zk+1 xi=xmax xi=x0zi=zmax midlfip gradv findnearest Yes Yes No No Newton Raphson (ui) (uy) idx jy≤h/2, jy>h/2 coord[y] (yw) Figure 13: Flowchart of the UML. The script ends when there are no more snapshots to import - not depicted in the flowchart for simplification reasons -, and the dataset is saved in both .txt and .h5 files. The shaded modules correspond to the previously described functions. Own elaboration.
3.7 Multi Expression Porgramming Page 30 of 53 For each loaded snapshot, each column is selected by delimiting a cube defined by their xand zcoordinates. Each column, which consists of a set of nodes from y= 0 to y=h, contains two boundary layers o profiles, one per half channel. The lower half is firstly studied; the shear stress at the lower wall uτ0, the wall model stress at y+= 100 is also retrieved (u∗ τ0|y+=100), and so is the local Reynolds Rel0|y+=100. The same analysis is repeated for the upper half, with the respective upper wall shear stress, wall model stress at y+= 100 and local Reynolds at such wall distance; uτh, u∗ τh|y+=100, Relh|y+=100. Once all the selected columns have been analyzed and all the associated results have been stored in the dataset, the next instant is imported and the spatial iteration begins again. The first column, the wall resolved shear stress, will be used as testing data as it is the magnitude to be predicted. The remaining columns, wall modeled shear stress and local Reynolds, will be the training data. 3.7 Multi Expression Porgramming Multi Expression Programming (MEP) is a type of genetic algorithm that yields a mathematical expression that describes a dataset. It uses an evolutionary approach such that multiple possible solutions are encoded in the same chromosome, a feature unique to MEP given that other genetic algorithms encode one solution per chromosome only. The MEP tool used has been created by the Barcelona Supercomputing Center Large-scale Computational Fluid Dynamics research group based on the MEPx library. Each chromosome is built with a set of functions - possible operations that the solution may contain - and a set of terminals - parameters or variables -. They are structured in genes such that a gene may be a terminal or a combination of other genes and a function, and the first gene must always contain a terminal. The classic example by Oltean (2006) [10] is: F={+,∗} T={a, b, c, d}C= 1. a 2. b 3.+,1,2 4. c 5. d 6.+,4,5 7.∗,3,5 8.+,2,6 Where Fand Tare the functions and terminals sets respectively, and Cis the chromosome with 8 genes. And the chromosomes are created by the MEP algorithm itself by initializing a population that meets parameters stated by the user. Therefore, the basic structures of a Multi Expression Programming are the chromosome, population, and functions or operations. For the used MEP tool, these can be defined, and have been set, as follows. •Chromosome As stated, it contains both the operations and variables and constants that form the search space. It also has storage space to allocate the values of the complexity of each expression as well as various types of fitness of the chromosome. Fitness are values that define the performance of the diverse expressions within the chromosome. Fitness Fitness can be measured as an overall value - best fitness given a certain complexity -, best fitness - best fitness regardless of the complexity - or best complexity - best fitness with minimum complexity -. Complexity The complexity of an expression is dependent on the number of constants, variables, and operations that it contains, as well as the type of operations - complexity associated with each operation and variable is set by the user -. This means the user can favor the appearance of certain variables or operations in the final expression. For the performed study, the following parameters are set: •Population The population contains all the chromosomes that will be taken into account when running MEP. As
3.7 Multi Expression Porgramming Page 31 of 53 many as possible are defined as testing or target chromosomes. The own algorithm can generate this population given come parameters defined by the user. •Functions The functions or operations are the possible mathematical operators that can appear in the resulting expression. For this study, these sets are the following: F={+,−,∗,ˆ, /, log,exp}T={u∗ τ, Rel} After the creation of a random population of individuals (chromosomes), the following loop of the MEP algorithm is repeated until reaching a maximum number of generations previously defined [10]. 1. Two parents are selected using a standard selection procedure, tournament selection in the used tool 2. Four offspring are obtained by recombining the parents based on the global fitness - two offsprings - and the best overall fitness of the chromosome - two remaining offsprings -. For the case of a single target, the shear velocity results in the same procedure. 3. The offspring are mutated. In the used MEP tool this is done by modifying the MEP chromosome based on probabilities set by the user. Uniform or one-point crossover is employed, the former is used in this work. 4. If the offsprings are better than the worst individuals in the population, they replace them for the next generation.
Results and discussion Page 32 of 53 4 Results and discussion This section covers the obtained results. Firstly, subsection 4.1 discusses the velocity gradients, shear velocity errors, and group selection while contrasting the previously defined methodology and concepts with the numerical data behavior. Next, subsection 4.2 narrates both the data obtained and how it changes with wall distance and the final inputs added to the dataset fed to MEP. Finally, subsection 4.3 describes the results provided by the genetic algorithm. 4.1 Profile analysis All the results presented for the profile analysis have been calculated for three wall distances, y+={1,20,100}, to contrast the results at different parts of the boundary layer. 4.1.1 Gradients The use of the finite difference approximation to calculate de gradient leads to the results discussed in this section. For this finite difference approximation, function gradv is called for each studied pair of profiles. As detailed in the previous Methodology section, the gradient is calculated using the second-order finite difference approximation. There are two nodes per column array in which this is not possible; the two wall nodes. The wall nodes are not operated with the second-order approximation, because they are the boundary nodes, this is, the last nodes of the column. Hence, a first-order approximation is used for them. The function computes the gradient and returns a column of gradient values. The plots in Figure 15 show the disparities between the values obtained with the finite difference approximation. The little dispersion for the node at the wall is fruit of the lack of noise resultant from turbulent structures; the node is at the very wall. The expected curve is precisely the mathematical tendency of a squared or a rooted value and shows the accuracy of the finite difference approximation. If we take into consideration the local Reynolds, it is easy to see why at the first off-wall node this curve should be that of a quadratic function: uτW M =µru1 y1 uτW M =µsRel1ν y2 1 On the other hand, this curve will gain dispersion as the node of study is at a greater distance from the wall. This is because performing finite differences at great wall distances fails to consider the activity of the flow between the wall and the selected wall distance. In between, turbulent structures interact with each other, thus adding noise or dispersion. This is again, clearly seen in Figure 15. Gradients for a non-dimensional wall distance of 20 and 100 present more dispersion. It is also worth highlighting in Figure 15 that the range of values for the gradient narrows - i.e. the difference between the maximum and minimum values the gradient takes - as the distance from the wall increases and that the maximum values also diminish as the wall distance is greater. This is expected behavior and is rooted in the principle behind the boundary layer. The typical streamwise velocity profile along with the wall distance already indicates that the ratio of change for said velocity is maximum at the near-the-wall region, and decreases as the fluid is more distant. For a flat plate, the upstream velocity - outside the boundary layer - could have a gradient equal to zero. 4.1.2 Shear velocities errors Theoretical, wall resolved and wall modeled Recall that there are three ways to calculate the shear velocity, and each of them is tied to each own limitations, which define their accuracy too. Notably, one of the most telling characteristics is how much data is needed to calculate each of them. For the theoretical one, only the setting parameters for the simulation are needed. The wall resolved one is solved by the LES solver and can be computed with averaged or local and instant inputs, and the wall model one requires inputs that again can
4.1 Profile analysis Page 33 of 53 0.04 0.06 0.08 0.1 0.12 10 20 30 uτW R ∂u/∂y y+=1 0.04 0.06 0.08 0.1 0.12 −5 0 5 10 uτW R y+=20 0.04 0.06 0.08 0.1 0.12 −1 0 1 2 uτW R y+=100 Figure 15: Velocity gradient results. From left to right, smaller to greater wall distance are shown, in concordance with the least to most dispersion for the gradient. Theoretical Wall resolved Wall model Value 0.0638 0.0676 0.0248 - 0.0685 Error 5.62% - 63.31% - 1.33% Table 4: Realtive errors between shear velocities. The results for each shear velocity are shown along the respective error. For the wall modeled, the most and least faulted results are displayed. Note that the first one, with an error of 63.31%, corresponds to the wall model result for the first off-wall node (y+= 1), whereas the second one, with an error of 1.33%, is for a node in the log-law region. This is because the log-law only holds in the logarithmic region of the law of the wall. be averaged or instant. Understanding that the wall resolved is the most accurate of the three and that the wall model might be accurate for the log-law region but not for the whole boundary layer domain is crucial. Table 4 shows this. Note also that the theoretical shear velocity uτis a reference to check whether the results are coherent or not, but it is important to understand that this is not only dependent on the performed analysis but also on other parameters as the mesh, the correct resolution of the viscous sublayer - node distance can be key - and more. This means that the error associated with the theoretical shear velocity is, in fact, no real error itself, but it can be used to measure with higher precision how accurate the wall model result is. Wall resolved and wall modeled This issue represents the basic idea behind this work. Although all shear velocities should ideally be in agreement, this is especially not true for local and instant values. This is precisely the problem to be solved with the correcting expression. However, even when averaged there is still a persistent error between the values. This is inherent in the different shear velocities’ nature. For instance, when studying the error between the wall resolved and wall model one, the error presented at each wall distance is intrinsically tied to the region of the boundary layer where each point rests. Figure 16 depicts a velocity profile calculated with the wall resolved shear velocity, and the associated error for its respective wall model shear velocity, thus illustrating this. This Figure 16 shows that the error is smaller in the greater wall distances, in the log law region. This is logical, as the wall model is, in essence, the log law. Recall that it only holds for non-dimensional wall distances such that 30 < y+<1000 approximately. Consequently, the expected error for distances below the exchange location is considerable. Furthermore, Figure 16 shows the error between these shear velocities for a spatio-temporal averaged pro-
4.3 Multi Expression Programming analysis Page 40 of 53 4.3 Multi Expression Programming analysis Some runs have been performed to obtain a correcting expression, all with similar results. Figure 21 depicts one of them, with the scattered data as well. No final expression has been found, as the obtained results are not entirely satisfactory and more work should be revised before running several more processes. As seen in Figures 18 for a non-dimensional wall distance of 100, there is notable dispersion in the used data. The used algorithm for Multi Expression Programming tends to smooth this by providing an expression that fits especially for the average values of the cluster, although it can also treat the dispersion as noise. In this case, the obtained expressions share some of the characteristics worth remarking: • Are mostly non-dependant on local Reynolds or the sehar velocity The expressions use only one variable, even when both the shear velocity and local Reynolds are proposed. This is due to the previously stated relationships, which are responsible for the data positioning in a plane-like surface, as seen in Figure 20. This also suggests that the multi-expression programming algorithm should be forced to use only one of them. As the wall resolved shear velocity is dependent on the local Reynolds, it seems suitable to use this. However, in relation to forcing the algorithm to use either one or the other, some runs have been performed with only the shear velocity or only the local Reynolds as a possible variable. The former have yielded results that differ very little in mathematical expression from the ones obtained when both variables were proposed as possible. The latter has yielded results that behave almost identically to Figure 21, but expressed only in terms of the local Reynolds, logically. Regardless of which of these two variables could be selected over the other, a possible conclusion is that more non-dimensional groups should be used, as it seems that none of these two quite fit the expression. A second discussion on the groups would be needed. • Is fitted in the average value of the target shear velocity The expression, when plotted, reveals itself to fit the middle of the cluster, almost like an average wall resolved shear velocity. This is because the MEP algorithm tends to smoothen the expressions it creates. This characteristic, along with the dispersion of the used data, quite limits the behavior of the proposed solutions. Overall, the obtained expression, as seen in the scatter plot, does a poor job at precisely predicting the corrected shear velocity. More work needs to be done to obtain a successful expression that can correct the local data. Although the multi-expression programming result presented is not entirely satisfactory, all the previously presented results are well defined and based on strong foundations, and compose a solid basis for obtaining a successful expression. Yielding a valid correcting expression from data is not trivial. Result analysis of the velocity profiles, gradient behavior, error, and relationship between the selected variables, all at different wall distances is key to improving the MEP results. Analysis of the distribution in the dataset and the plane-like dispersion is extremely relevant, too. All these results and the conclusions extracted from the obtained multi-expression programming expressions are useful for new iterations of the problem.
4.3 Multi Expression Programming analysis Page 41 of 53 0.04 0.06 0.08 0.10.12 1,500 2,000 0.04 0.06 0.08 0.1 0.12 WM shear velocity u∗ τLocal Reynolds Rel WR shear velocity uτ Data MEP expression Figure 21: Scatterplot of the correcting expression. The plot contains all the data used and the correcting expression evaluated for each sample, too.
Conclusions and further work Page 42 of 53 5 Conclusions and further work This final section serves as a closure, briefly reviewing the work carried out, then discussing the conclusions extracted from this project, and finally proposing future lines of work. 5.1 Summary This work has had as its main focus the creation of a dataset suitable for being exploited by a Multi Expression Programming algorithm and the obtaining of a correcting expression. To build such a dataset, previous treatment of the given data and calculations have had to be performed, as well as a discussion on what dimensional groups to include in the dataset. The previous calculations have started with a study of the correlation of the given data through Pearson’s coefficient values, firstly. From this, spatial and temporal steps have been selected to determine the sampling pattern used so that spatio-temporal independence between data is ensured. Then the selected data has been filtered and sorted following said pattern, and calculations regarding the fluid have been performed. These have included the theoretical shear velocity, velocity gradient using the second-order finite difference approximation, wall stress, wall resolved shear velocity, local Reynolds, and wall modeled nodal shear velocity, among other intermediate results. Except for the theoretical shear velocity, all have been calculated both locally and instantaneously. Once this has been done, a posterior study of how the extracted results are related and their dependence on the wall distance - this is, the sublayer to which they belong - has been carried out along a discussion on these parameters’ impact predictability, or lack thereof. Finally, a dataset has been built for the MEP algorithm to use. The algorithm’s parameters have been adjusted to limit the search space and performance factors, and an expression has been obtained and discussed. All these steps have lead to the main goals of this work. Using machine learning tools like multi-expression genetic programming, the first objective of this research has been to investigate such relevant groupings across the dataset. The second objective has been to be able to use this to find a corrective expression that, given other data and the designated groups, produces the wall shear stress. 5.2 Conclusions The study and analysis of the channel flow for a WRLES have yielded the creation of a dataset suitable for its exploitation and the obtention of a correcting expression capable of returning the value of the wall shear stress given the shear velocity and local Reynolds number at a certain wall distance. It is worth highlighting that this study has been carried out with local and instant values, which is in itself part of the mismatch between the shear stresses calculated and, therefore, the motivation behind the synthetization of a correcting expression. The inherent error of wall modeled results with respect to wall resolved ones is also to be corrected with the resulting expression. The dataset and corrective expression have only been calculated for a specific wall distance and Reynolds number. Describing the problem has yielded the conditions under which the solution must be searched for. In this sense, the understanding that not all the considered shear velocities (theoretical, wall resolved, and wall model) admit non-averaged values, and which of them are susceptible to more error depending on diverse factors, is key. This is, the wall resolved data is the most accurate, the wall modeled data is susceptible to correction only in the log-law region (where it holds, it would lack physical sense in other regions), and the theoretical data should only be used as reference. Along these, computational limitations have also been described. Newton Raphson convergence and initial guess requirements include that the initial guess is greater than zero but smaller than the expected solution. The computational implementation also lay the groundwork for posterior script development such as descriptions of the obstacles when working with a discrete set of data (disposal of nodal information) or the need to compute the velocity gradient via the finite difference approximation.
5.3 Further work Page 43 of 53 Once completed the analysis of the values and errors of each shear velocity, the next discussion had its focus on the selection of the groups to propose for the dataset used by the genetic algorithm as variables and its yielding. This arised the need of establishing a sampling pattern as well as describing the workflow of the algorithm to perform such analysis. The Pearson’s coefficients indicated the minimum steps to use for the pattern to ensure independent samples. The discussion on the group selection involved the proposal of two variables only - to obtain rather a simple solution -, and the setting of the wall modeled shear velocity (the entity to be corrected) and local Reynolds as selected groups. A posterior discussion on the local Reynolds relationship with other log-law variables gave a set of expected trends in the created dataset. The results, shown for three different wall distances, prompted the discussion on some shown behaviors, see Figure 18. Firstly, in regards to the wall model shear velocity values concerning the wall resolved ones, the most relevant issue was with the range of values difference. For greater wall distances the wall modeled velocities were higher, but their ranges were notably more narrow than that of the wall resolved shear velocities. The trends seen were different for each wall distance considered; the wall model shear velocity increasing with the wall resolved one for the first off-wall node and a wall distance of y+= 20, but remaining somewhat constant for a wall distance of y+= 100, in the log-law region. For the local Reynolds, the behavior is quite similar to the shear velocity (obtained with the wall model), but with different scaling. This is rooted in the fact that the local Reynolds number can be interpreted as a way of expressing the streamwise velocity, especially for constant wall distance and kinematic viscosity. Under these conditions, the log law can be written in the function of the local Reynolds instead of the streamwise velocity. This yields that for a constant wall distance the wall model shear velocity is a function of the local Reynolds. Overall, both of these magnitudes show the same behavior and yield a similar correcting expression. This is, they provide the same information. A new discussion on groups and variable selection should be done. Given these characteristics of the final dataset, the data exploited by the MEP algorithm can be analyzed accordingly. As seen both in Figure 18 and the histogram in Figure 19, the aforementioned difference in value range between the wall resolved shear velocity and the wall model one is notable in the created dataset and has had an impact on the regression. The final model obtained with the MEP algorithm returns an acceptable value of the wall resolved shear velocity that is coherent with the used data. However, because of the dispersion, the result is not very accurate, and the expression tends to return a value of the center of the cluster. It is not simple to generate a corrective expression from data. Improving the MEP results requires careful investigation of the velocity profiles, gradient behavior, error, and relationships between the chosen variables at various wall distances. This result analysis, which has been successfully done in this work, is a strong basis for future improvements. For this, the study of the dataset’s distribution and the plane-like dispersion is also crucial. All of these findings and the conclusions that have been drawn from the multi-expression programming expressions are helpful for future iterations of the issue. 5.3 Further work Along with both the description of the methodology and the results’ discussion, diverse questions have arisen. The following list presents a set of issues derived from this work that could be worth investigating further: •Expansion to other non-dimensional wall distances The decision of analyzing data from only one non-dimensional wall distance limits the generality of the obtained expression, as it is only applicable at a said distance. This process could be repeated for various other wall distances - it would be a rather automatic process as the scripts are already done -, or the variable of wall distance could be added to the terminal set in the genetic algorithm. The latter would require allocating more space when creating and saving the dataset. This would be an interesting addition to the model because it would widen its domain of use. This was tested in this work before setting a maximum of two variables for the correcting expression, and
5.3 Further work Page 44 of 53 the MEP algorithm was unable to fully perform any analysis. If implemented, this should be taken into consideration. •Verification of the model with other geometries The model has been created exploiting data from a channel, but should be verified firstly with other channels with the same settings as the used one, and also with other geometries. Of course, the number of applications would increase enormously if it held for other geometries too. If the model were only to be useful for channels, the repetition of this study for other geometries could yield a set of models. Both the use of a single model for more than one geometry or a set of models for a set of geometries would broaden this work’s result applicability and use. •Verification of the model with other Reynolds numbers For different Reynolds numbers, the flow would be more turbulent, was the Reynolds to be higher, or less turbulent, if the Reynolds were to be lower. This would directly affect the dispersion associated with behavior at greater wall distances, presumably making it more predictable for a lover Reynolds and less predictable for higher ones. The study and possible corrections of the model for other Reynolds numbers would broaden the applicability of the model, too. •Expansion to other non-dimensional groups Relative to both previous items, the proposal of groups could be redone to test how other variables characterize the shear velocity correction. In addition, for other geometries such as pipes or nozzles, other non-dimensional groups like the Nusselt number could be relevant.
REFERENCES Page 45 of 53 References [1] Joan Calafell Sandiumenge. Efficient wall modeling for large eddy simulations of general non-equilibrium wall-bounded flows. 2019. [2] Karthikeyan Duraisamy, Ze J Zhang, and Anand Pratap Singh. New approaches in turbulence and transition modeling using data-driven techniques. In 53rd AIAA Aerospace sciences meeting, page 1284, 2015. [3] Joel Guerrero. Les equations - filtered navier-stokes equations. [4] Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools for Python, 2001–. [5] Soshi Kawai and Johan Larsson. Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy. Physics of Fluids, 24(1):015105, 2012. [6] Johan Larsson, Soshi Kawai, Julien Bodart, and Ivan Bermejo-Moreno. Large eddy simulation with modeled wall-stress: recent progress and future directions. Mechanical Engineering Reviews, 3(1):15–00418, 2016. [7] Adri´ an Lozano-Dur´ an and Hyunji Jane Bae. Self-critical machine-learning wall-modeled les for external aerodynamics. arXiv preprint arXiv:2012.10005, 2020. [8] Timofey Mukha. Modelling techniques for large-eddy simulation of wall-bounded turbulent flows. PhD thesis, Acta Universitatis Upsaliensis, 2018. [9] Franck Nicoud, JS Baggett, Parviz Moin, and William Cabot. Large eddy simulation wall-modeling based on suboptimal control theory and linear stochastic estimation. Physics of fluids, 13(10):2968–2984, 2001. [10] Mihai Oltean. Multi expression programming–an in-depth description. arXiv preprint arXiv:2110.00367, 2021. [11] Stephen B Pope. Turbulent flows. Cambridge university press, 2000. [12] Hans Reichardt. Vollst¨ andige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ ur Angewandte Mathematik und Mechanik, 31(7):208–219, 1951. [13] DB Spalding. A single formula for the law of the wall. Journal of Applied Mechanics, 28(3):455–458, 1961. [14] Hendrik Tennekes, John Leask Lumley, Jonh L Lumley, et al. A first course in turbulence. MIT press, 1972. [15] XIA Yang, Suhaib Zafar, J-X Wang, and Heng Xiao. Predictive large-eddy-simulation wall modeling via physics-informed neural networks. Physical Review Fluids, 4(3):034602, 2019.