Xinlei Zhang
email: zhang.xinlei@hotmail.com

Thomas Gomez

Olivier Coutier-Delgosha

Bayesian optimisation of RANS simulation with ensemble-based variational method in convergent-divergent channel

Keywords: RANS, Bayesian optimisation, EnVar, Venturi, Bump, inlet velocity inference, model correction inference

This paper investigates the applicability of a hybrid data assimilation approach, namely ensemble-based variational method (EnVar), to optimise Reynolds Averaged Navier-Stokes (RANS) simulations in convergent-divergent channel from the perspective of Bayesian inference. Concretely, the ensemble-based variational method is applied to infer the inlet velocity and turbulence model corrections by assimilating Direct Numerical Simulation (DNS) results or limited experimental data. The approach is first adopted to infer the inlet velocity profile for the WallTurb Bump and Venturi geometry. The improvement can be achieved near the inlet region for the bump, but for Venturi in light of the view field limited in adverse pressure gradient region, the observation space is not sensitive to the perturbation of inlet condition. In a second step, the model corrections in k -ω SST model are investigated by assimilating the limited sparse experimental data. With the inferred model corrections, the predictions on both velocity and turbulent kinetic energy (TKE) get improved. The results indicate that the ensemble-based variational method is efficient in inferring unknown quantities of both low dimension (D = 20) and high dimension (D = 2400) with small ensemble size robustly and non-intrusively. This approach could prove very useful for Bayesian inference or optimisation in CFD problems.

Introduction

Despite significant development of Computational Fluid Dynamics (CFD) for several decades, the high-fidelity resolution Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) is still computationally intractable for most applications, especially with high Reynolds numbers. Reynolds-Averaged Navier-Stokes (RANS) simulation will remain dominant for industrial applications in the near future. However, it has been noted that RANS model cannot make accurate predictions for the turbulent flows in the convergent-divergent channel where there exist mean curvature and pressure gradient, due to the ambiguous boundary condition and the inadequate turbulence model. On the other hand, experimental investigations have to be faced with the challenges of the sparse measurements in the limited observable field, noise contamination and insufficient resolution for the small scale flow. To address these issues, data assimilation (DA) approach widely used in the oceanography and geography gains the spotlights over the past few years in the turbulence community. Based on Bayes theorem, data assimilation can integrate the low-fidelity RANS calculation with the high-fidelity resolution from LES/DNS or sparse experimental measurements to infer the unknown boundary condition or underlying model information and thus optimise the RANS predictions.

Data assimilation can be sorted by variational data assimilation method, Ensemble Kalman Filter (EnKF) method, and hybrid methods. All these methods leverage the Maximum a Posteriori (MAP) estimation, while in the variational method, MAP is formulated as the minimisation of the cost function through adjoint optimal least square techniques. A body of works based on the variational method has demonstrated its inferential performance and robustness to replicate the flow status. Gronskis et al. [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] established a data assimilation framework based on the variational method incorporating with adjoint optimisation method to generate the inflows condition for DNS. Foures et al. [START_REF] Dpg Foures | A data-assimilation method for reynolds-averaged Navier-Stokesdriven mean flow reconstruction[END_REF] applied the variational method to minimise the discrepancy between the time-averaged velocity fields of a DNS resolution and an incompressible RANS simulation for the two-dimensional flow past a circular cylinder at a Reynolds number of Re = 150. Symon et al. [START_REF] Symon | Data assimilation of mean velocity from 2D piv measurements of flow over an idealized airfoil[END_REF] applied the variational method to reconstruct the flows with relatively high Reynolds number of 13,500 around an idealised airfoil by assimilating the mean velocity field from time-averaged Particle Image Velocimetry (PIV) measurements.

