2. Knowledge-integrated additive learning of the wall model

Figure 1. Schematic of the proposed KIA approach.

The model assumes that the near-wall region of various flow states can be represented by an ODE for attached flows, augmented with a forcing term to account for non-equilibrium effects in flows with pressure gradients and separated flows. The proposed knowledge-integrated additive (KIA) learning approach learns the forcing term in a case-after-case manner (figure 1). It begins with a baseline ODE for attached flows, for which the viscous term dominates, and additively learns the forcing term to account for the neglected effects. The case-after-case continual learning without catastrophic forgetting is enabled by the encapsulation and fusion capabilities of the proposed approach:

(i) The encapsulation capability allows a modular-form model learned for each case, which consists of the predefined basis functions and neural networks for their coefficients. The basis functions define how the non-equilibrium effects are represented, while the neural networks model how different basis functions act when simulating a certain flow state.

(ii) The fusion capability synthesises models learned from various flow datasets by introducing a distance function. This distance function measures how effectively the existing model represents the new flow dataset, ensuring that only the previously unrepresented data are learned and the model’s previously acquired capabilities are preserved.

The employed TBL-s equation with a forcing term $f(\boldsymbol{u},p)$ where $\mathbf{u}$ is velocity vector, $p$ is the pressure is in the following form,

(2.1) \begin{equation} \mathcal{F}(\boldsymbol{u},p) = \frac {\partial }{\partial {y}}\left [\left (\nu + \nu _t\right ) \frac {\partial {u}}{\partial {y}}\right ] - f(\boldsymbol{u},p) =0, \end{equation}

where $\mathbf{u}$ is velocity vector, $p$ is pressure, $\nu$ is kinematic viscosity, and $u$ is the velocity parallel to the wall. The eddy viscosity $\nu _t$ in (2.1) is modelled using Prandtl’s linear mixing length model with van Driest damping function (Pope Reference Pope2000),

(2.2) \begin{equation} \nu _t=\nu \kappa {y^*}\left (1-{\textrm e}^{-y^*/A}\right )^2, \end{equation}

where Kármán constant $\kappa =0.4$ , $A=17$ , $y^*={y}u_{\tau {p}}/\nu$ , y denotes the wall-normal coordinate, and $u_{\tau {p}}=\sqrt {u_{\tau }^2 + u_p^2}$ with $u_p = |({\nu }/{\rho })({\partial {p}}/{\partial {x}})_h|^{1/3}$ the characteristic velocity defined using the pressure gradient in the tangential direction (Manhart, Peller & Brun Reference Manhart, Peller and Brun2008), with $h$ for the sampling height of the wall model (which is defined at the center of the first cell off the wall in this work).

The forcing term $f$ is approximated using the data-driven counterpart $f^J\!$ (where $J$ denotes the number of flow datasets additively employed for training), with the corresponding ODE denoted by $\mathcal{F}^J$ . In this work, the $f^J$ is in the following form:

(2.3) \begin{align} &f^J = \underbrace {\sum _{n=0}^{N}\phi _n(y)C_n^1(\mathcal{D}^{1}, \mathcal{F}^{0})}_{\text{learned from flow dataset }1} + \underbrace {\sum _{n=0}^{N}\phi _n(y)C_n^2(\mathcal{D}^{2}, \mathcal{F}^{1})}_{\text{learned from flow dataset } 2} + \underbrace {\sum _{n=0}^{N}\phi _n(y)C_n^3(\mathcal{D}^{3}, \mathcal{F}^{2})}_{\text{learned from flow dataset }3} + \dots \end{align}

(2.4) \begin{align} &= \underbrace {\phi _0(y)\!\sum _{j=0}^{J}C_0^j(\mathcal{D}^{j}\!, \mathcal{F}^{j-1})}_{\text{basis function 1}} + \underbrace {\phi _1(y)\!\sum _{j=0}^{J}C_1^j(\mathcal{D}^{j}\!, \mathcal{F}^{j-1})}_{\text{basis function 2}} + \underbrace {\phi _2(y)\!\sum _{j=0}^{J}C_2^j(\mathcal{D}^{j}\!, \mathcal{F}^{j-1})}_{\text{basis function 3}} + \dots, \end{align}

