A Machine Learning Approach to Improve Turbulence Modelling from DNS Data Using Neural Networks^†

Round 1

Reviewer 1 Report

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
p.3: the value of the Reynolds number should be given for the plane channel flow. Also, the precise definition of the Reynolds numbers should be clearly stated, for both cases (which velocity?, channel half width as reference length?).

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.:
p.2: ANN acronym should be introduced.
p.2, eq.(1): this is an incompressible expression. The incompressible assumption is certainly valid in the present context, but it should be clearly stated.
p.3: the signification of p should be recalled.
p.3-4: when the simulations are introduced, it should be clearly specified that DNS is introduced first. Also, the span-wise dimension of the computational domains should be given.
p.4: "triangular grid", "quadrilateral grid". I think this refers to the 2D grid in the centre plane. This should be clearly stated.
p.4: are the k and w equations taken from Wilcox (1988)?
p.5: "...that can be set by the user and are ultimately are the result of a parameter search." Check phrasing.
p.8, Fig.5: the colour-scale should be explained. How many clusters are used?
p.11: "...in order to compensate for the lower values of k and still achieve..." An 's' is missing in 'achieves".
p.11: "For the serpentine the difference to the DNS in the length of the first separation bubble went from 2.98..." All the numerical results should be introduced before the conclusions. Are these values normalized? This should be clarified.

Author Response

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
    p.3: the value of the Reynolds number should be given for the plane channel flow. Also, the precise definition of the Reynolds numbers should be clearly stated, for both cases (which velocity?, channel half width as reference length?). -> This has been addressed specifying that “both cases are run at a Reynolds number of 5600, based on channel height and bulk mean velocity.”

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.:
    p.2: ANN acronym should be introduced. -> added
p.2, eq.(1): this is an incompressible expression. The incompressible assumption is certainly valid in the present context, but it should be clearly stated. -> rewritten in its general form as the formulation is not limited to incompressible
p.3: the signification of p should be recalled. -> added
p.3-4: when the simulations are introduced, it should be clearly specified that DNS is introduced first. Also, the span-wise dimension of the computational domains should be given. -> additional details are provided for the width of DNS and RANS domain
p.4: "triangular grid", "quadrilateral grid". I think this refers to the 2D grid in the centre plane. This should be clearly stated. -> clarified
p.4: are the k and w equations taken from Wilcox (1988)? -> correct, reference added
p.5: "...that can be set by the user and are ultimately are the result of a parameter search." Check phrasing. -> the reason for running a random search is explained more in detail later in the paper, so this sentence was removed, as it can be misleading
p.8, Fig.5: the colour-scale should be explained. How many clusters are used? -> added details
p.11: "...in order to compensate for the lower values of k and still achieve..." An 's' is missing in 'achieves". -> thank you!
p.11: "For the serpentine the difference to the DNS in the length of the first separation bubble went from 2.98..." All the numerical results should be introduced before the conclusions. Are these values normalized? This should be clarified. -> added a table to introduce the numbers before the conclusions and specified that all lengths are in non-dimensional units as per definition of the geometry

