{'fraction_non_alphanumeric': 0.024844894865388396, 'fraction_numerical': 0.007607871761039836, 'mean_word_length': 4.789371038517796, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A PDE-free, neural network-based eddy viscosity model coupled with RANS equations

Most turbulence models used in Reynolds-averaged Navier-Stokes (RANS) simulations are partial differential equations (PDE) that describe the transport of turbulent quantities. Such quantities include turbulent kinetic energy for eddy viscosity models and the Reynolds stress tensor (or its anisotropy) in differential stress models. However, such models all have limitations in their robustness and accuracy. Inspired by the successes of machine learning in other scientific fields, researchers have developed data-driven turbulence models. Recently, a nonlocal vector-cloud neural network with embedded invariance was proposed, with its capability demonstrated in emulating passive tracer transport in laminar flows. Building upon this success, we use nonlocal neural network mapping to model the transport physics in the k-epsilon model and couple it to RANS solvers, leading to a PDE-free eddy-viscosity model. We demonstrate the robustness and stability of the RANS solver with a neural network-based turbulence model on flows over periodic hills of parameterized geometries. Our work serves as a proof of concept for using a vector-cloud neural network as an alternative to traditional turbulence models in coupled RANS simulations. The success of the coupling paves the way for neural network-based emulation of Reynolds stress transport models.


Introduction
Turbulence is a physical phenomenon that is ubiquitous in natural and industrial flows. It is prevalent no matter in internal flows such as fluid inside pipes and channels or exterior flows like the air surrounding airplanes and vehicles. Turbulent flows are characterized by their chaotic nature and the wide range of length and time scales [1], which poses a great challenge in understanding and predicting them. Modelling turbulent flows has thus been an important research topic over the past half-century.
The exact description of fluid flows is given by the Navier-Stokes equations when the continuum assumption applies. However, the computational cost of directly solving the Navier-Stokes equations (DNS) is prohibitively high as it grows cubically (Re 3 ) with regards to Reynolds number [2]. Simulating turbulent flows at tractable costs requires additional modelling efforts. The most commonly used turbulent flow simulation techniques are Reynolds-averaged Navier-Stokes (RANS) models. The idea is to decompose the turbulent flow into the mean flow field described by the RANS equations (primary equations) and the fluctuating velocity field whose influence on the mean flow field is modelled by turbulence closures. The construction of closure models usually involves a combination of physical understanding and a calibration procedure accommodating simple canonical flows. This procedure, however, is limited by the capability of incorporating various flow configurations [3]. Recent development in machine learning techniques has opened up new avenues for solving the long-standing closure problems [4]. Thanks to the high expressive power of machine learning models such as neural networks, it is possible to make use of much larger datasets and obtain more generalizable and/or accurate turbulence models [5][6][7][8].

Invariances in machine-learned turbulence models
It has been observed that the prediction performance of machine-learned turbulence models can be significantly improved by embedding the known physics [9][10][11]. Most works along this line aim at embedding the symmetry of the turbulence flows, which refers to the invariance properties of the flow under transformations of the coordinate system. Specifically, the invariance properties include invariance under translation, rotation, and uniform motion (Galilean invariance) of the reference frame. Another property of turbulence flows that is seldom incorporated in data-driven turbulence models is its nonlocality. In most data-driven models, the inputs are only local flow variables (e.g., velocity gradient and its linear combinations) at the point of interest [10,12,13]. This simplification is only reasonable under the assumption that there is a local balance between the production, redistribution and the dissipation of the Reynolds stress [1]. However, as shown in transport equations of Reynolds stress and other turbulence quantities, turbulence quantities at one point are not only determined by the local flow field; instead, they are greatly influenced by the upstream flow structure and boundary conditions [14].
A vector-cloud neural network (VCNN) has been recently proposed [15] to incorporate the nonlocality.
The idea of the VCNN framework is to use the information from a set of points surrounding the point of interest as the input. One challenge for using a group instead of a single point is to guarantee permutation invariance [16]. This invariance comes from the fact that the elements in the input set have no intrinsic ordering, and thus the outputs should depend on the set as a whole but not on the specific ordering of the elements. The permutation invariance brings another challenge in addition to the invariance properties required by the flow physics.
In this work, we drew inspirations from the VCNN and developed the neural network shown in Figure 1 to emulate the k-ε turbulence model, ensuring all the invariance properties are embedded. The network aims to establish a mapping from mean flow properties in the neighbourhood to turbulence quantities at the point of interest. The network consists of an embedding network and a fitting network, and they are connected by a linear transformation. In the embedding network, higher-level representations of the input feature set are extracted by embedding functions. The linear transformation between two networks guarantees the invariance properties. Finally, the fitting network takes in the invariant feature matrix and provides the final predictions of the targeted variables. Details of the methodology and proof of invariance properties are provided in Section 2. The advantages of the proposed framework are its nonlocality, strict invariance-preserving properties and good scalability [17]. Most of the current machine-learning-based turbulence models are local as mentioned earlier. In comparison, our model takes in information from the neighbouring region and is able to capture more physics in the real flows. Compared to the nonlocal models obtained with other approaches, such as Convolutional Neural Network (CNN) [18][19][20][21], the proposed model strictly ensures invariances. One of the major drawbacks of CNN is a lack of geometric invariance [22]. General practice for reducing CNN's deviation of invariance is augmenting the dataset [23], in which the input data is transformed by translation and rotation. The enlarged dataset is then fed into the network such that the network acquires certain robustness towards variations of the input data under different coordinate systems. However, as implied by the procedure, there is no guarantee of the exact invariance properties of the model. Minor violation of the symmetries may not influence the classification problem of image recognition, for which CNN was initially developed [24]. Nevertheless, in a physics modelling scenario, it may lead to intrinsically wrong results or even breakdowns of the simulation [17,25]. Considering the stringent requirement of physics modelling, having solid invariance properties is a major advantage of our framework. Besides, CNN is rarely employed to deal with data on nonuniform or non-Cartesian mesh, which are usually used in fluid simulations for local refinement or adapting complex geometries.
Whilst the previous work on VCNN successfully embedded the nonlocality and frame-independency into the network, there are limitations on the architecture and the application scenario. One line of development is to modify the architecture such that it can be compatible with tensor-based outputs. This is crucial in eventually using VCNN to emulate Reynolds stress transport equations. Such adaptation for equivariant tensor outputs of the network has already been addressed in a follow-up work [26]. Another line of devel-opment stems from the limitation that the performance of the network was only evaluated in predicting the transport of passive scalar in a laminar flow. Its performance in predicting turbulent flow fields is unknown considering the complexity of turbulent flows compared to laminar flows. Furthermore, the final goal of a turbulence model is to serve as a part of the fluid solver and provide stable predictions to the RANS equation solver. Therefore, the neural network model needs to be evaluated in a coupled setup with the RANS solver.
This paper aims to address these critical issues.

Major challenges in coupling neural network model to RANS Solvers
Compared to the previous work, the current work of predicting turbulence quantities in a coupled setup is more challenging, due to the interaction between the neural network predictions and the RANS equations.
The lack of stability is a common potential problem for neural network models. When coupled with the RANS equations and evaluated in the iteration process of the solver, the neural network model must deal with non-converged flow fields most of the time. Such inputs of mean fields typically are not present in the training data. Consequently, the neural network is required to perform extrapolation tasks during the iterations. The behaviour of the networks in such circumstances can be unpredictable. When the predictions provided by the neural network are fed back to the RANS equations, the instability is transported to RANS equations.
Studies have shown that the RANS equation can be very sensitive to the Reynolds stresses [27], which poses an additional challenge to the closure modelling. An alternative strategy for improving the robustness is to incorporate the RANS solver in the training process, but that would introduce major challenges in the training process by requiring adjoint solver [28] or ensemble simulations [29].
Apart from the difficulty of extrapolation, the turbulence quantities are more complicated than the concentration field in terms of flow physics. The turbulence quantities are described by two partial differential equations (PDEs) coupled with the source term. In comparison, the concentration is described by a single convection-diffusion equation. In terms of their distribution, the range of magnitude of the turbulence quantities is much larger than that of the concentration field, which makes them numerically more difficult to deal with. Besides, the mean flow field is much more complex with a higher Reynolds number. The existence of a thinner boundary layer indicates more rapidly varying velocities spatially, which can also bring challenges to the neural network model.

Contribution of present work
The present work has two major contributions. First, it extends the previous work on vector-cloud neural networks [15] from laminar flows to turbulent flows. Instead of predicting the hypothetical passive concentration field, this work directly emulates the transport equations for turbulent kinetic energy (k) and dissipation rate (ε), which are of physical significance in turbulence modelling. From the perspective of machine learning, the distributions of turbulence quantities are more challenging due to their large range (for example, in the cases considered in this work, the dissipation rate ε ranges from 10 −4 to 2). Special transformation and normalization of the output variables were applied to resolve such a problem. Second, we evaluate the trained neural network model in a coupled setup, which is rarely studied in previous works.
Specifically, we couple the neural network-emulated turbulence model with the RANS equation solver to investigate its robustness and stability.
In this work, the neural network-based turbulence model is trained using the data generated from the k-ε model [30] and is subsequently evaluated against it. Its accuracy and applicability are thus restricted by the training set itself. We emphasize that the objective here is neither to replace nor to surpass the k-ε model with neural networks. Rather, this work serves as a proof of concept regarding the applicability of the VCNN to turbulence modelling problems. By combining present work with the VCNN architecture with equivariance [26], the framework can be further extended to predict Reynolds stress tensor anisotropy.
Potential advantages of the integrated framework over traditional models are, for instance, the flexibility in incorporating training data, and less burden of modelling unclosed terms such as the pressure-strain-rate tensor.
The remaining sections of the paper are organized as follows. Section 2 describes the problem to be solved and the methodology used, including the selected input data matrix, the neural network architecture, and the coupling procedure. Section 3 presents the main results obtained in numerical simulations on periodic hill flow cases with discussions. Section 4 concludes the paper.

Problem statement
In RANS simulations, the flow field is described by the mean flow equations and a closure model. For incompressible flows of constant density ρ, the mean flow equation is as follows: in whichD Dt := ∂ ∂t + U · ∇ is the mean substantial derivative, U j is the mean velocity, u i and u j are the fluctuating velocities, p is the mean pressure, and ν is the kinematic viscosity. The unclosed term in Eq. (1) u i u j (the covariance of fluctuating velocities) is also referred to as the Reynolds stress tensor. In the k-ε turbulence model, the unclosed Reynolds stress tensor follows the Boussinesq hypothesis: where δ ij is the Kronecker delta and ν T is the eddy viscosity. The eddy viscosity is defined by where C D = 0.09 is a model constant. The turbulent kinetic energy k and turbulent kinetic energy dissipation rate ε are described by the transport equations: in which σ k = 1.00, σ ε = 1.30, C ε1 = 1.44, and C ε2 = 1.92, and P = 2ν T S ij S ij is the production of k with S ij the strain rate tensor.
In this work, the goal is to develop a neural network-based turbulence model that plays the role of transport equations (4) and (5)  model. Therefore, as with the linear or nonlinear eddy viscosity models, the weak equilibrium assumption is also made in the neural network-emulated k-ε model [1,31], since only the non-equilibrium in the magnitude of the Reynolds stress is accounted for, and the Reynolds stress anisotropy is still modelled based on local strain rate as shown in Eq. (2).

Data generation and processing
The generation and processing of data follow four steps: (1) perform simulations with the RANS solver with the k-ε turbulence model, (2) determine the neighbourhood on which the central point is dependent, (3) calculate and normalize features and labels, (4) stack stencils for all the sample points for training or testing the neural network.
They are described in detail in the rest of the section.

Flow simulation
The flow is simulated in a channel with periodically appearing hills on the bottom [32]. A family of such periodic hill geometries are parameterized by a slop parameter α as shown in Figure 2 [33]. The slope parameter α is defined as w 1.93H where H is the height of the hills and w is the width of the hills. It is clear that the smaller the slope parameter, the steeper the bottom of hills. L y remains the same while L x is larger for larger slope parameters. A fragment starting with a downslope following a flat region and ending with an upslope is extracted from the periodically varying geometry as the simulation domain. The inlet and outlet are considered as periodic boundaries and the top and bottom of the channel are non-slip walls. The testing dataset comprises 10 extrapolation flow cases and 10 interpolation flow cases. Half of the extrapolation cases have slope parameters larger than 2 and another half have slope parameters smaller than 1. The interpolation test cases are between two neighbouring training cases. The composition of the training and testing set are summarized in Table 1.

Influence region
According to the nonlocality of the turbulence physics and the nature of the transport equations for k and ε, each point in the flow is affected by the entire history of the fluid particle, as well as pressure-coupling with other regions. The VCNN incorporates part of the nonlocality by considering the stencil containing the neighbouring points of the point of interest. It is presumed that the stencil has the shape of an ellipse with its center located at the central point and its major axis aligned to the mean velocity at that point. The idea is that the diffusion is isotropic whilst the advection is strongest along the local flow velocity direction. The lengths of the semi-major and semi-minor axes are determined according to the following equations [20]: in which ν l = 0.1, ζ l = 3.0 are constants related to the flow, and = 0.2 is the tolerance of approximation.
The first two variables ν and ζ correspond to the diffusion and dissipation coefficient in the 1D convectiondiffusion-reaction equation, from which the original derivation of Eq. (6) is made. We assume that such a simplified calculation is sufficient for our purpose and it is later proved by the parametric study of the size of the influence region in Section 3.
The orientation, shape, and size of the stencil vary according to the local mean velocity based on Eq. (6).
A visualization of the stencils at various locations for the case α = 1 is presented in  Since both the size of the stencil and the mesh density vary from location to location, the number of points in a stencil also varies. Specifically, stencils in the middle of the channel are larger while the mesh there is sparser. In contrast, stencils near walls are smaller while the mesh there is denser. By considering the stencil as a set, VCNN has the flexibility of processing stencils with a varying number of points. We took advantage of this property and adopted different strategies for training and validation. For the convenience of training, a random sampling procedure is applied to the training data set, in which 200 points are randomly drawn from the stencil at each mesh point. For the best prediction performance, we use full stencils for validation, meaning that all the data points within the stencils are used as the neural network inputs.

Input features and normalization
In this work, we select 11 features that are relevant to predicting turbulence quantities to build the in which n is the stencil size (number of sampled points in the stencil, for training dataset n = 200).
It is worth noting that when building the input dataset, we have already taken the numerical requirements of neural networks into consideration. For instance, 0 , r and r are introduced such that the feature value is suitable for neural network training (approximately within the range [0, 1]). However, one of the features (s) and both of the labels (k and ε) are not normalized yet. In this work, we additionally cap the strain rate magnitude to 3 for simplicity. The turbulent kinetic energy is normalized by the reference velocity k ∝ ko (k o represents the original turbulent kinetic energy in the flow simulation and k is the normalized value to be predicted by neural networks) and scaled by turbulence intensity I t = 0.2 such that the normalized value falls roughly between 0 and 1. The normalized k follows The normalization of k involves the reference velocity U ref . In our case it can be regarded as the relative velocity with regards to the channel and thus it is Galilean invariant.
The turbulent kinetic energy dissipation rate is transformed by taking the logarithm ε ∝ log 10 ε o considering its distribution in a broad range and its positivity: in which ε o represents the original turbulence kinetic energy dissipation rate and ε is the normalized dissipation rate to be predicted by neural networks. It is further scaled by a constant 0.5 to attain the proper range that is suitable for training.

Neural network architecture and training
The VCNN architecture is composed of two multilayer perceptron (MLP) sub-networks. An overall schematic of the network architecture is provided in Figure 5 and major steps are numbered. The first sub-network is referred to as the embedding network, which aims at extracting higher-level features that are correlated to the turbulence quantities (step 1 ○). The embedding network has a m-dimensional output and it represents m one-dimensional embedding functions φ 1 , . . . , φ m . These embedding functions are evaluated at each cloud center i using the frame-independent scalar features c i , i = 1, 2, . . . , n as inputs and provide G ij = φ j (c i ), j = 1, 2, . . . , m. Through training, the embedding functions will embody important nonlinear As a nonlocal model, the information contained in the single stencil points needs to be integrated into features of the entire stencil. This is achieved through averaging the input features q i according to their embedded weighs G ij (step 3 ○): Eq. (10) can be regarded as m forms of weighted summation of the input features. The operation is permutation invariant since the weights are dependent of single stencil points G ij = φ j (c i ) as emphasized earlier.
It follows directly: Each of the elements in L is a summation over functions of each sample in the stencil. Switching orders of any two of the samples will not influence the results.
Upon introducing q i (thus x and u) into the calculation of L, rotation invariance is violated. This can be solved by applying a pairwise projection (step 4 ○): The pairwise projection QQ means projecting each feature vector to every other feature vector. In this way, the resultant matrix captures the correlation among samples that is independent of any rotation of the reference frame. As both the averaging based on the embedded weights and the pairwise projection are linear transformations, the associative law of linear transformation applies. The matrix D is both permutation invariant and rotation invariant. As the input of the fitting network, D also guarantees the desired invariant properties of the network. In the fitting network, the model's final outputs, the turbulence quantities k and ε are predicted (step 5 ○).
In our implementation, a submatrix G * ∈ R n×m of the full embedding matrix G ∈ R n×m is used in the linear transformation in order to reduce the dimension of D and save computational cost while all the invariance properties are still preserved: We choose m = 4 so that G * contains the first 4 columns of G. The exact hyperparameters of the network used in this work are shown in Figure 6. Besides, the ultimate target turbulence viscosity ν T in Eq. (3) is also included as a supplement term for guiding the network towards better prediction performance. The complete loss function then follows where λ is a weighting factor balancing the two normalized error terms of k and ε with the error term of ν T ; k,ε andν T represent predictions made by the neural network model. Adam optimization algorithm [36] is used for training neural network parameters with an initial learning rate 10 −3 . The learning rate decreases by a factor of 0.7 for every 600 steps. The network is trained until the training error reaches the lowest and the validation error stays at a low level. The number of epochs is around 2000. The construction of the neural network and its training are implemented using the PyTorch machine learning library [37].

Coupling with RANS equation solver
The neural network described in the previous section should work as a usual turbulence model in flow simulation. In this work, the neural network-based turbulence model is coupled with the SimpleFOAM solver from the OpenFOAM package. The network requests mean flow solutions (mean velocities u and mean strain rate s), geometry data (relative coordinates x , cell volume θ, boundary and wall distance b and η, proximity to cloud center r and r ) from the OpenFOAM solver and returns turbulence quantities k, ε.
The pipeline of data exchange is presented in Figure 7 and the steps for a coupled solver are as follows: (1) The mean flow field and the turbulence quantities are initialized by the initial condition. The absolute wall distance η o of all the mesh points is computed.
(2) SimpleFOAM solver reads the initial data or results from the last time step and solves for the mean flow quantities, e.g., flow field, pressure. Data required by the neural network such as the strain rate magnitude, cell volume, etc., are also calculated. (3) The data is read and processed as described in Section 2.2, including collecting the stencil point, calculating features and normalization. The processed data is put into the input matrix Q and fed into the trained neural network.
(4) The data is processed first through the embedding network, then transformed to achieve the permutation invariance and frame-independence. Finally, the prediction on turbulence quantities is given by the fitting network and is written into OpenFOAM format.

Performance in uncoupled and coupled setup
We test the performance of the trained neural network model in both uncoupled and coupled setups.
Specifically, in the uncoupled setup, the neural network model uses the converged mean flow field obtained from traditional flow simulation to predict the turbulence quantities, which is similar to how the network is trained. In a coupled setup, the neural network model is coupled with a traditional RANS equation solver.
The coupled solver starts iteration from a given initial condition, that is, an unconverged flow field. During the simulation, the neural network model makes predictions on turbulence quantities in each iteration. The simulation runs until it reaches convergence. In both cases, the neural network prediction is compared with that of the standard k-ε turbulence model. It was found that, in the uncoupled setup, the neural network model performed very well in interpolation flow cases, but the prediction worsened drastically in extrapolation cases. In the coupled setup, the network model performance was comparably good as it was in the uncoupled setup.

Uncoupled results
In the interpolation testing case α = 1.45, the prediction of the neural network is highly consistent with the results of k-ε model. Although the flow case α = 1.45 is not included in the training set, the contour of ν T with k-ε ( Figure 8a) and that with neural network model (Figure 8b) are nearly identical in terms of their pattern and magnitude. When comparing their cross-section profiles (Figure 8c) along the channel, it is also clear that the two profiles overlapped each other. The neural network prediction is only inferior to that of k-ε model in terms of its smoothness. The neural network model performed less satisfactorily in extrapolation cases. The performance deteriorated rapidly in the testing cases with slopes further from the training set. Simulation results in Figure 9 show the prediction capacity of the neural network model on mild extrapolation cases. In the flow cases α = 0.9 and α = 2.1, the network prediction is close to that of the reference results provided by the k-ε For a more quantitative analysis of the model performance, an overall metric of validation (prediction) error rate is introduced: whereν T is the neural network prediction, and ν T refers to the results of k-ε model (regarded as the ground truth). The errors of all mesh points are summed and the total error is normalized by the sum of squared ground truth results.
The neural network model shows good performance in the interpolation regime but in the extrapolation regime, the performance is less satisfactory. This finding is clearly shown in Figure 10, where the prediction error rate is plotted against the slope parameters. Within the training set range, the prediction error of the neural network model is as low as the training error (around 3%). The prediction error rate grows rapidly when the slope parameter is reduced from 1.0 to 0.5. Similarly, the error rate increases when the slope parameter goes beyond 2.0, although the error rate of the case α = 2.5 is slightly lower than that of the case α = 0.5.

Coupled results
In the coupled setup, the neural network model performs comparably well as in the uncoupled setup in terms of the model accuracy. The converged fields in one interpolation and one mild extrapolation flow case of the coupled solver with the neural network-based turbulence model and the k-ε model are presented in Figure 11. The cross-section profiles of the eddy viscosity ν T and the mean velocity in x and y direction (u and v) are shown. We note that the difference between the neural network-based turbulence model and the baseline k-ε model is very small. Especially in terms of the mean flow field (the final results of RANS simulation), the difference is barely discernible.
More importantly, the network has good stability and the coupled solver is capable of reaching a steady state when initialized with an unconverged flow field at the early stage. The initial condition of the simulation is shown in grey dashed lines in Figure 11, and corresponds to the baseline simulation at the 500 th iteration.
There is a large discrepancy between the initial condition and the final converged results, especially in terms of the eddy viscosity. Such an unconverged flow field is not included in any training cases. Even in an unseen flow field, the neural network-based turbulence model is capable of providing reasonable prediction without causing divergence of the primary RANS equation and brings the simulation to a convergence that is very close to the k-ε model.
The starting point of 500 th iteration is found through grid search, i.e., we tested the neural networkcoupled RANS solver starting from the converged flow field all the way down to the 1 st iteration. Tests starting from iterations earlier than the 500 th hardly succeeded. This is because the difference between the mean flow fields in the initial stage and those in the final stage is too large, and consequently the network is not capable of giving valid predictions. Receiving incorrect turbulence quantities, the simulation diverges very soon. How to achieve even better stability of the neural network-based turbulence models remains an important question to be further investigated.

Scale of the influence region
The parametric study on the scale of the influence region shows a rapidly diminishing marginal benefit of using a large influence region in our case. In Figure 12, the validation error is plotted against the coefficient describing the relative size of the influence region as compared to the baseline. This coefficient follows , in which l 1 is the test semi-major axis of the influence region, and l 10 is the original semi-major axis of the influence region based on Eq. (6). The sem-minor axis l 2 is scaled by the same constant C I . Therefore, the ratio between the area of the test and baseline influence region is I/I 0 = C 2 I . The validation error shows an overall declining tendency with regards to the enlarging influence region.
However, there are clearly two stages that can be observed. The first interval is from C I = 0 to around C I = 0.1, in which the improvement of performance is the most evident. In the second interval (C I > 0.1), the performance stays roughly at the same level.
As indicated in Figure 13, when coefficient C I reaches 0.1, the stencil in most parts of the flow domain contains neighbouring points both along and orthogonal to the streamwise direction. By including these critical points, the network performance is greatly enhanced. Influence regions larger than this can be considered as containing the most important set of information for predicting turbulence quantities. It is noticed that in the recirculation region, C I = 0.2 can also be filed under the first interval. Because in the near-wall region, only until C I = 0.2 are the neighbouring points along the main flow direction included.
In the second interval, by including points in a larger neighbourhood, the prediction capacity of the network increases only by a small amount. There is possibly even a decline of performance in flow case α = 0.5 when the influence region increases from C I = 2.0 to C I = 3.0. This may be due to the fact that noise consisting of points too far away is introduced and undermines the prediction ability.

Flow physics learned by the neural network
The neural network model has strong boundary dependence as shown in the visualization of one of the embedding functions. We set m = 1 for a better interpretation of the network. In this way, only the first embedding function is multiplied twice in the linear transformation step and it is attached with more importance than other embedding functions. The embedded weights G i1 defined by the embedding function φ 1 where i denotes the index of spatial points is visualized in Figure 14. Across stencils at different locations, the contours of the weights showed a consistent pattern. The gradient of the embedded weights is well aligned with the wall-normal direction, which is especially evident for stencils near the periodic hills. The pattern of weights agreed well with the physics that the flow inside a periodic hill channel is dominated by the wall.

Conclusion
In this work, we apply vector-cloud neural network [20] to address the long-standing turbulence modelling problem. It is a network architecture incorporating the nonlocality of turbulence quantities and guaranteeing the invariance properties of turbulent flows. This work is a proof of concept for the application of VCNN to turbulence modelling in a couple setting and pave the way for the framework to predict Reynolds stress tensor anisotropy.
The neural network-based turbulence model was embedded into the RANS solver and has demonstrated good stability in the coupled setup. We tested the stability and robustness of the coupled solver by applying it at different stages (e.g., 500 th iteration, 100 th iteration) of the simulation. The neural network coupled solver could bring the simulation to a convergence (with good accuracy) starting from an early stage (500 th iteration) as shown in our testing cases. The stability of neural network models is usually not guaranteed in physical simulations. Thus, our numerical test results demonstrate the promising future of VCNN in other application scenarios.
When adapting the network from laminar to turbulent flows, the wide distribution of turbulence quantities k and ε brings about challenges for machine learning. We apply a transformation to each of the inputs and outputs. It is shown in numerical tests that the training efficiency has largely improved after transformation.
Besides, based on the parametric study on the scale of the influence region, there is a diminishing marginal utility when increasing the influence region. Whether this discovery is universal to other flow cases requires further investigation. Finally, the physics learnt by the network showed a strong boundary relevance which was consistent with the expected flow physics for channel flows. It suggested that the vector-cloud neural network was able to learn the correct nonlocal physics from the provided data.