where $\phi _n$ is the predefined basis function, $\mathcal{D}^{j}$ represents the $j$ th flow dataset, $C_n^j$ is the coefficient for the basis function $\phi _n$ at $j$ th learning. In equations (2.3) and (2.4), the forcing term is written in two different forms, (i) case-wisely in the first equation, and (ii) basis-function-wisely in the second equation. The case-wise form indicates that the forcing term is expanded according to flow datasets. The basis-function-wise form signifies that the forcing term consists of a set of basis functions, with coefficients learned from various flow datasets. Different forms of basis functions can be employed. In this work, the basis function is of a polynomial form as follows:

(2.5) \begin{equation} \phi _n=\left (\frac {y}{h}\right )^n, \end{equation}

where $h$ represents the wall-normal distance for the first off-wall grid node. Although each basis function cannot be related to a specific flow regime because of the complexity of the flow, we argue that the polynomial-form basis functions provide a universal way for approximating the non-equilibrium effects. The coefficient $C_n^j$ for the basis function is computed via

(2.6) \begin{equation} C_n^j = c_n^jD_n^j, \end{equation}

where the coefficient $c_n^j$ models the added flow data, and the distance function $D_n^j$ measures how well the new data are represented using the existing model. The distance function $D_n^j$ in (2.6) equals zero if the new flow data can already be described by the model previously learned. In this work, the distance function $D^j_n$ is in the following form:

(2.7) \begin{equation} D^j_n = 1-{\textrm e}^{-{\left (B_n^j\right )}^2\left |\dfrac {\tau ^{j}_{w, r}-\tau ^{j-1}_w\left (\mathcal{D}^{j-1}, \mathcal{F}^{j-1}\right )} {A_h U_b}\right |^2}, \end{equation}

where $\tau ^{j}_{w, r}$ is the real wall shear stress from flow data set ${\mathcal{D}}^j$ , $B_n^j$ is the coefficient to be learned, $A_h=\rho /\int _0^h {\textrm {d}y}/{\nu +\nu _t}$ , and $U_b$ is the characteristic velocity for normalisation.

In these formulations, the undetermined coefficients $c_n^j$ and $B_n^j$ are modelled using a neural network with 5 hidden layers, each containing 20 neurons. In this work, the input features are $\ln ({hU_b}/{\nu } )$ , wall curvature $\kappa _w{h}$ , ${\nu {u_h}}/{hU_b^2}$ , ${u_h^2}/{U_b^2}$ , ${u_h{v_h}}/{U_b^2}$ , the tangential pressure gradient ${h}/{\rho {U_b^2}}{\partial {p}}/{\partial {x}}|_h$ , and the normal pressure gradient ${h}/{\rho {U_b^2}}{\partial {p}}/{\partial {y}}|_h$ . Here, $\kappa_w$ is the wall curvature,  $u_h$ is the velocity parallel to the wall, and $v_h$ is the velocity normal to the wall. The activation function used in the hidden layers is the hyperbolic tangent function. In the input features, a single value of $U_b$ , i.e. the bulk velocity, is employed at all spatial locations for the simulated canonical flows. For flows with complex configurations, a single value of $U_b$ may become inappropriate, and a position-specific characteristic velocity is preferred. Dividing the flow into building blocks of canonical flows is one possible way (Lozano-Durán & Bae Reference Lozano-Durán and Bae2023). The other possible way is to use the local outer flow information, which, however, requires an unambiguous definition.

With the learned forcing term, (2.1) can be solved numerically to obtain the tangential velocity in the near-wall region using the velocity boundary conditions at the wall and at $y=h$ from the outer flow. In this work, the wall-shear-stress boundary condition is employed in the outer flow simulations. It is not necessary to have the velocity distribution in the near-wall region. Taking advantage of the polynomial-form basis functions, an explicit expression for the wall shear stress $\tau _w$ can be easily obtained by integrating (2.1) twice from the wall to $h$ above the wall, and is given as follows:

(2.8) \begin{equation} \tau _w = A_h\left [U_h - \sum _{j=1}^J\left ( h{C_0^j}\int _0^h{\frac {\dfrac {y}{h}{\textrm {d}}y}{\nu + \nu _t}} + \frac {h}{2}{C_1^j}\int _0^h{\frac {\left(\dfrac {y}{h}\right)^2{\textrm {d}}y}{\nu + \nu _t}} + \dots \right ) \right ]. \end{equation}

3. A priori tests on the performance of KIA learning