Also, the EnKF method where the state statistics are updated with an ensemble of realisations, have been intensively applied to quantify and reduce the uncertainty in RANS simulation. Colburn et al. [START_REF] Colburn | State estimation in wall-bounded flow systems. part 3. The ensemble Kalman filter[END_REF] used EnKF to estimate the near-wall turbulent flow based on the wall information from DNS resolution. Kato and Obayashi [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] applied EnKF to infer the optimal parameters in the Spalart-Allmaras turbulence model for zero-pressure gradient flat plate boundary layer at Mach number of 0.2 and Reynolds number of 5 × 10 6 . They [START_REF] Kato | A data assimilation methodology for reconstructing turbulent flows around aircraft[END_REF] also used the Ensemble Transform Kalman Filter to integrate the CFD and experimental fluid dynamics (EFD) to replicate the transonic turbulent flows over RAE 2822 airfoil and ONERA M6 wing through estimating the proper angle of attack, Mach number, and turbulent viscosity. Heng et al. [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed bayesian approach[END_REF] introduced uncertainty in Reynolds stress directly and adopted an iterative ensemble Kalman method to reduce the model-form uncertainty in kmodel for the flow over periodic hills and the flow in a square duct by assimilating very sparse observations. Despite both the variational data assimilation method and EnKF have demonstrated their applicability for the Bayesian optimisation of CFD problems, the intrusiveness of variational method makes it difficult to implement, and the high sensitivity to the prior statistics for EnKF causes the inference results prone to be inferior. Recently, a hybrid method, namely ensemble-based variational method(EnVar), gains attention to solve the inverse problem for CFD simulation with both robustness and non-intrusiveness. Mons et al. [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF] first explored the applicability of this kind of data assimilation technique into unsteady flows with coherent gust and assessed the robustness of variational data assimilation method, Ensemble Kalman Smoother and EnVar method to reconstruct the flows around a cylinder. They demonstrated that the ensemble-based variational method could be robust as the variational method and circumvent the efforts on the adjoint model by estimating the prior statistics with ensemble realisations. Meanwhile, the inferential performance can be very satisfactory and is not quite sensitive to prior and observation statistics comparing to Ensemble Kalman Smoother.

The current study investigates the applicability of the ensemble-based variational method to optimise the RANS simulation from the perspective of inferring improved inlet boundary condition and underlying model corrections by incorporating with DNS resolution or sparse experimental observation. The test cases from low dimension (D = 20) to high dimension (D = 2400) of input parameters all demonstrated the merits of the proposed approach. It is worth noting that the limited observation from the PIV experimental measurements is also integrated with RANS simulation to recover the flow status.

The rest of this paper is structured as follows. The ensemble-based variational scheme is presented in Section 2. The CFD code and the practical implementation of the data assimilation framework are described in Section 3. The applications of this framework to infer inlet velocity and model corrections are presented in Sections 4 and 5 respectively. Section 6 is dedicated to the conclusion and perspectives.

Data assimilation framework

Ensemble-based variational scheme

The ensemble-based variational method is a hybrid data assimilation approach that combines the variational data assimilation method with EnKF. Compared to the variational method, the formulation of EnVar uses a Monte Carlo ensemble to estimate the prior statistics thus circumventing the efforts on the adjoint operator. This method equates MAP to the optimal control problem where the optimal control vector is searched based on gradient-decent optimisation to update prior information. Meanwhile, the prior distributions are re-estimated at each DA iteration thus counteracting the influence from the improper prior, while in EnKF scheme the Kalman gain matrix to update the prior distribution, is directly constructed from the approximated prior statistics.

The state vector is the input parameters, representing the quantities to be inferred. For simplification, the vector is assumed to be Gaussian distributed. The mean of this normal distribution corresponds to the initial guessed or prior state vector, and the vector can be expressed as

α = α (e) + E β, ( 1
)
where β is the control vector with dimension of N en . N en is the ensemble size. E = (α (1) α (e) , α (2) -α (e) , . . . , α (N en ) -α (e) )

Regarding the observation, it is the reference data (e.g. friction/pressure coefficient, velocity fields, and so on) from high-fidelity simulation or experiment that have a substantial influence on the inferred quantities. Random measurement noise in observation is assumed to be Gaussian distributed with zero mean, uncorrelated, and are characterised by the relative standard deviation σ o . The observation y can be prescribed as

y = h(α + ) + , ( 2
)
where h is the observation operator that maps the state space to observation space, α + is reference trajectory projecting to the observation and is the possible random measurement error.

Based on the Bayes theorem, to maximise the posterior is equivalent to minimise the cost function J

p(α | y) ∝ p(α)p(y | α) ∝ e -J → J = 1 2 α -α (e) 2 B -1 + 1 2 d -y 2 C -1 , ( 3
)
where d is the RANS realisations in observation space based on state vector α, y is the observation from DNS solutions or the experimental measurements, * 2

B -1 = * T B -1 * ; * 2 C -1 = * T C -1 * .
As in (3), the cost function J is composed of two parts: the difference between the prior and its realisations, and the difference between RANS prediction in observed quantities and the observation, weighed by the background covariance B and observation error covariance C, respectively. d can be realised by the linearisations as in the following formula: (1) )h(α (e) ), h(α (2) )h(α (e) ), . . . , h(α (N en ) -α (e) )).