Reviewer 2 Report

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
- Page 2. Please introduce the full name for ANNs (artificial neural network). The reader not familiar with ANNs will have difficulty to understand this abbreviation.
- For the channel flow the turbulent shear stress obtained with c_1=-1 (fig. 4) falls in between the turbulent shear stress reproduced using DNS and the SST model. In contrast, the mean velocity profile reproduced with c_1 = -1 is completely wrong (different mass flow rate?). It is difficult to trust this result. What is a reason for this difference?
- Page 5 (top). It is not clear how the invariants and nondimensional functions provided at the end of Page 4 were applied to construct the ANN model. Please detail this. Motivate the selection of nondimensional functions.
- Page 6. It is said ” It is not necessary (and sometimes counter-productive) to include high strain regions and separation in the training set; it is important instead to capture the shear layers where most of the mixing happens” . Please explain better why the high strain regions have to be excluded from the training set. Is it due to technical difficulties in taking into account the high strain regions or different flow physics? I understand that the high strain regions in front part of the bubble do not necessarily generate high turbulent shear stress. Explain better the “training regions” in figure 5. Are they showed by dark grey or light grey zones ?
- Page 10. It is said that the ability to predict the shape of the separation bubble is influenced by model ability to predict the boundary layer separation point. I think it is also important to predict the reattachment point. In the rear part of the bubble the flow is characterized by a strong turbulence production (even higher than in the attached turbulent boundary layer). Without capturing the turbulent mixing in the rear part of the bubble the flow will not reattach. This is also visible in figure 9, where much stronger mixing is observed in DNS in the rear part of the bubble. This mixing is not fully captured by the “frozen field’. Provide better explanation why including the flow separation region is not beneficial for the model training. Is it beneficial to construct the ANN model separately for attached and separated flows?
- Use different line styles in figures 4 and 8. It is difficult to see difference between CH2 and DNS in figure 8.

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.:
- Page 10. It is said ”While for the channel this assumption is verified, for the serpentine there are some significant quantitative and qualitative differences (Figure 9): the levels of k are lower for the frozen field by a factor of 1.7”. I understand that in this case the ANN model was trained on the channel flow and it was applied for modeling of more challenging case (serpentine). In the “frozen field” model it is assumed that the mean velocity field, the mean velocity gradients and the turbulent shear stress are taken from DNS and all the rest should come from the k-transport equation. The Authors assumed the “closed” form of the dissipation and the turbulent diffusion terms in the k-equation. The dissipation is obtained by solution of the other transport equation (omega). The turbulent diffusion is modelled using the gradient diffusion hypothesis. But the DNS data are available. Why the “frozen field” result in figure 9 was constructed without taking into account the “exact” values of the dissipation and the turbulent diffusion terms? The agreement between DNS and “frozen field” might be better in figure 9 with other terms taken explicitly. In other words, the error which we see in figure 9 might not be only due to misrepresentation of the turbulent shear stress, but might also be due to an incorrect modelling of the dissipation and the turbulent diffusion terms.
- Page 11. It is said “In this paper we demonstrated how we can use high fidelity DNS data and a TensorFlow Machine Learning framework to generate Reynolds stress anisotropy formulations that replace the Boussinesq approximation.” I do not see the stress anisotropy in the present work. Only the single shear stress component is showed in fig. 4. Such a conclusion can be formulated if all Reynolds stress components are validated against DNS. Please reformulate this conclusion.

Author Response

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
    -    Page 2. Please introduce the full name for ANNs (artificial neural network).  The reader not familiar with ANNs will have difficulty to understand this abbreviation. -> addressed
-    For the channel flow the turbulent shear stress obtained with c_1=-1 (fig. 4) falls in between the turbulent shear stress reproduced using DNS and the SST model. In contrast, the mean velocity profile reproduced with c_1 = -1 is completely wrong (different mass flow rate?). It is difficult to trust this result. What is a reason for this difference? -> added a sentence to provide more details on why this is happening “It is worth noting that the way velocity profiles deviate from the linear behaviour near the wall depends on how close to the wall tau_12 starts to grow significantly.”. I hope this helps clarifying. Mass flow rate differences have been verified to be relevant re-running at the same pressure gradient and Reynolds number.
-    Page 5 (top). It is not clear how the invariants and nondimensional functions provided at the end of Page 4 were applied to construct the ANN model. Please detail this. Motivate the selection of nondimensional functions. -> added details on how the invariants and nondimensional functions constitute the input layer of the NN. As the choice of which features to use is better explained in a later part of the paper (they are one of the parameters of the random search) this has been rephrased in a more general sense.
-    Page 6. It is said ” It is not necessary (and sometimes counter-productive) to include high strain regions and separation in the training set; it is important instead to capture the shear layers where most of the mixing happens” . Please explain better why the high strain regions have to be excluded from the training set. Is it due to technical difficulties in taking into account the high strain regions or different flow physics? I understand that the high strain regions in front part of the bubble do not necessarily generate high turbulent shear stress. Explain better the “training regions” in figure 5. Are they showed by dark grey or light grey zones ? -> details on the training regions added; the reason to remove some of the clusters from the training is the result of the random search; this has been explained in more detail.
-    Page 10. It is said that the ability to predict the shape of the separation bubble is influenced by model ability to predict the boundary layer separation point. I think it is also important to predict the reattachment point. In the rear part of the bubble the flow is characterized by a strong turbulence production (even higher than in the attached turbulent boundary layer). Without capturing the turbulent mixing in the rear part of the bubble the flow will not reattach. This is also visible in figure 9, where much stronger mixing is observed in DNS in the rear part of the bubble. This mixing is not fully captured by the “frozen field’. Provide better explanation why including the flow separation region is not beneficial for the model training. Is it beneficial to construct the ANN model separately for attached and separated flows? -> rephrased as “the ability to predict the shape of the separation bubble mainly relies on having a good estimate of the turbulent mixing in the flow reaching the separation point and, more importantly, just outside the rear part of the bubble, as this strongly influences the reattachment point.”; having zonal models is an option but has not been investigated for now; this statement has been added.
-    Use different line styles in figures 4 and 8. It is difficult to see difference between CH2 and DNS in figure 8. -> agreed and addressed

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.:
    -    Page 10. It is said ”While for the channel this assumption is verified, for the serpentine there are some significant quantitative and qualitative differences (Figure 9): the levels of k are lower for the frozen field by a factor of 1.7”. I understand that in this case the ANN model was trained on the channel flow and it was applied for modeling of more challenging case (serpentine). In the “frozen field” model it is assumed that the mean velocity field, the mean velocity gradients and the turbulent shear stress are taken from DNS and all the rest should come from the k-transport equation.  The Authors assumed the “closed” form of the dissipation and the turbulent diffusion terms in the k-equation. The dissipation is obtained by solution of the other transport equation (omega). The turbulent diffusion is modelled using the gradient diffusion hypothesis. But the DNS data are available. Why the “frozen field” result in figure 9 was constructed without taking into account the “exact” values of the dissipation and the turbulent diffusion terms? The agreement between DNS and “frozen field” might be better in figure 9 with other terms taken explicitly.  In other words, the error which we see in figure 9 might not be only due to misrepresentation of the turbulent shear stress, but might also be due to an incorrect modelling of the dissipation and the turbulent diffusion terms. -> That is correct and the main purpose of the frozen field is indeed to obtain a k and omega field for the training based on the values that the RANS k-w equations would naturally converge to, when the flow field is the one obtained from DNS.  The reason for this is that basically we need the k-w  calculated from the standard k-omega SST (Wilcox 1988) equations: the solution to these equations are needed for the training because the standard k-w SST is still the back-bone of all revised anisotropy models in the present work. For this reason we cannot freeze other terms of the k-omega model to the DNS values. A different approach could be to look at the k and omega transport equation instead of the Reynolds stress – strain relationship and correct the way production and turbulent transport are modelled in the k-omega equations. This approach has been considered and requires some alterations to the current framework.