The wall-resolved LES data of the periodic hill flow are employed for training and testing. The data employed for model training are from cases with Reynolds numbers $Re=2800, 5600, 10\,595, 19\,000$ and the hill slope $\alpha =1.0$ . In total, approximately 78 000 data pairs from the upper wall region and the lower hill region are utilised for model learning, which span wall-normal distances from $y/H=0.005$ to $y/H=0.1$ , where $H$ represents the height of hill, ensuring a comprehensive representation of the near-wall flow dynamics. The learning capability of the proposed approach is tested in a priori tests. The cases for model testing include (i) the periodic hill flow at a higher Reynolds number $Re=37\,000$ with the hill slope $\alpha =1.0$ , and (ii) the two cases with different hill slopes $\alpha =0.5, 1.5$ at Reynolds number $Re=10\,595$ . To test the additive learning capability, we separate the training data into two sets, i.e. the data near the upper wall and the data near the lower hills. Near the upper wall, the flow is featured by pressure gradients without flow separation, while it is characterised by a flow separation bubble near the lower hills. The baseline model is denoted as ODE_Zero. In the first additive learning, the flow data near the upper wall are employed, with the learned model dubbed as ODE_Force4_add1; in the second additive learning, the flow data near the lower hills are added to obtain the model ODE_Force4_add2. The number $4$ in the model’s name indicates that the employed basis functions are $\phi _0 = 1$ , $\phi _1 = {y}/{h}$ , $\phi _2 = ({y}/{h} )^2$ , $\phi _3 = ({y}/{h} )^3$ . More terms have been tested, showing no further improvements.

Figure 2. Performance of the models additively learned for one time (ODE_Force4_add1) and two times (ODE_Force4_add2) tested using the periodic hill cases. Time-averaged skin friction coefficients along the upper wall and lower hills for (a) $Re=10\,595$ with $\alpha =1.0$ , (b) $Re=37\,000$ with $\alpha =1.0$ , (c) $Re=10\,595$ with $\alpha =0.5$ and (d) $Re=10\,595$ with $\alpha =1.5$ . WRLES denotes wall-resolved large-eddy simulation.

The model’s performance is demonstrated in figure 2. It is seen that the baseline model, ODE_Zero, faces difficulty in predicting the wall shear stress on both the upper wall and the lower hills as expected. After the first additive learning using the upper wall data, the ODE_Force4_add1 model successfully predicts the time-averaged skin friction coefficients along the upper wall. At the lower hills, improvements were obtained in predicting the friction coefficient in the uphill region. Further improvements are achieved after the second additive learning using the data near the bottom hills. Overall, it has been shown that the model ODE_Force4_add2 additively learned twice has successfully retained its ability to predict wall shear stress at the upper wall and acquired the new capabilities to forecast $C_f$ at the surface of the lower hills.

The final note in this section concerns the impact of the learning order on model accuracy. An additional test was conducted, where the model was first trained using the lower hill data and then the upper wall data. The results indicate no further improvements after the second learning phase, as the model trained on the lower hill data already performs well for both the upper wall and the lower hill in terms of the friction coefficient.

4. Applications to wall-modelled LES

One important feature of the learned model is that it can retain the capabilities previously learned. Specifically, the model additively learned using the periodic hill data should be able to revert to the ODE without a forcing term when it is applied to wall-bounded turbulent flows with the logarithmic law satisfied. This certainly happens in the a priori tests, but becomes complicated in the a posteriori tests, for which the inputs are from the instantaneous flow fields, and the discretisation errors and the subgrid-scale-model errors also play a role. Last, the model performance needs to be tested for flow configurations other than periodic hills.

In this section, we apply the learned model to wall-modelled LES of turbulent channel flows, flows over periodic hills, and the flow over a two-dimensional (2-D) Gaussian bump, with comparisons to the reference data from wall-resolved LES or experiments and the results from the traditional wall models. In the turbulent channel flow cases with the friction Reynolds number (which is defined using the friction velocity $u_\tau$ and the half-width of the channel $\delta$ ) $Re_{\tau }=1000$ , 2000, 5200. The computational domain is $L_x \times L_y \times L_z=7.0\delta \times 2.0\delta \times 3.5 \delta$ , where $\delta$ is the half-width of the channel, in the streamwise, vertical and spanwise directions, respectively. The number of grid nodes is $N_x\times N_y \times N_z=33\times 33\times 33$ . The height of the first off-wall grid cell is $\Delta y_1 = 0.0625\delta$ . The parameters for the periodic hill cases are shown in table 1, for which systematic tests were carried out. In the Gaussian bump case, the simulation set-up follows our prior work (Zhou et al. Reference Zhou, Zhang, He and Yang2025). The geometry is defined by the equation $y(x) = h\exp [ -(x/x_0)^2 ]$ , where the bump height $h = 0.085L$ , the constant $x_0 = 0.195L$ , and $L$ represents the spanwise width of the three-dimensional bump configuration (Slotnick Reference Slotnick2019). The Reynolds number, based on $L$ and the upstream reference velocity $U_\infty$ , is $Re_L = 2\times 10^6$ . The computational domain is $L_x \times L_y \times L_z = 2.8L \times 0.5L \times 0.04L$ and is discretised using a grid with $N_x \times N_y \times N_z = 920 \times 91 \times 51$ nodes in the streamwise, vertical and spanwise directions, respectively. The first grid spacing in the vertical direction is $\Delta y_1 = 0.03h$ .