d = h(α (e) ) + H β, H = (h(α
(

) 4
Using B (1/(N en -1))E E T and substituting α with (1), the cost function in (3) can be rewritten as a quadratic equation on control vector β:

J = 1 2 (N en -1)ββ T + 1 2 h(α (e) ) + H β -y 2 C -1 . ( 5 )
The gradient and hessian of J can be derived straightforwardly as

∂J ∂β = (N en -1)β + H T C -1 (h(α (e) ) + H β -y), ( 6
)
∂J 2 ∂ 2 β = (N en -1)I + H T C -1 H . ( 7
)
In order to minimise the cost function [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF], one iteration of Newton CG method is performed with ( 6) and [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed bayesian approach[END_REF]. The obtained β is used to update the state vector α according to [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. This iterative process is continued until the converge criterion is reached. The procedure of the ensemble-based variational method is shown in Figure 1 and the details can be summarised as follows:

Step 1. Give a first guessed or prior state vector α (e) , and prescribe the prior and observation statistics respectively;

Step 2. Realise N en samples around the mean vector α (e) ;

Step 3. Propagate each sample of state vector to observation space through solving RANS equation;

Step 4. Analyze the control vector β by minimising cost function with (6) and (7);

Step 5. Update the mean of state vector with the analyzed control vector β based on (1); Step 6. Return to Step 2 and resample around the updated ensemble mean until the fitting criterion or maximum iteration is reached.

Fitting criterion

The goal of the data assimilation approach is to best fit the numerical prediction with DNS or experiments. Hence in order to evaluate the optimal performance, it is first necessary to calculate the ratio of cost function value before and after the data assimilation process in order to assess the extent to which the cost function is reduced. The ratio of cost function J can be defined as:

r J = J end J 0 . ( 8
)
Moreover, the cost function is constituted of two parts as illustrated in the formula [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] and thus reduction of the cost function may not directly reflect the decrease of the discrepancy between the numerical prediction and the observation in the quantities of interest. Therefore, the formula ( 9) is introduced to measure this distance:

J ob = 1 2 y -h(α (e) ) 2 , ( 9
)
where * is Euclidean norm. In this work, the data assimilation process is terminated as J ob < 10 -3 . The ratio of the assimilated and initial J ob is also introduced to evaluate the efficiency of the optimisation process

r J ob = J ob end J ob 0 . ( 10
)

Numerical setup

CFD solver

A 2D steady incompressible solver is used to perform the numerical simulations on twoequation RANS turbulence model. The governing equations can be expressed as:

div(F c -F v ) = S, (11)
with

F c = ⎛ ⎝ ρu ρv ρu 2 + p ρuv ρuv ρv 2 + p ⎞ ⎠ and F v = ⎛ ⎝ 0 0 2μS xx + τ xx 2μS xy + τ xy 2μS xy + τ xy 2μS yy + τ yy ⎞ ⎠ ,
where ρ is density; u and v is stream-wise and wall-normal mean velocity; F c and F v denote the convective and viscous flux densities; μ, S and τ represent molecule viscosity, mean strain rate, and Reynolds stress respectively; S is the source term.

The Reynolds stress term τ is modelled by Menter's k -ω SST model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. And the code uses the HLPA (Hybrid Linear / Parabolic Approximation) non-oscillatory second-order scheme for the convective term. The SIMPLE (Semi-implicit Method for Pressure-Linked Equations) algorithm is applied to solve the coupled mass and momentum conservation equations, on two-dimensional structured curvilinear-orthogonal meshes.

The non-slip condition is used on the wall. The first grid for all cases in this paper is well placed in the viscous layer and the wall function is not incorporated therefore to eliminate the effects of assumptive wall boundary condition.

Data assimilation implementation

Firstly, an initial guess of the input parameter is given as prior, and then N en samples are constructed around the first guess based on Gaussian process. To ensure the smoothness of the obtained samples, a non-diagonal covariance matrix to describe the prior distribution is prescribed as:

i,j = σ (x i )σ (x j )b 2 exp - x i -x j 2 l 2 c , (12)
where the variance σ (x) is constructed based on the discrepancy between the prediction from initial RANS simulation and the observation in order to inform where large or small perturbations are expected; b refers to the characteristic values, for inlet velocity reconstruction it is based on inlet bulk velocity while for model correction inference, it is taken as 1; l c is constant correction length based on the height of the channel at inlet; x i refers to the position of the ith control volume. With mean vector α (e) and covariance matrix i,j , the N en samples around the first guess could be realised. The observation error is assumed to be the uncorrelated Gaussian distribution with zero mean. For simplification, the covariance matrix C is constructed as a diagonal matrix based on the observation error which is defined as T . The resolution in the observation space is mapped from the state vector forwardly with CFD solver. The posterior is obtained by solving the Bayesian optimisation problem based on the ensemble-based variational scheme. For the optimisation method, Newton-CG is used to minimise the cost function.

Inlet velocity inference

The proper inlet velocity has critical implications for the performance of numerical prediction, and many forward methods have been proposed to generate the inflow condition for DNS and LES, such as the recycling-based method and synthetic turbulence generators [START_REF] Wu | Inflow turbulence generation methods[END_REF], whereas high computational efforts are required in these methods. In this section, the ensemble-based variational method is explored to infer the inlet velocity profile based on the limited quantities from DNS resolution or experimental measurements. Two test cases are used to evaluate the performance of the data assimilation scheme on the reconstruction of the inlet velocity profile. The first one is turbulent flow in the WallTurb Bump where DNS dataset is available, and the other one is non-cavitating flow in Venturi where the experimental X-ray measurements can be considered as the reference data. A summary of data assimilation experiments for inlet velocity inference is given in Table 1.

WallTurb Bump

Flow configuration

The turbulent flow in WallTurb Bump [START_REF] Marquillie | Direct numerical simulation of a separated channel flow with a smooth profile[END_REF] is one canonical flow and widely used to verify the performance of numerical methods [START_REF] Jesus | Large eddy simulations and rans computations of adverse pressure gradient flows[END_REF][START_REF] Kuban | Les modeling of converging-diverging turbulent channel flow[END_REF][START_REF] Laval | Large-eddy simulations of control of a separated flow over a 2D bump by means of pulsed jets[END_REF]. The Reynolds number for this flow is 12,600. The computational domain is -5.22 < x < 7.34; 0 < y < 2. In order to simulate this flow, the structured curvilinear-orthogonal mesh is generated with 125 cells in streamwise direction and 60 cells in the normal-to-wall direction. The y + of first mesh adjacent to the wall is around 1. The mesh of the WallTurb Bump is shown in Figure 2.

For the setup of data assimilation, the input parameters are the lower half of inlet velocity with the dimension of 30, and the ensemble size N en is set to be consistent with the dimension of input space. The variance field σ is based on the discrepancy field of the prior inlet velocity and the DNS velocity with a multiplication of 10 -3 . The value of multiplication need some trials and errors since large variance will increase the spatial extent of subspace where the optimal solution is searched, and thus result in the instability of the optimisation process, while too small variance will make the process robust but quite slow.

In order to infer the inlet velocity profile for this case, the friction coefficient C f is regarded as the observation in consideration of the strong correlation between the velocity and the friction coefficient. The skin friction coefficient is defined by:

C f = τ w 0.5ρU 2 ref
where τ w is the wall shear stress τ w = μ(du/dy)| y=0 and U ref is reference velocity, which is taken as the dimensionless inlet bulk velocity 1. Since presumably the friction coefficient is not sensitive to the inlet condition in the region with the bump where the flow encounters the strong pressure gradient, the observation is confined near the inlet region and exclude the part adjacent to the inlet. Thus, the observation includes the friction coefficient on the bottom wall from x = -4.2 to -1.9, and the standard deviation of observation error is taken as 5 × 10 -4 . The Newton-CG method is applied to minimise the cost function. As in this case the ensemble size is quite small, the Hessian of the cost function can be explicitly expressed, and the Newton-CG is more robust compared to the quasi-Newton method with approximated Hessian matrix such as BFGS. The comparison of Newton CG and implicit BFGS method to infer inlet velocity for the bump is presented in Appendix 1.

Physical constraints

Due to the non-uniqueness of the optimal solution, data assimilation process may result in the inferred inlet velocity losing physical meaning, for instance too high velocity adjacent to the wall and non-symmetric inlet velocity profile for this case, which may also bring about the divergence problem for CFD solver. Therefore, three constraints are given based on physical knowledge. Firstly, the dimensionless velocity at the first grid to the bottom wall is fixed at 2 × 10 -3 via interpolation from DNS resolution to avoid the steep velocity gradient near the boundary and have a reasonable y + at the inlet. Secondly, the reconstructed inlet velocity profile may lead to the variation of flux and accordingly change the flow condition. Therefore, the flux at the inlet is corrected after each data assimilation iteration, by multiplying the ratio of the updated flux and reference flux to ensure the flux constant. Besides, in this case, the inflow into the channel should be developed turbulent flow and presumably symmetric. On the other hand, only the friction coefficient on the bottom wall is considered as the observation; hence the inlet velocity near the top wall is not able to be recovered due to lack of information. Thus only the bottom half profile is taken as the input parameter, while another half is constructed by symmetric projection.

Results

Two different first guess is imposed as the prior inlet velocity: the parabolic and the flat profile. The assimilated results of inlet velocity and observed friction coefficient are presented in Figure 3 with the comparison to the prior and DNS. The inferred inlet velocity profiles are both quite close to the DNS data as shown in Figure 3(a,c). Accordingly, the predictions in terms of friction coefficient presented in Figure 3(b,d) also have a good agreement with DNS data in the observed region. Generally, both cases can obtain a good inference on inlet velocity and improve the prediction on C f by assimilating the friction coefficient from DNS.

Figure 4 shows the evolution of cost function to the iteration number. For both cases, the cost function is reduced significantly while with the prior of parabolic profile the minimisation is more efficient and converged in the first five iterations. In concrete, the cost function J reduces from 32,386 to 80 in the case with parabolic prior, and the ratio r J ≈ 0.0025, while for the prior with flat inlet profile, the cost function J reduces from 57,607 to 114 and r J ≈ 0.002. The norm of discrepancy between numerical prediction and reference with first guessed parabolic inlet velocity is reduced from 0.0155 to 0.0031 and r J ob ≈ 0.20 while the discrepancy for the case with the flat inlet velocity profile can be reduced from 0.0191 to 0.0030 and r J ob ≈ 0.157.

The contour plots of the velocity U field with first guessed parabolic velocity are presented in Figure 5. The visible improvements can be observed near the inlet and the favourable pressure gradient region. However, in the adverse pressure gradient region, the prediction is not improved with optimal inlet velocity, which is not surprising since it is a consensus that RANS model cannot give good predictions with confidence when it encountered strong adverse pressure gradient. A large-scale separation after the summit of the bump can be seen in RANS calculation, whereas in the DNS solution the reattachment occurs in the downstream near the starting point of separation and thus there is no noticeable separation. In other words, the prediction of velocity U in the adverse pressure gradient region is insensitive to the inlet velocity, and the poor prediction may be due to the RANS model inadequacy. This section is mainly to explore the applicability of data assimilation to optimise inlet velocity, so the uncertainty in the model is not concerned.

Venturi

The second case is for Venturi channel, which is extensively used in the investigations of turbulent cavitating flows [START_REF] Charrière | Numerical investigation of periodic cavitation shedding in a venturi[END_REF][START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF]. The field information near the throat is quite challenging to be captured by current RANS model due to abrupt curvature change. In this work, the data assimilation approach is applied to infer the proper inlet velocity in non-cavitating flow by assimilating one velocity profile from experimental data.

Flow condition

The experiments were conducted by Khlifa et al. [START_REF] Khlifa | Fast x-ray imaging of cavitating flows[END_REF]. Through applying ultrafast X-Ray imaging into the turbulent flows, the velocity was measured within a Venturi-type test section with 18 • convergence angle and 8 • divergence angle. The cross-section in the entry of the Venturi is a rectangle of 17 mm × 4 mm, and the height of the throat is 15.34 mm. In the experiment, the flux rate is 55.5 L/min, representing the entry bulk velocity at 13.6 m/s. Reynolds number is 2.3 × 10 5 . The time-averaged velocity u is obtained by averaging all the instantaneous stream-wise velocity processed from the high speed photography images. For the numerical setup, the structured mesh is generated with 200 grid in the stream-wise direction and 70 grids in the normal to the wall direction. The y + of the first grid near the wall is ranged from 1 to 3. The mesh of the venturi-type section is shown in Figure 6.

However, since the measured area is in the adverse pressure gradient region where the velocity field is not sensitive to the inlet condition, the data assimilation experiment within Venturi turned out to be a failure. The results are presented in Appendix 2. In view of this, we simplified the Venturi geometry as a divergent channel with the resolution domain starting from the throat of the Venturi section. Thus the measured field is quite close to the inlet, and the flow status can be very sensitive to the inlet condition. The structured mesh is generated with 131 cells in the x-wise direction and 70 cells in the y-wise direction.

Data assimilation setup

In the Venturi-type section, the measured window is placed near the bottom wall and quite small compared to the height of the channel. Therefore, only the inlet velocity profile adjacent to the bottom wall is considered as the input parameters, and the velocity in the other area is fixed at constant 1.1 to keep the flux consistent. The dimension of the input parameters is 20. The velocity at the first grid is also fixed at 2 × 10 -3 to have a reasonable y + .

As for the observation, the velocity field can be captured by the PIV experimental measurements while other turbulent quantities of interest on the wall such as skin friction coefficient are not straightforward to measure. Hence, we regarded the velocity profile at x = 0.0008 m as the observation which is quite close to the inlet boundary, thereby ensuring that the velocity at this specific position can be affected flexibly by inlet condition. The ensemble size, in this case, is set as 20. The variance σ (x) is constructed based on the discrepancy of the RANS simulation and experimental velocity profile at the observed position and the added multiplication is 10 -5 . The standard deviation of observation is σ o = 10 -6 , which represents high confidence in the experimental data since the observed position and the inlet is almost linearly correlated, and thus the inverse problem is well posed for this case.

Results

Figure 7 shows the assimilated and prior results for inlet velocity and the observed velocity profiles. Figure 7(a) presents the evolution of inferred inlet velocity profile every 5 iterations, while Figure 7(b) presents the improved prediction on the velocity profile at x = 0.0008 m. It can be seen that EnVar method can reconstruct the inlet velocity to fit quite well with the reference velocity profile at the specific position near the inlet.

Figure 8 shows that the evolution of cost function with respect to the iterations. The cost function J is reduced significantly from 5.2 × 10 11 to 1.4 × 10 8 after 20 iterations with ratio r J of 2.7 × 10 -4 and r J ob of 2 × 10 -5 .

Figure 9 presents the contour plots of velocity u of prior, posterior and experiment. Noticeable improvements can be seen in the observed position x = 0.0008 m comparing to the prior. However, in the other areas, the flow status is entirely different, since substantially the flows in Venturi cannot be represented with the divergent channel.

Model correction

As illustrated in Section 4.2, in Venturi-type section the measured field is mainly in the adverse pressure gradient region, where the velocity is insensitive to the inlet condition, and the RANS method is incapable of making accurate prediction due to the convex curvature. Therefore, in this section, the RANS model-form uncertainty is considered through the field inversion approach [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]. The conventional k and ω transport equation in k -ω SST model is inadequate especially in the presence of adverse pressure gradient. Thus, the underlying source term is introduced in the k -ω SST model by three different means depending on where the correction terms are inserted (TKE equation or ω equation). The sensitivity of these corrections concerning the observation is analyzed, and accordingly, the correction fields are inferred through solving the inverse problem with EnVar approach. The summary of DA experiments for model correction inference is presented in Table 2.

Correction in k equation

Data assimilation setup

In the TKE transport equation of the k -ω SST model, the sum of turbulent-transport term and pressure diffusion term is represented through Boussinesq assumption in analogy to molecular transport process [START_REF] Wilcox | Turbulence modeling for CFD[END_REF] which probably result in poor prediction for complex flows. The correction variable β c is introduced in front of the production term in TKE transport equation to account for the model uncertainty, which is equivalent to introduce a source term (1 -β c ) × P in the TKE transport equation. Thus with data assimilation, the flow status can be recovered through finding out the optimal correction field to have a good agreement with reference data.

∂k ∂t + ∂(u j k) ∂x j = β c P -β * ωk + ∂ ∂x j (ν + σ k ν t ) ∂k ∂x j ( 13
)
Since the measurement region is only near the throat, in order to reduce the dimension of input space and guarantee that the inferred correction and observation are locally correlated, the range of input parameters are confined in the area with extension to the observed region instead of the whole computational domain. The correction variables are well imposed on the mesh grids. The dimension of input parameters is 1200 with 40 points along the x-wise direction and 30 points along the y-wise direction.

As for the observation, four profiles on the velocity u and TKE are concerned. For each profile, there are 24 measurement points. Given that the numerical solution is imposed on the mesh points while the experimental data is on a uniform Cartesian grid, to compare the numerical resolution with the measurements, each experimental data is interpolated on the numerical mesh points. Besides, since the initial discrepancy of velocity and TKE causes the weight of the two observed quantities on the cost function to be different, TKE is normalised by multiplying the ratio of the initial discrepancy between prediction and reference data in velocity to that in TKE, thereby keeping them in the consistent range.

For the data assimilation setup, as it is not practical to draw all the samples for high dimension case, we drew 50 samples with truncated Singular Value Decomposition(SVD) which can capture more than 99 percent of the variance. The standard deviation of observation σ o is 10 -1 . It is noted that the large observation error represents the experimental instinct uncertainty from both measurements and post-processing [START_REF] Khlifa | Fast x-ray imaging of cavitating flows[END_REF]. Moreover, because the model error is not considered in this work, it can be regarded as the whole process error including the data error as well as the model error. Also, the strong non-linearity for this case will increase the ill-posedness of the optimisation problem, hence for the ensemble-based methods, the inflation on the observation error covariance matrix is commonly adopted to regularise the problem. The large observation error here can be taken as the implicitly inflated one. The variance field as shown in Figure 10 is constructed through interpolation based on the discrepancy between RANS results and very sparse experimental data in TKE where the added correction term has a direct impact. And the multiplication of the variance is taken as 1 × 10 -3 . The first guessed input parameter is set as 1 in the domain of interest.

Results

Figure 11 shows the reduction of cost function values with respect to iterations. It can be seen that the data assimilation process is well converged and the cost function is reduced significantly in the first 20 iterations from 5.1 × 10 5 to around 2.2 × 10 5 with the descending ratio of 0.43, while J ob is decreased with a ratio of around 0.17 at the end of optimisation.

Figure 12 is the comparison of TKE and stream-wise velocity u among the prior, posterior and experimental data along profiles. The profiles of prior TKE are quite distant from the reference profiles. The RANS model cannot capture the high TKE near the wall and throat region. After the data assimilation process, a noticeable improvement can be seen comparing to prior, and the prediction in TKE has a good agreement with experimental data especially at the first two positions as shown in Figure 12(b). However, from Figure 12(a) the prediction on velocity u is not quite improved. The correction embedded in the TKE transport equation may not have substantial effects on the observed velocity, that is to say, the perturbations on the corrections can impact significantly on the prediction of TKE obviously, but have no sufficient influences on velocity u.

The contour plots of posterior TKE and velocity u with correction in the k transport equation are shown in Figure 22(c,d). The apparent improvement can be found in TKE, and the region near the wall with high TKE can be recovered. However, for the velocity, no apparent improvements can be seen, and the relatively low value in the region near the wall cannot be captured even though with the correction in the k equation.

Correction term in ω

Data assimilation setup

The specific dissipation ω transport equation in k -ω SST turbulence model is heavily modelled with an ad-hoc form. The underlying source term in this equation is also probably responsible for the poor predictive performance on the velocity and TKE. In this subsection the correction β c is introduced in ω transport equation as the following formula: The range of input parameters and observations are the same in section 5.1 as well as the ensemble size. The standard deviation of observation is also taken as 10 -1 . The variance σ (x) is taken based on the discrepancy between RANS and experimental data on velocity u since the correction in ω transport equation can have strong correlations with velocity based on our prior study. The multiplication on variance σ (x) is 1 × 10 -3 . The variance field is shown in Figure 14. The first guessed input parameter is set as 1 in the inferred domain as well.

∂ω ∂t + ∂(u j ω) ∂x j = β c γ μ t P -βω 2 + ∂ ∂x j (ν + σ ω ν t ) ∂ω ∂x j + S sst (14)

Results

Figure 15 shows that the convergence curve of the cost function with correction in the ω equation. It can be seen that the ratio of cost function J can reduce to 0.13 and for J ob it can decline to 0.25. And after approximately 30 iterations, no visible improvement can be reached.

Figure 16 presents the evolution of RANS prediction in observation space. With the correction in the ω transport equation, the stream-wise velocity u can be significantly improved even though TKE becomes inferior. With the velocity approaching the experimental data, the gradient of velocity near the wall is gradually reduced which results in the production term in the TKE transport equation related to the velocity gradient decrease accordingly. That is why it can be observed that the improvements in velocity u and TKE are mutually impeded: Once the prediction on velocity u is improved, the velocity gradient near the wall will become reduced, which results in the reduction of TKE and the further departure from the experimental measurements. The contour plots of the posterior with correction in the ω equation are presented in Figure 22(e,f). The velocity field is well reconstructed comparing to the experiments, especially in the upstream.

Correction in k and ω equation

Data assimilation setup

Because of the corrections in the TKE and special dissipation transport equation having dominant effects on the prediction of TKE and velocity u respectively, in this subsection the correction terms are introduced simultaneously in TKE and ω transport equation. The input space has 2400 dimensions including correction variables in both TKE and dissipation transport equation. Other parameters are uniform with the former subsections.

Results

Figure 18 shows that the results of the cost function with corrections in both k and ω equations. It can be seen that the value of cost function reduces significantly within 45 iterations with ratio r J of 0.087 and r J ob of 0.13. Compared with the previous cases, the efficiency of EnVar method decrease as the dimension of input parameters is increased. Figure 19 presents the evolution of the velocity u and TKE profiles by comparison with prior and experimental data. The predictions in both u and TKE is improved; however, the improvements of the predictions in the two observation are mutually restrained, which leads to that the further optimisation of the velocity would deteriorate the prediction on TKE.

From the contour plots in Figure 22(g,h), the apparent improvement on the prediction of velocity u and TKE can be seen comparing to prior. However, there is still a significant departure from the experiments. The reason why the optimised results cannot get further close to the reference data may be due to two aspects: first the experimental data has its instinct uncertainty including the measurement noise, insufficiency resolution and so on; secondly, it has been noted that the primary source of uncertainties in the RANS model is from Reynolds stress. The correction in the scalar k and ω equation is still under the frame of linear eddy-viscosity assumption; hence it can only concern the magnitude of Reynolds stress tensor but cannot take the orientation of Reynolds stress into considerations. In other words, the impact of these corrections on the observations may be not sufficient to represent the structural uncertainty in the RANS model. To this end, this work can be extended to the framework in [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed bayesian approach[END_REF] to infer the uncertainties directly in the Reynolds stress term. Moreover, the measured TKE may be difficult to be replicated by the modelled TKE since the modelled TKE cannot be equivalent to the 'true' TKE in a real flow. With the production term in k equation increasing, the resolution on the TKE increase as well. Thus, the prediction on k in the near wall region can be improved significantly comparing to the initially underestimated value. However for the region away from the wall, the production term itself becomes trivial, and thus TKE will reduce the sensitivity to the multiplicative correction of the production term; hence the assimilated results are almost similar to prior, especially for the case with corrections in both k and ω transport equation. Also, the resolution on ω gets increased with the optimal correction in the ω equation, which leads to that the modelled Reynolds stress tend to decrease, while the velocity near the wall is reduced accordingly and get close to the experimental measurements. Nevertheless, in the area far from the wall, the inferred correction is also increased but has few effects on velocity since the region is in the outer layer where the Reynolds stress is not dominant.

Conclusion

The ensemble-based variational method is presented to optimise RANS simulation by inferring improved inlet velocity and underlying model corrections in k -ω SST model. This approach combines the variational data assimilation method and ensemble Kalman method by transforming the MAP as the optimal control problem and meanwhile using the ensemble technique to estimate prior statistics so that it can sustain the advantages of both robustness and intrusiveness.

Firstly, two representative flows in the convergent-divergent channel, the turbulent flow in the WallTurb Bump where DNS resolution is available and non-cavitating flow in Venturi-type section with the X-rays experimental observation, were tested with the proposed data assimilation approach to infer the ambiguous inlet velocity profile. With improved inlet velocity, the field of velocity u can be reconstructed in good agreement with reference data mainly near the inlet region, while the velocity u in the APG region is not quite sensitive to the inlet condition. Further, the underlying model corrections in k -ω SST model were inferred for the non-cavitating flow in Venturi. The sensitivity of the correction term in the k and the ω equation is analyzed respectively. The predictions in velocity u and TKE both can be improved with corrections in the k and ω equation but still have large discrepancies comparing to the experiments which may be due to the insufficiency of the correction in the frame of Boussinesq hypothesis. The robustness of ensemble-based variational method for the inverse problem in complex turbulent flows is demonstrated.

In light of the limitation of the RANS model-form uncertainty under the frame of linear eddy viscosity assumption, current work is being conducted to explore the applicability of this data assimilation scheme to quantify and reduce uncertainties in Reynolds stress directly. Also, the EnVar method used in this work is the standard incremental ensemble version, where the background error covariance B is estimated with low-rank linear approximation and independent between consecutive data assimilation iterations. The ensemble update with consideration of the associated error covariance update will be investigated in future work [START_REF] Yang | Enhanced ensemble-based 4DVar scheme for data assimilation[END_REF]. Besides, the method utilises the limited ensemble realisations to estimate the sensitivity matrix H which may result in that the optimisation process is prone to diverge especially for the high dimension problem. Hence, the regularisation technique will be explored to be introduced in this data assimilation scheme to address this issue. Moreover, the machine learning technique will be explored to extract the model knowledge based on the inferred source term field which can be expected to be used for industrial applications [START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF].

It can be seen that the cost function cannot be further reduced after three iterations and the inferred velocity get almost stagnant, which concludes that the velocity at the region near the throat of Venturi is not sensitive to the inlet velocity especially to the velocity near the wall.

Figure 1 .

 1 Figure 1. Schematic illustration of EnVar method.

Figure 2 .

 2 Figure 2. Mesh of WallTurb Bump.

Figure 3 .

 3 Figure 3. Data assimilation results about the inferred inlet velocity and the prediction in C f : (a,b) for the prior of parabolic velocity; (c,d) for the prior of flat velocity.

Figure 4 .

 4 Figure 4. Data assimilation results of cost function J and J ob : (a,b) for the prior of parabolic inlet velocity; (c,d) for the prior of flat inlet velocity.

Figure 5 .

 5 Figure 5. Velocity U field with the first guess of parabolic velocity for prior(a), posterior(b) and DNS(c).

Figure 6 .

 6 Figure 6. Mesh of the Venturi-type section.

Figure 7 .

 7 Figure 7. Data assimilation results of velocity at inlet(a) and velocity at x = 0.0008 m(b).

Figure 8 .

 8 Figure 8. Data assimilation results of cost function J(a) and J ob (b).

Figure 9 .

 9 Figure 9. Velocity u field for prior(a), posterior(b) and experiment(c).

Figure 10 .

 10 Figure 10. Variance σ (x) of β c in TKE transport equation.

Figure 11 .

 11 Figure 11. Data assimilation results of cost function J(a) and J ob (b).

Figure 12 .

 12 Figure 12. Comparison in velocity u (a) and TKE (b) along profiles among prior, posterior and experiment.

Figure 13 .

 13 Figure 13. Inferred β c profiles.

Figure 14 .

 14 Figure 14. variance σ (x) of β c in ω transport equation.

Figure 15 .

 15 Figure 15. data assimilation results in cost function J(a) and J ob (b).

Figure 16 .

 16 Figure 16. Comparison in velocity u (a) and TKE (b) along profiles among prior, posterior and experiment.

Figure 17 .

 17 Figure 17. Inferred β c profiles.

Figure 18 .

 18 Figure 18. Data assimilation results of cost function J(a) and J ob (b).

Figure 19 .

 19 Figure 19. Comparison in velocity u (a) and k (b) along profiles among prior, posterior and experiment.

Figure 20 .

 20 Figure 20. Inferred β c profiles in k equation.

Figure 21 .

 21 Figure 21. Inferred β c profiles in ω equation.

Figure 22 .

 22 Figure 22. Contour plots of velocity u (first column) and TKE (second column): (a,b) prior; (c,d) posterior with correction in k transport equation; (e,f) posterior with correction in ω transport equation; (g,h) posterior with correction in both k and ω transport equation; (i,j) experiment.

Figures 13, 17

 17 Figures 13,[START_REF] Khlifa | Fast x-ray imaging of cavitating flows[END_REF][START_REF] Yang | Enhanced ensemble-based 4DVar scheme for data assimilation[END_REF], and 21 are the inferred model correction profiles at the corresponding position for each case. Generally, the production term in TKE and dissipation transport equation are both increased by the introduced corrections. With the production term in k equation increasing, the resolution on the TKE increase as well. Thus, the prediction on k in the near wall region can be improved significantly comparing to the initially underestimated value. However for the region away from the wall, the production term itself becomes trivial, and thus TKE will reduce the sensitivity to the multiplicative

  Figures 13,[START_REF] Khlifa | Fast x-ray imaging of cavitating flows[END_REF][START_REF] Yang | Enhanced ensemble-based 4DVar scheme for data assimilation[END_REF], and 21 are the inferred model correction profiles at the corresponding position for each case. Generally, the production term in TKE and dissipation transport equation are both increased by the introduced corrections. With the production term in k equation increasing, the resolution on the TKE increase as well. Thus, the prediction on k in the near wall region can be improved significantly comparing to the initially underestimated value. However for the region away from the wall, the production term itself becomes trivial, and thus TKE will reduce the sensitivity to the multiplicative

Figure A3 .

 A3 Figure A3. data assimilation Results to infer the inlet velocity in Venturi.

Table 1 .

 1 Data assimilation results for inlet velocity inference.

	Geometry	α	Prior	dim(α)	y	N en	r J	r J ob
	Bump		parabolic	30	C f	30	0.0025	0.20
	Venturi	inlet velocity	flat parabolic	20	u	20	0.002 2.7 × 10 -4	0.157 2 × 10 -5

Table 2 .

 2 Summary of DA experiments for model correction inference.

	α	dim(α)	y	N en	r J	r J ob
	β c in k equation	1200			0.43	0.17
	β c in ω equation	1200	TKE and u	50	0.25	0.13
	β c in k and ω equation	2400			0.087	0.13

Funding

This work is supported by the scholarship from China Scholarship Council (CSC) under the Grant CSC NO. 201607040047.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendices Appendix 1. Comparison between BFGS and Newton-CG

The performance of the minimisation method BFGS and Newton-CG is compared in the case of inlet velocity inference for Bump geometry. Figure A1 presents the evolution of cost function J to the iteration with different prior inlet velocity (parabolic or flat velocity profile). Even though both methods can result in a similar reduction in cost function after 16 iterations for the case with parabolic inlet velocity, Newton CG is faster and more robust compared to BFGS. And for the case with flat one, the BFGS method cannot reach similar results as Newton-CG within 30 iterations.

Figure A2 presents the inference of the inlet velocity and the predictions in C f accordingly. Comparing to the results with Newton CG in Figure A3, it is apparent that the eventually inferred velocity with Newton CG has a better agreement with DNS especially in the case with the flat velocity.

Appendix 2. Test case of inlet velocity inference in Venturi

We conducted a test case to infer inlet velocity in Venturi-type section. The input parameters are the inlet velocity at the first 30 grids adjacent to the bottom wall. The prior is given as the parabolic curve. The results are shown in Figure A3.