-    Page 11. It is said “In this paper we demonstrated how we can use high fidelity DNS data and a TensorFlow Machine Learning framework to generate Reynolds stress anisotropy formulations that replace the Boussinesq approximation.” I do not see the stress anisotropy in the present work. Only the single shear stress component is showed in fig. 4. Such a conclusion can be formulated if all Reynolds stress components  are validated against DNS. Please reformulate this conclusion. -> added in equation (1) the decomposition in isotropic and anisotropic part to later justify the use of the term anisotropy; for the turbulent channel in figure 4 tau_12 is the only term that plays a role in the dynamics of the velocity (as per eq.3); RANS calculations from which a posteriori results are taken have been run replacing the entire anisotropy formulation as per eq.2, hence the wording in the conclusion, where it is claimed that the anisotropy formulation is replaced, not validated.

Reviewer 3 Report

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
Review of Manuscript: ETC14-paper-704

Title: A Machine Learning Approach to Improve Turbulence Modelling from DNS Data Using Neural Networks

Authors: Y.F. Marioni, E. Alvarez de Toledo Ortiz, A. Cassinelli, F. Montomoli, P. Adami and R. Vazquez

The manuscript investigates by numerical means the capabilities of a more tolerant, viz. sophisticated, Reynolds stress – strain relation to predict two “simple” turbulent flows. The constitutive relation used relies on the effective viscosity hypothesis of Pope, proposed back in 1975, which allows to identify the Reynolds stresses exclusively through local quantities, in a more complex fashion than through a simple proportionality with the rate of strain. The authors are certainly aware of the fact that the complex effective viscosity hypothesis relies on severe restrictions on the flow features, that may vanish the remarkable additional effort (which is why nobody uses it …). I am here specifically referring to the necessity for the flow to be 1) high Reynolds 2) homogeneous (or nearly homogeneous) 3) respect certain proportionalities of all relevant scales 4) independent from the boundary conditions, except for the scaling parameters (your reference n. 7, Pope, JFM, V72, N2, pp 331-340, 1975). Therefore, the generality of the presented results dealing with two low Reynolds number cases is moderate, to say the least. Of course, the resolution requirements are what they are, and this is what it could be done with the present-day computational artillery. I appreciate that at least one of the test cases has only one homogeneous direction, yielding added value compared to the standard channel flow. Before entering a number of specific objections, let me say that I do not believe in this route, for the very simple reason that the possible correct prediction of the Reynolds stresses is not sufficient for a correct prediction of the deformation tensor and of the velocity field, given the non-linearity. Further, the a-priori training (I would call it doping) of the ANN with the DNS data originated from a simple, specific, low Reynolds, homogeneous flow can hardly serve more general cases. Finally, Reynolds stresses are not local quantities.
My proposal to this good team of researchers would be to follow the classical route, that is to calibrate the more advanced stress-strain relation, with your proposed system (or another one, it wouldn’t matter), on the basis of the standard building block flows, viz. channel, boundary layers, wakes and jets, for which a set of good quality DNS or q-DNS can be retrieved or built. This would drastically enhance the generality of your study.
Specific comments:
Page 4, item 2.
I understand the idea of solving the k and omega equations using the DNS frozen velocity field. More instructive would be to do the a-priori testing of all terms of the transport equations, as well as the k and omega variables themselves, to explore the magnitude of the approximations embodied in the closure. This has already been done for the channel, but certainly not for the return channel.
Also, the a-posteriori testing of the RANS models is scarcely presented, while many useful information could be readily withdrawn from the return channel case.
Page 8.
I am not sure I understand your channel calculations.
Usually, the equilibrium calculations are carried out with the help of a forcing term added to the streamwise momentum equation. Inspecting Figure 4, I can readily conclude that the bulk velocities of the calculations considerably differ, and so thus does the Reynolds number. Thus, you are comparing calculations at different Reynolds numbers. This is weird. I understand that one can either impose the mass flow and compute the friction, or vice versa, so that bulk and friction Reynolds numbers cannot simultaneously agree with the DNS data. But here the differences are huge, far more than allowed. This tells me that the way boundary conditions are imposed differ. The shear stress in figure 4, should be shown using inner scaling, and not outer scaling. Adding the total stress would help to quantify the deviation from the exact linear behavior, which can be inferred from the curly turbulent stress distribution. It also very weird that the k-omega prediction fails to match the DNS data by such a large amount. This result largely disagrees with Wilcox results (see figure 4.15 at page 134 of his book).
On the serpentine case it would have been instructive to compare the predicted Reynolds stresses together with DNS data.
Overall, this is a good paper which, with small additional efforts can become worth archiving.

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.:
/