Table 1. Set-up of the periodic hill cases. Here, $Re_H$ is based on the bulk velocity and the height of the hill $(H)$ , $\Delta {y_{1}}$ is the height of the first off-wall grid node. For each case, two grid resolutions are employed: the baseline and the coarse.

The turbulent channel flow tests show good performance for all three Reynolds numbers $Re_\tau =1000, 2000, 5200$ . In figure 3, we show the turbulence statistics predicted by wall-modelled LES using the learned wall model for $Re_\tau =5200$ . As seen, the results from the learned model (i.e. ODE_Force4) show good agreement with the ODE_Zero model, the ODE_Px model (in which the forcing term $f$ in Equation (2.1) only includes the tangential pressure gradient term), and the direct numerical simulation (DNS) data for the mean streamwise velocity. About the primary Reynolds shear stress, all the wall models closely align with the DNS data for $y/{\delta }\gt 0.2$ , but underestimate the magnitude for $y/{\delta }\lt 0.2$ . All the wall models exhibit good agreement with the DNS data for the Reynolds normal stresses in the outer region. Some discrepancies are observed in the near-wall region, especially for the vertical component ( $\left \langle {v'}{v'}\right \rangle ^+$ ). Overall, the learned model aligns well with the classical models (i.e. ODE_Zero, ODE_Px), demonstrating that it can retain the capabilities previously learned. Furthermore, similar performance is observed for the models learned with different numbers of basis functions.

Figure 3. A posteriori tests using the turbulent channel flow case at Reynolds number $Re_{\tau }=5200$ for (a) mean streamwise velocity, (b) primary Reynolds shear stress $\langle u'v'\rangle^+$ and wall-normal Reynolds normal stress $\langle v'v'\rangle^+$ , (c) streamwise Reynolds normal stress $\langle u'u'\rangle^+$ and spanwise Reynolds normal stress $\langle w'w'\rangle^+$ , where $u^{\prime}$ , $v^{\prime}$ , and $w^{\prime}$ represent the velocity fluctuations in the streamwise, wall-normal, and spanwise directions, respectively. The DNS data are from Lee & Moser (Reference Lee and Moser2015). The numerical value following the term ‘Force’ in the legend indicates the number of basis functions utilised.

Figure 4. A posteriori tests using the periodic hill case at Reynolds number $Re=37\,000$ higher than the training datasets for (a) time-averaged streamwise velocity, (b) time-averaged vertical velocity, (c) primary Reynolds shear stress, and (d) turbulence kinetic energy at Reynolds number $Re=37\,000$ . The reference wall-resolved LES data are from Zhou et al. (Reference Zhou, Whitmore, Griffin and Bae2023a ).

Figure 5. A posteriori tests using the periodic hill cases with hill slopes different from the training datasets and Reynolds number $Re=10\,595$ for (a,e) time-averaged streamwise velocity, (b,f) time-averaged vertical velocity, (c,g) primary Reynolds shear stress and (d,h) turbulence kinetic energy. The reference wall-resolved LES data are from Zhou et al. (Reference Zhou, Whitmore, Griffin and Bae2023a ). The legend is the same as figure 4.

Figure 6. A posteriori tests using the periodic hill cases. Time-averaged skin friction coefficients along the upper wall and lower hills for (a) $Re=10\,595$ with $\alpha =1.0$ , (b) $Re=37\,000$ with $\alpha =1.0$ , (c) $Re=10\,595$ with $\alpha =0.5$ and (d) $Re=10\,595$ with $\alpha =1.5$ .

The capability of the learned model for wall-modelled LES of separated flows is demonstrated via the periodic hill flow cases for Reynolds numbers higher than the training data and different hill slopes. The performance is first examined for the $Re=37\,000$ case (case PH37 000- $\alpha$ 1.0 in table 1), for which the Reynolds number is higher than the training flow data. Although some minor discrepancies exist in the separation bubble region for the turbulence kinetic energy as shown in figure 4, the learned models show an overall good agreement with the wall-resolved LES data for the time-averaged streamwise velocity ( $\langle u \rangle$ ), vertical velocity ( $\langle v \rangle$ ), primary Reynolds shear stress ( $\langle u^{\prime } v^{\prime }\rangle$ ) and turbulence kinetic energy ( $k$ ) for both grid resolutions. Figure 5 shows the time-averaged results of flow over periodic hills with hill slopes $\alpha =0.5$ and $\alpha =1.5$ (cases PH10 595- $\alpha$ 0.5 and PH10 595- $\alpha$ 1.5 in table 1), being different from the training flow data. As seen, in comparison with the baseline model, the ODE_Force4 model follows well the reference data for the time-averaged streamwise velocity and vertical velocity as well as the Reynolds shear stress and turbulence kinetic energy. The underprediction of the magnitude of $k$ at certain streamwise locations is due primarily to the grid resolution employed in wall-modelled LES. The level of agreement improves when compared with the $k$ values from filtered wall-resolved LES results, as demonstrated in our prior work (Zhou et al. Reference Zhou, Zhang, He and Yang2025). The time-averaged skin friction coefficients from the wall-modelled LES are compared with the reference data in figure 6. While some discrepancies are observed at certain streamwise locations, the ODE_Force4 model achieves an overall good agreement.

Figure 7. A posteriori tests using the periodic hill cases. Relative errors computed using (4.1) for (a) time-averaged streamwise velocity, (b) time-averaged vertical velocity, (c) primary Reynolds shear stress and (d) turbulence kinetic energy.

To quantitatively assess the prediction accuracy of the proposed wall models, we compute the relative error at each streamwise location as follows:

(4.1) \begin{equation} err_q = \frac {\sum |q_{{ref}} - q_{{wm}}|\Delta {y}}{\sum |q_{{ref}}|\Delta {y}}, \end{equation}

where $q$ denotes a flow quantity with the subscripts ‘ref’ and ‘wm’ for the reference and wall-modelled LES data, respectively. The prediction accuracy of the models is evaluated using the periodic hill cases, for which results are available across various models and grid resolutions. Specifically, the errors of the vertical profiles, shown in figures 4 and 5, are first computed at each streamwise location and then averaged for each case. The obtained errors are shown in figure 7. As observed, the errors are significantly reduced when comparing the results of the ODE_Force4 model with those of the ODE_Zero model for both grid resolutions. For the baseline grid, the average relative errors for the ODE_Force4 model are below 5 % for the mean streamwise velocity and under 20 % for the other three quantities, respectively.

The results from the 2-D Gaussian bump case are presented in figures 8 and 9. The four streamwise locations selected for comparison, i.e. $x/L=-0.4,-0.1,0.1,0.2$ , correspond to regions with a mild adverse pressure gradient (APG), a strong favourable pressure gradient (FPG), a strong APG, and flow separation, respectively. As seen, the three models perform reasonably well at locations with APG and FPG. However, at the location downstream of the bump ( $x/L=0.2$ ), the ODE_Zero model fails to accurately predict the flow separation. For the ODE_Force4 model, although some discrepancies are observed around $y/L=0.07$ , it accurately predicts the vertical distribution of $\langle u \rangle$ in the flow separation region.

Interesting trends are observed in the comparison between the ODE_Force4 and $\text{FEL}_{\tau _w}$ predictions. Specifically, the ODE_Force4 model performs better in the upstream of the bump but is less accurate in the separation region, while the opposite is true for the $\text{FEL}_{\tau _w}$ model. This indicates that the history effect from the upstream flow (Bobke et al. Reference Bobke, Vinuesa, Örlü and Schlatter2017; Parthasarathy & Saxton-Fox Reference Parthasarathy and Saxton-Fox2023; Yang et al. Reference Yang, Chen, Zhang and Kunz2024) plays a non-negligible role in the streamwise flow evolution. However, this effect is not explicitly included in either model and should be accounted for in future model development.

Lastly, the computational overhead is assessed for the ODE_Force4 model in comparison with the ODE_Zero model. The increase in the wall-clock time is found to be modest, with the ratio ranging from 1.02 for the channel case to 1.16 for the 2-D Gaussian bump case.