Author Response

Mandatory Request Changes:Mandatory Changes: Requested changes which are essential for the understanding and completeness of the paper. Paper of author(s) who have not complied with these requests may be rejected.:
    Review of Manuscript: ETC14-paper-704

Title: A Machine Learning Approach to Improve Turbulence Modelling from DNS Data Using Neural Networks

Authors: Y.F. Marioni, E. Alvarez de Toledo Ortiz, A. Cassinelli, F. Montomoli, P. Adami and R. Vazquez

The manuscript investigates by numerical means the capabilities of a more tolerant, viz. sophisticated, Reynolds stress – strain relation to predict two “simple” turbulent flows. The constitutive relation used relies on the effective viscosity hypothesis of Pope, proposed back in 1975, which allows to identify the Reynolds stresses exclusively through local quantities, in a more complex fashion than through a simple proportionality with the rate of strain. The authors are certainly aware of the fact that the complex effective viscosity hypothesis relies on severe restrictions on the flow features, that may vanish the remarkable additional effort (which is why nobody uses it …). I am here specifically referring to the necessity for the flow to be 1) high Reynolds 2) homogeneous (or nearly homogeneous) 3) respect certain proportionalities of all relevant scales 4) independent from the boundary conditions, except for the scaling parameters (your reference n. 7, Pope, JFM, V72, N2, pp 331-340, 1975). Therefore, the generality of the presented results dealing with two low Reynolds number cases is moderate, to say the least. Of course, the resolution requirements are what they are, and this is what it could be done with the present-day computational artillery. I appreciate that at least one of the test cases has only one homogeneous direction, yielding added value compared to the standard channel flow. Before entering a number of specific objections, let me say that I do not believe in this route, for the very simple reason that the possible correct prediction of the Reynolds stresses is not sufficient for a correct prediction of the deformation tensor and of the velocity field, given the non-linearity. Further, the a-priori training (I would call it doping) of the ANN with the DNS data originated from a simple, specific, low Reynolds, homogeneous flow can hardly serve more general cases. Finally, Reynolds stresses are not local quantities.
My proposal to this good team of researchers would be to follow the classical route, that is to calibrate the more advanced stress-strain relation, with your proposed system (or another one, it wouldn’t matter), on the basis of the standard building block flows, viz. channel, boundary layers, wakes and jets, for which a set of good quality DNS or q-DNS can be retrieved or built. This would drastically enhance the generality of your study. -> This feedback is really appreciated, there are a number of restrictions in place and the authors are aware of them. However, the first step in this journey has been to put in place a framework that can gradually be expanded and has its goal in improving the design process in the context of the turbomachinery industry. It is agreed that based on what is shown in this paper the generality of the models can be questioned; nonetheless, working with these two test cases served the purpose of a proof of concept that has already revealed many technical features and procedural aspect of the method. This will be surely extended, first to additional canonical or simpler test cases (2D mixing, square duct, axi jet) and then move to more industry-representative cases (turbine 2D profile, wing tip vortex, 3D cascade, chevron nozzle). Ultimately the goal is not to identify a closure that can work for all fluid dynamics cases, but rather have a process in place to adapt the closure to a family of problems (e.g. turbine blades). For the near future we believe on the necessity to improve models based on the eddy viscosity hypothesis: this is still behind most of the CFD run in the industry and more complex formulations of the anisotropy such as the one in “An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows” from Wallin et al. are part of the backbone of the CFD capability in aerospace companies.

Specific comments:
Page 4, item 2.
I understand the idea of solving the k and omega equations using the DNS frozen velocity field. More instructive would be to do the a-priori testing of all terms of the transport equations, as well as the k and omega variables themselves, to explore the magnitude of the approximations embodied in the closure. This has already been done for the channel, but certainly not for the return channel. -> Understanding the role of the uncertainty of each of the terms in the transport equation is without doubt very instructive. It requires however an even more flexible framework than the one that has currently been put in place, where each term can be either frozen, or modelled based on standard RANS, or modelled with an additional correction term. So far this has been done only for one of the elements, viz the constitutive law and to some extent, although not presented in the paper, for the production term in the k equation. This would effectively be a separate branch of this research, which for now is mainly attempting to target industry-relevant problems with a step-by-step methodology.
Also, the a-posteriori testing of the RANS models is scarcely presented, while many useful information could be readily withdrawn from the return channel case. -> I totally agree with this, unfortunately the format of ETC conference papers only allows 12 pages (with very generous margins and blank spaces) and, not being able to reference previous work as this is the first paper we write on the subject, I felt like priority has been given to the presentation of the framework. More details about a posteriori results will follow in the future with more comparisons included.
Page 8.
I am not sure I understand your channel calculations.
Usually, the equilibrium calculations are carried out with the help of a forcing term added to the streamwise momentum equation. Inspecting Figure 4, I can readily conclude that the bulk velocities of the calculations considerably differ, and so thus does the Reynolds number. Thus, you are comparing calculations at different Reynolds numbers. This is weird. I understand that one can either impose the mass flow and compute the friction, or vice versa, so that bulk and friction Reynolds numbers cannot simultaneously agree with the DNS data. But here the differences are huge, far more than allowed. This tells me that the way boundary conditions are imposed differ. The shear stress in figure 4, should be shown using inner scaling, and not outer scaling. Adding the total stress would help to quantify the deviation from the exact linear behavior, which can be inferred from the curly turbulent stress distribution. It also very weird that the k-omega prediction fails to match the DNS data by such a large amount. This result largely disagrees with Wilcox results (see figure 4.15 at page 134 of his book). -> Channel calculations are run at a prescribed pressure gradient and depending on the Reynolds stress -strain constitutive law variations in velocity profiles can be very significant. Following the comment, calculations have been rerun now adjusting the kinematic viscosity to account for the difference in bulk velocity in order to achieve the same Reynolds number (5600). This has been added to the paper as “The RANS calculation was run imposing a pressure gradient along the duct, adjusting the molecular viscosity for c1 = -1 and k-omega SST in order to keep the same Reynolds number, as the bulk mean velocity differed more than 2% from the DNS.” Plots in figure 4 have been attempted in inner scaling, but they appear much harder to interpret without adding any significant piece of information. Probably due to having a different version of Wilcox’s book I am not able to check against figure 4.15 at page 134; however, most discrepancies are in or start from high curvature regions, for which correction terms have been investigated as a means to account for the lack in TKE production. We do not see the difference in prediction to DNS as unexpected and this class of flows is widely recognised as difficult for standard closures (e.g. U-bend and similar internal flows with separations), with RANS exhibiting a noticeable deficit in TKE.

On the serpentine case it would have been instructive to compare the predicted Reynolds stresses together with DNS data. -> As previously stated
Overall, this is a good paper which, with small additional efforts can become worth archiving.

Recommended Requested Changes:Recommended Changes: Changes will improve the quality of the paper. Authors are strongly encouraged to comply with these requests.: