1. Introduction
Reynolds-averaged Navier–Stokes (RANS), supported by turbulence closure models, remains a predominant computational fluid dynamics (CFD) approach in industry. The RANS approach maintains a balance between computational efficiency and accuracy for predicting mean flow quantities, making it critical for applications across aerospace, automotive and energy sectors (Rumsey Reference Rumsey2014; Bush et al. Reference Bush, Chyczewski, Duraisamy, Eisfeld, Rumsey and Smith2019). The main difficulty, however, lies in closing the Reynolds stress tensor, especially for complex flows such as boundary layer transition, separation and secondary flows. Traditional closures are based on physical arguments, empirical tuning or, more recently, data-driven calibration (Launder & Sharma Reference Launder and Sharma1974; Spalart & Allmaras Reference Spalart and Allmaras1992; Wilcox Reference Wilcox2008; Ling, Kurzawski & Templeton Reference Ling, Kurzawski and Templeton2016).
The demand for improving turbulence models has led to the use of machine learning methods (Duraisamy, Iaccarino & Xiao Reference Duraisamy, Iaccarino and Xiao2019; Brunton, Noack & Koumoutsakos Reference Brunton, Noack and Koumoutsakos2020). Many machine learning methods have been introduced into turbulence modelling, including Bayesian optimisation (Xiao et al. Reference Xiao, Wu, Wang, Sun and Roy2016), ensemble Kalman filter (Zhang et al. Reference Zhang, Xiao, Luo and He2022), adaboost decision trees (Ling & Templeton Reference Ling and Templeton2015), random forests (Wang, Wu & Xiao Reference Wang, Wu and Xiao2017), tensor basis neural networks (Ling et al. Reference Ling, Kurzawski and Templeton2016), and symbolic regression (Weatheritt & Sandberg Reference Weatheritt and Sandberg2016; Schmelzer, Dwight & Cinnella Reference Schmelzer, Dwight and Cinnella2020; Zhao et al. Reference Zhao, Akolekar, Weatheritt, Michelassi and Sandberg2020), among others. In these studies, the modelling strategies vary considerably, from correcting the parameters of linear eddy viscosity models to forming completely new nonlinear eddy viscosity models, with the optimisation targets varying from predicting transition (Duraisamy, Zhang & Singh Reference Duraisamy, Zhang and Singh2015; Zhang et al. Reference Zhang, Chen, Zeng, Tang and Yan2023) to characterising boundary layer separation (Zhu et al. Reference Zhu, Zhang, Kou and Liu2019; Wu & Zhang Reference Wu and Zhang2023). These models, which are typically trained to improve predictive accuracy for specific cases, can thus be termed as expert models.
Despite considerable advances, achieving robust generalisation for data-driven turbulence models remains not only an active research area but also a significant challenge. Based on the accumulated experience from industrial turbulence simulations, evidence suggests that a universal, simple and local turbulence model may be difficult to achieve (Rumsey & Coleman Reference Rumsey and Coleman2022). To enhance generalisation for traditional turbulence models, one approach is to perform parameter calibration in various application scenarios. The GEKO model, for example, provides empirical parameter combinations for the
$k{-}\omega$
model under a wide range of cases (Menter & Matyushenko Reference Menter and Matyushenko2025). This cross-case generalisation capability, i.e. the capability to make accurate predictions in diverse flow scenarios, is particularly important and even more challenging for models trained by machine learning methods (Sandberg & Zhao Reference Sandberg and Zhao2022).
To achieve generalisation across diverse flow scenarios, researchers have explored numerous approaches that leverage model architecture design, training data diversity, feature engineering and physical constraints. These efforts can be broadly categorised into three distinct paradigms. The first strategy focuses on parameter calibration using data-driven approaches to enhance adaptability. For instance, Bin, Huang & Yang (Reference Bin, Huang and Yang2023) constructed a unified model by fine-tuning the parameters of the Spalart–Allmaras (SA) model, thereby extending its applicability across various flows. The second approach seeks to improve generalisation through multi-case training strategies, where the model is exposed to distinct geometries and physical characteristics during the learning phase. Notable examples include the symbolic-regression-based training developed by Fang et al. (Reference Fang, Zhao, Waschkowski, Ooi and Sandberg2023), the multi-case surrogate optimisation by Amarloo et al. (Reference Amarloo, Rincón, Reclari and Abkar2023), and the progressive data-augmented framework to incrementally enhance cross-case generalisation by Rincón et al. (Reference Rincón, Reclari, Yang and Abkar2025). The third category involves model aggregation or ensemble methods, which utilise machine-learning-based weighting functions to blend multiple baseline models. In a series of recent studies (de Zordo-Banliat et al. Reference de Zordo-Banliat, Dergham, Merle and Cinnella2024; Cherroud et al. Reference Cherroud, Merle, Cinnella and Gloerfelt2025; Oulghelou et al. Reference Oulghelou, Cherroud, Merle and Cinnella2025), machine learning methods were employed to train weighting functions, successfully leveraging the strengths of different data-driven models to achieve improved accuracy across multiple cases. Parallel to these RANS-based developments, significant progress has also been made in wall-modelled large-eddy simulations (WMLES) (Lozano-Durán & Bae Reference Lozano-Durán and Bae2023; Arranz et al. Reference Arranz, Ling, Costa, Goc and Lozano-Durán2024), in which a model trained on several canonical flows exhibits good agreement with reference data for both canonical tests and realistic aircraft configurations.
It should be noted that most existing data-driven models,enhanced by either multi-case joint training or model aggregation, are considered reliable mostly within the vicinity of their training distributions, while applying to completely unknown flow regimes remains a critical challenge. In practical CFD applications, users are frequently confronted with novel flow configurations and operating conditions distinct from prior experiences. Directly training the existing model with newly obtained data often degrades its original performance, which is known as catastrophic forgetting (French Reference French1999). A common approach is to retrain the model jointly on both pre-existing and new data, but this approach might face a series of challenges. As new regimes arise continuously in realistic applications, joint retraining can be hindered by complex optimisation tasks, data accessibility issues, and performance degradation on previously well-modelled cases, all of which motivate continual learning. In order to achieve continual learning for WMLES, Zhang et al. (Reference Zhang, Zhou, Yang and He2025) proposed an additive framework for training wall models, showing promising success in a series of cases. Nevertheless, novel frameworks with the capability of continual learning are still needed for developing generalisable RANS models.
To achieve a framework capable of sustainable generalisation, we propose a novel approach termed progressive mixture-of-experts (PMoE) in the present study. The PMoE method is founded on the mixture-of-experts (MoE) architecture (Jacobs et al. Reference Jacobs, Jordan, Nowlan and Hinton1991), which has been the subject of much research in the field of machine learning (Shazeer et al. Reference Shazeer, Mirhoseini, Maziarz, Davis, Le, Hinton and Dean2017; Fedus, Zoph & Shazeer Reference Fedus, Zoph and Shazeer2022) and has become very popular in the field of large language models due to its scalability (Du et al. Reference Du2022; Dai et al. Reference Dai2024). Building on the typical architecture consisting of a router and a group of experts, the MoE offers a distinct advantage for turbulence modelling, i.e. the decomposition of a complex, high-dimensional physical problem into manageable sub-tasks. Central to this architecture is the mechanism of sparse activation. Unlike monolithic models, the router in an MoE network activates only a specific subset of expert models for a given local flow state, leaving the majority silent during the inference process. This property is crucial for industrial CFD applications, as it does not incur the prohibitive computational cost usually associated with large-scale neural networks, while achieving sufficient representational expressiveness in order to deal with various scenarios featuring different flow physics.
One important novel feature of the present PMoE framework, which is missing for the generic MoE structure, is the continual learning capability. This is in order to achieve sustainable generalisation, enabling the model to incorporate newly introduced datasets without full retraining. Accordingly, the PMoE is designed to be equipped with a modular autoencoder-based router. The PMoE framework is applied to develop a turbulence model by successively introducing the data of four flow regimes, including aerofoil wake, channel, periodic hill and square duct flows. The trained PMoE model is then extensively assessed, validating its performance in training cases and unseen cases with various operating conditions and configurations. The objective is mainly to show the ability of PMoE to incorporate new flow scenarios without retraining or degrading existing experts.
Before introducing the framework in detail, it is important to highlight three key issues that need to be addressed to develop a PMoE turbulence model generalisable for various flow regimes. First, the PMoE framework needs to classify diverse flow scenarios in an unsupervised manner using only RANS data, and this requirement brings a significant challenge regarding the selection and extraction of appropriate physical features for the pre-processing procedure to serve as effective model inputs. Second, the router architecture must support continual learning capabilities, determining whether a new case belongs to a known flow regime or represents a novel one. For the latter, the architecture must allow for the rapid addition of new modules without the need for retraining. Finally, the distinct nature of turbulence across different flows implies that a unified model form may be insufficient (Duraisamy et al. Reference Duraisamy, Iaccarino and Xiao2019), and different flow regimes often require fundamentally different correction forms. Consequently, the design of specific expert architectures tailored to these varying physical characteristics is crucial for the success of a truly generalisable model. These topics will be addressed in detail in the following sections.
This paper is organised as follows. The structure of the new framework is discussed in § 2. The training process and the setting of cases are introduced in § 3. The model assessment, validation and computational cost analysis are given in § 4. Finally, conclusions are offered in § 6.
2. Methodology
The MoE framework, which has the potential to develop generalisable turbulence models, is introduced in § 2.1. Moreover, the router part of a generic MoE framework is typically implemented by multi-layer perceptron (MLP), which enables multi-task learning but requires further improvement for continual learning. In § 2.2, we introduce the PMoE structure, demonstrating how the newly proposed framework facilitates continual generalisation of turbulence modelling.
2.1. Introduction to the generic MoE framework
As shown in figure 1, the generic MoE consists of a set of
$n$
experts
$\{E_1, \ldots ,E_n\}$
and a gating network (also named a router)
$G$
, and the key advantage of this framework is the decomposition of a complex, high-dimensional problem into manageable sub-tasks dealt by different experts. For turbulence modelling, each expert
$E_i$
can be viewed as a specialised closure model tailored for a specific flow regime, while the router
$G$
functions as a flow regime classifier.
Schematic of the generic MoE framework with an MLP Softmax gating network. (a) The whole process of the MoE framework. (b) Zoom-in view of a gating network based on MLP.
Figure 1. Long description
The image illustrates the generic Mixture of Experts (MoE) framework featuring a Multi-Layer Perceptron (MLP) Softmax gating network. In part (a), the overall process of the MoE framework is depicted. It shows an expert group consisting of multiple experts labeled E1, E2, ..., En, each contributing to the output. A gating network, represented by a computer and server, processes the input x and generates gating values G(x)i for each expert. These values determine the contribution of each expert to the final output. Part (b) provides a detailed view of the gating network based on MLP, showing a neural network structure with interconnected nodes. The gating network-MLP uses a Softmax function to produce the gating values.
For a set of
$n$
experts, the output
$\boldsymbol{y}$
for a given input flow state
$\boldsymbol{x}$
is the weighted sum of expert predictions:
\begin{equation} \boldsymbol{y} = \sum _{i=1}^{n} G(\boldsymbol{x})_i\, E_i(\boldsymbol{x}). \end{equation}
Here,
$G(\boldsymbol{x})_i$
is the
$i$
th component of the output vector of
$G$
, and represents the probability (or confidence) that the corresponding flow
$\boldsymbol{x}$
belongs to the regime governed by expert
$i$
. A simple choice of the gating function (Jordan & Jacobs Reference Jordan and Jacobs1994), as shown in figure 1(b) for example, is to multiply the input
$\boldsymbol{x}$
by a trainable weighting matrix
$\unicode{x1D652}_g$
and then apply the
$\textrm{Softmax}({\cdot })$
function, as
\begin{equation} G_{\sigma }(\boldsymbol{x})=\textrm{Softmax}\left (\boldsymbol{x} \unicode{x1D652}_{g}\right ), \end{equation}
where the subscript
$\sigma$
denotes the Softmax nonlinearity.
The key advantage of the MoE framework lies in its sparse activation mechanism, which dramatically reduces computational costs while maintaining model capacity. Typically, sparse MoE employs a top-
$k$
routing strategy where only the most relevant
$k$
experts are activated for each input. For example, we set
$k=1$
in our current implementation, thus the router selects the one most relevant expert per input token, ensuring that the other
$n-1$
experts remain completely dormant during model inference. With this feature of sparse activation, the computational cost of the MoE model, despite having a large pool of experts
$\{E_1, \ldots ,E_n\}$
, scales only with the size of a single expert rather than the entire ensemble. Consequently, the inference cost remains nearly constant regardless of the total number of experts, enabling the deployment of extremely large-scale models with billions or even trillions of parameters with relatively low computational requirements (Shazeer et al. Reference Shazeer, Mirhoseini, Maziarz, Davis, Le, Hinton and Dean2017; Fedus et al. Reference Fedus, Zoph and Shazeer2022).
While the generic MoE offers a path towards modular turbulence modelling, applying it to industrial CFD requires a capability that standard architectures lack, i.e. continual learning. A truly generalisable turbulence model must be able to assimilate new flow physics (e.g. transitioning from simple shear flows to complex corner separations) without ‘forgetting’ previously learned regimes. However, standard MLP-based routers in figure 1 are highly susceptible to catastrophic forgetting (French Reference French1999), as updating the weights to accommodate a new flow regime often degrades the classification accuracy for existing ones, highlighting the urgent need to extend the framework with continual learning capability.
In the broader machine learning literature, continual learning is attracting more and more interest, and the corresponding strategies can be generally classified into three paradigms (Wang et al. Reference Wang, Zhang, Su and Zhu2024): (i) regularisation-based methods, which constrain weight updates to preserve important parameters (e.g. Kirkpatrick et al. Reference Kirkpatrick2017); (ii) replay-based methods, which retain or generate representative samples from earlier tasks for joint retraining (e.g. Rolnick et al. Reference Rolnick, Ahuja, Schwarz, Lillicrap and Wayne2019); and (iii) architecture-based methods, which dedicate separate model components to different tasks to prevent interference (e.g. Rusu et al. Reference Rusu, Rabinowitz, Desjardins, Soyer, Kirkpatrick, Kavukcuoglu, Pascanu and Hadsell2016). Among architecture-based approaches, Aljundi, Chakravarty & Tuytelaars (Reference Aljundi, Chakravarty and Tuytelaars2017) proposed Expert Gate, a lifelong learning framework in which each task is associated with a dedicated autoencoder gate and expert network. At test time, the autoencoder with the lowest reconstruction error routes the input to its corresponding expert, while the reconstruction errors also quantify inter-task relatedness to guide knowledge transfer when training new experts. This design demonstrated effective sequential task learning for image classification. The present PMoE framework builds on this broad principle of modular autoencoder-based routing, while introducing several novel strategies essential for turbulence modelling, as discussed in the following sections.
2.2. The PMoE framework for turbulence modelling
As shown by the schematic in figure 2, the PMoE framework consists of three key components, including the procedure of pre-processing to extract flow features as model inputs, the router structure responsible for classifying various types of flows, and a group of experts each trained for a specified flow scenario. In particular, the PMoE replaces the MLP router with a modular bank of autoencoders, each dedicated to recognising a specific flow signature. This ensures that the introduction of a new expert and its corresponding router module does not interfere with the existing modules. Accordingly, the PMoE framework is expected to develop progressively with new datasets introduced successively, and the details are discussed in the following.
Schematic of the proposed PMoE framework. The framework consists of three parts: (i) a pre-processing stage that extracts physics-informed features from the baseline RANS field; (ii) a modular autoencoder-based router that classifies the incoming flow regime by comparing reconstruction errors across components, and triggers model expansion when no existing component achieves sufficient confidence; and (iii) a sparsely activated group of heterogeneous experts, each tailored to a specific flow regime. Here, SR means symbolic regression, and FI-NN means field inversion neural network.
Figure 2. Long description
The diagram illustrates the proposed PMoE framework, which includes three main parts: pre-processing, a router, and a group of experts. The pre-processing stage extracts physics-informed features from the baseline RANS field. The router, which is a modular autoencoder-based system, classifies the incoming flow regime by comparing reconstruction errors across components and triggers model expansion when no existing component achieves sufficient confidence. The experts are a sparsely activated group of heterogeneous models, each tailored to a specific flow regime. The diagram shows the flow of information from the pre-processing stage to the router and then to the experts.
2.2.1. Feature extraction and physics-informed sampling
During the pre-processing stage, the input to the PMoE framework relies on local flow features extracted from a baseline RANS calculation. In the present study, we use the SA model (Spalart & Allmaras Reference Spalart and Allmaras1992) as the baseline due to its wide applications. Following previous data-driven RANS efforts (Ling & Templeton Reference Ling and Templeton2015; Wang et al. Reference Wang, Wu and Xiao2017), we construct an input vector
$\boldsymbol{x} = [q_1, q_2, \ldots , q_M]$
consisting of
$M=7$
flow features, as detailed in table 1. Note that the listed variables are normalised using
$A/({\left |A \right |+\left |B \right |})$
to ensure consistent scaling across different flow regimes.
Flow features used as inputs in the framework:
$\unicode{x1D652}$
is the rotation-rate tensor,
$\unicode{x1D64E}$
is the strain-rate tensor,
$\rho$
is the fluid density,
$\nu$
is the fluid viscosity,
$\tilde {\nu }$
is the SA modified turbulent eddy viscosity,
$d$
is the distance to the nearest wall,
$L_{\textit{ref}}$
is the reference height, and
$U_i$
is the mean velocity. The normalised feature
$q_i$
is obtained by normalising the corresponding raw features value
$\hat {q}_i$
with normalisation factor
$q_i^*$
according to
$q_i={\hat {q}_i}/({ |\hat {q}_i |+ |q_i^* |})$
.
Table 1. Long description
The table presents flow features used as inputs in a framework, listing seven key features with their descriptions, raw features, and normalization factors. The features include vorticity magnitude, ratio of excess rotation rate to strain rate, ratio of pressure normal stresses to shear stresses, distance to the nearest wall, SA viscosity, pressure gradient along streamline, and non-orthogonality between velocity and its gradient. Each feature is detailed with its corresponding mathematical representation and normalization factor to ensure consistent scaling across different flow regimes.
We note that the features listed in table 1 are not universally effective across all flow regimes. Prior studies (Yin et al. Reference Yin, Yang, Zhang, Chen and Fu2020; de Zordo-Banliat et al. Reference de Zordo-Banliat, Dergham, Merle and Cinnella2024; Cherroud et al. Reference Cherroud, Merle, Cinnella and Gloerfelt2025) suggest that a subset of these variables may often be sufficient for specific cases, and methods such as permutation feature importance (PFI) (Wu et al. Reference Wu, Zhang, Xu and He2025) can be employed to identify the relevant features for a specific case. Nevertheless, to accommodate future applications that may involve a broad range of flow regimes with distinct critical features, we adopt the full set as the model input.
A critical challenge in training a generalisable router is ensuring that the training data represent the dynamically active regions of the flow, rather than being dominated by the free stream or quiescent zones. Therefore, we employ a physics-informed sampling strategy during the pre-processing stage. Rather than uniformly sampling the entire computational domain, the sampling region is restricted based on the underlying flow physics. Grid points are selected only where either the strain-rate magnitude
$\|{\unicode{x1D64E}}\|$
or the vorticity magnitude
$\|\unicode{x1D652}\|$
exceeds a specific threshold relative to the domain statistics. Specifically, we define the sampling mask
$\mathcal{M}$
such that a point is included if
\begin{equation} \|\unicode{x1D652}_j\| \gt \epsilon \times\textrm{median}( \|\unicode{x1D652}\|) \quad \textrm{or} \quad \displaystyle \|{\unicode{x1D64E}}_j\| \gt \epsilon \times\textrm{median}( \|{\unicode{x1D64E}}\|), \end{equation}
where
$\|\unicode{x1D652}_j\|$
and
$\|{\unicode{x1D64E}}_j\|$
denote the magnitudes of vorticity and strain-rate at the
$j$
th sampled point, respectively, and
$\epsilon$
is a thresholding parameter set to
$0.05$
in the present study. A sensitivity test has been performed, and the results are robust within
$\epsilon \in \{0.01, 0.05, 0.1\}$
. This criterion ensures that the selected samples are concentrated in boundary layers, shear layers and separation zones, etc., effectively filtering out grid points of limited physical significance such as those in free-stream regions. From this masked region,
$N$
points are randomly sampled to form the input matrix
${\unicode{x1D653}} \in \mathbb{R}^{N \times M}$
. The masking criterion in (2.3) is applied only during the sampling stage for the PMoE router. In regions with limited physical significance, such as free-stream regions, the baseline RANS predictions are generally adequate to provide accurate predictions, thus no extra corrections are needed.
2.2.2. Modular autoencoder-based router for continual learning
The core of the PMoE framework is its ability to determine whether the incoming case belongs to a known scenario or represents a novel one that requires model expansion. As illustrated in figure 3, the PMoE router features a modular architecture, in which each router component adopts a standard autoencoder architecture following the design of Hinton & Salakhutdinov (Reference Hinton and Salakhutdinov2006), aiming to compress and reconstruct the feature space of a specific flow regime. The idea of learning internal representations through minimising input reconstruction traces back to the back-propagation framework of Rumelhart, Hinton & Williams (Reference Rumelhart, Hinton and Williams1986). The network in this study consists of an encoder
$f_{\theta }$
mapping the input
$\boldsymbol{x}$
to a latent representation
$\boldsymbol{z}$
, and a decoder
$g_{\phi }$
reconstructing it as
$\hat {\boldsymbol{x}}$
. The network parameters of the encoder and the decoder shown in figure 3(a) are represented by
$\theta$
and
$\phi$
, respectively. The training objective for the autoencoder is to minimise the reconstruction loss over the
$N$
sampled points, as
\begin{equation} \boldsymbol{z}=f_{\theta }(\boldsymbol{x}), \quad \boldsymbol{\hat {x}}=g_{\phi }(\boldsymbol{z}), \quad \theta , \phi =\arg \min _{\theta , \phi } \mathcal{L}(\boldsymbol{x}, \boldsymbol{\hat {x}}). \end{equation}
Here, the reconstruction error
$\mathcal{L}(\boldsymbol{x},\boldsymbol{\hat {x}})$
is assessed as
\begin{equation} \mathcal{L}(\boldsymbol{x}, \boldsymbol{\hat {x}})=\frac {1}{N}\sum _{j}^{N}L_j=\frac {1}{N}\sum _{j}^{N}\sqrt {\frac {1}{M}\sum _{i}^{M} (\hat {x}_{i,j}-x_{i,j}) ^2 }, \end{equation}
where for the
$j$
th point,
$x_{i,j}$
is the
$i$
th component of its input feature vector, and
$L_j = \sqrt {({1}/{M}) \sum _{i=1}^{M} ( \hat {x}_{i,j} - x_{i,j} )^2 }$
refers to its corresponding reconstruction error over the
$M$
features.
Structure of the router based on an autoencoder. (a) Neural network architecture design for the autoencoder component. (b) The modular architecture of the router.
Figure 3. Long description
The image consists of two parts labeled (a) and (b). Part (a) depicts the neural network architecture design for the autoencoder component. It shows layers labeled as Input layer, Encoder, Bottleneck layer, Decoder, and Output layer, with varying numbers of nodes in each layer. Part (b) illustrates the modular architecture of the router, featuring multiple encoders and decoders processing input features and calculating loss. The loss is represented as the difference between the predicted output and the actual input.
For every router component, a compact architecture is adopted, consisting of several hidden layers with moderate width to balance representational capacity and generalisation as shown in figure 3(a). Leaky rectified linear unit (Maas, Hannun & Ng Reference Maas, Hannun and Ng2013) activations are employed between layers to introduce nonlinearity, and the Adam optimiser (Kingma & Ba Reference Kingma and Ba2015) is used with initial learning rate
$10^{-2}$
. A step-based learning-rate scheduler is employed, where the learning rate is reduced by a factor
$0.98$
every
$100$
epochs to improve convergence stability. The number of training epochs is capped at
$5000$
based on empirical convergence studies, while an early-stopping strategy is employed to improve efficiency and mitigate overfitting. For all component trainings, the loss consistently reaches a plateau well before the maximum epoch count, and extending the training further does not yield noticeable performance improvement. After the training of the
$k$
th router component
$C_k$
is converged, we set the
$99.9$
th percentile of the reconstruction error of the training data as the corresponding threshold
$T_k$
, in order to eliminate the misleading influence of points with extreme distributions. A sensitivity test has been performed, and the results are robust with
$T_k\gt 99\,\%$
.
During the inference process, the router must determine if a new input flow belongs to a known regime or represents a novel physical scenario through the reconstruction error. For a new input case with
$N$
sampled points, a point
$j$
is considered recognised by component
$C_k$
if its local reconstruction error satisfies
$L_{j} \lt T_k$
. The global confidence level
$p_k$
of expert
$k$
regarding the current flow is defined as the relative frequency of recognised points
\begin{equation} p_k = \frac {N_k}{N}, \quad \text{where } N_k = \sum _{j}^{N} \mathbb{I}(L_j \lt T_k), \end{equation}
where
$\mathbb{I}({\cdot })$
is the indicator function.
The framework employs a winner-takes-all or top-
$1$
gating strategy. The router selects the expert
$E_{\mathcal{K}}$
corresponding to the component with the maximum confidence as
\begin{equation} \mathcal{K} = \arg \max _{k} p_k. \end{equation}
However, activation is conditional and governed by a global acceptance threshold
$T_{\textit{accept}}$
. When the maximum confidence
$p_{\mathcal{K}}$
exceeds a threshold
$T_{\textit{accept}}=90\,\%$
, the expert
$E_{\mathcal{K}}$
will be activated. Conversely, if
$p_{\mathcal{K}}$
falls below
$T_{\textit{accept}}$
, then the flow is deemed unknown, which triggers the continual learning process. In this case, a new autoencoder component
$C_{\textit{ne}w}$
and a new expert
$E_{\textit{ne}w}$
are initialised and trained on the new dataset. It should be emphasised that the expansion process is strictly additive. When a new regime is introduced, the corresponding router component and expert are trained, while all previously learned components and experts remain frozen. This modular design ensures isolation between old and new components during training, thus inherently overcoming catastrophic forgetting. Note that while a threshold
$T_{\textit{accept}}=90\,\%$
is found to be effective for the distinct regimes in the present study, future industrial applications may require this hyperparameter to be tuned based on the desired sensitivity to novel flow physics.
The rationale for employing an autoencoder for flow classification rests on the premise that the hierarchical importance of physical features varies significantly across different flow scenarios, which can be seen from the feature importance analysis in Appendix A. During the routing stage, an input feature matrix
${\unicode{x1D653}}\in \mathbb{R}^{N\times M}$
is compressed into a low-dimensional bottleneck representation
$\mathcal{H}_\theta = f_\theta ({\unicode{x1D653}})$
to extract the feature combinations that are most critical for reconstructing the local flow physics. Following the information bottleneck principle proposed by Tishby, Pereira & Bialek (Reference Tishby, Pereira and Bialek2000), the training objective is formulated as
\begin{equation} \min _{\theta }\left [I({\unicode{x1D653}};\mathcal{H}_\theta )-\gamma I(\mathcal{H}_\theta ;{\unicode{x1D654}})\right ], \end{equation}
where
$\theta$
denotes the encoder parameters,
$I(\boldsymbol{\cdot };\boldsymbol{\cdot })$
represents mutual information, and
$\gamma$
governs the trade-off between compression and relevance. For the reconstruction task, the output matrix
${\unicode{x1D654}}\in \mathbb{R}^{N\times M}$
is taken to be consistent with the input
$\unicode{x1D653}$
. The information bottleneck formulation implies the existence of a minimal informative feature subset
$S \subseteq \{1, 2, \ldots , M\}$
such that the sub-matrix
${\unicode{x1D653}}_S\in \mathbb{R}^{N\times |S|}$
satisfies
\begin{equation} I({\unicode{x1D653}}_S;{\unicode{x1D654}})=I({\unicode{x1D653}};{\unicode{x1D654}}),\quad I({\unicode{x1D653}}_{S\prime};{\unicode{x1D654}})\lt I({\unicode{x1D653}};{\unicode{x1D654}}),\quad \forall S^\prime\subset S . \end{equation}
This implies that the autoencoder effectively isolates feature combinations that are strictly relevant to flow reconstruction. Consequently, for distinct flow scenarios, the network relies on different minimal subsets. This mechanism enables flow identification by exploiting the dependence of the latent representation
$\mathcal{H}_\theta$
on flow-specific informative features. Note that the dimension of the latent bottleneck layer is a critical parameter that determines how much the input data can be compressed. We set the latent layer to have four nodes in the present study as shown in figure 3(a), which is based on a series of sensitivity tests.
2.2.3. Sparsely activated expert group with heterogeneous formulation
A distinct advantage of the PMoE framework is its explicit exploitation of heterogeneous expert formulations. Unlike data-driven models that typically enforce a single-form closure correction, e.g. a neural network predicting Reynolds stress anisotropy (Ling et al. Reference Ling, Kurzawski and Templeton2016), the PMoE allows each expert
$E_k$
to adopt the mathematical structure best suited for its specific flow regime. This flexibility ensures that simple flows can be modelled by interpretable algebraic corrections, while complex non-equilibrium flows can leverage high-capacity neural networks. We remark that although conventional MoE architectures can in principle accommodate different expert structures, standard implementations almost universally employ homogeneous experts due to the requirement of joint end-to-end training. In the PMoE framework, each expert is trained independently, removing this constraint and enabling fundamentally different closure strategies.
In the present implementation, we utilise the SA model (Spalart & Allmaras Reference Spalart and Allmaras1992) as the baseline model. The standard transport equation for the modified eddy viscosity
$\tilde {\nu }$
is given by
\begin{equation} \frac {\textrm{D} \tilde {\nu }}{\textrm{D} t} = \underbrace {C_{b1} \tilde {S} \tilde {\nu }}_{\mathcal{P}} + \underbrace {\frac {1}{\sigma }\big [\boldsymbol{\nabla }\boldsymbol{\cdot }((\nu +\tilde {\nu })\, \boldsymbol{\nabla }\tilde {\nu })+C_{b2}( \boldsymbol{\nabla }\tilde {\nu })^{2}\big ]}_{\mathcal{T}} - \underbrace {C_{w1} f_{w}\left (\frac {\tilde {\nu }}{d}\right )^{2}}_{\mathcal{D}}, \end{equation}
where
$\mathcal{P}$
,
$\mathcal{T}$
and
$\mathcal{D}$
represent the production, transport and destruction terms, respectively. The standard definitions for the model constants (such as
$C_{b1}$
,
$C_{b2}$
,
$\sigma$
and
$C_{w1}$
) and closure functions such as
$f_w$
are retained from Spalart & Allmaras (Reference Spalart and Allmaras1992), and we only present the terms that will be modified here for brevity. In particular, the turbulent eddy viscosity is computed from
$\nu _{t}=\tilde {\nu } f_{\nu 1}$
, and
$f_{\nu 1}$
is defined as
\begin{equation} f_{\nu 1}=\frac {\chi ^{3}}{\chi ^{3}+c_{\nu 1}^{3}}, \quad \chi \equiv {\tilde {\nu }}/{\nu }, \end{equation}
where
$c_{\nu 1}=7.1$
, and
$\chi$
represents the ratio between
$\tilde {\nu }$
and the molecular viscosity
$\nu$
. Moreover,
$\tilde {S}$
is the modified vorticity magnitude. Within the PMoE framework, we demonstrate the integration of three distinct expert formulations.
Type I: parameter correction for wall-attached flows. For equilibrium boundary layers where the baseline model structure is sound but parameter calibration is suboptimal, the expert provides a scalar correction to existing coefficients. For the wall-attached expert
$E_1$
, we employ a symbolic regression approach (Weatheritt & Sandberg Reference Weatheritt and Sandberg2016) to recalibrate the damping function
$f_{v1}$
as suggested by Bin et al. (Reference Bin, Huang and Yang2023). The expert predicts a corrected functional form that replaces the standard definition in (2.11) to better capture the logarithmic layer velocity profile.
Type II: correction of production term for flow separation. For flows dominated by separation and strong pressure gradients, the equilibrium assumption in the linear eddy viscosity hypothesis often fails. In such regions, conventional models may mischaracterise the eddy viscosity distribution near separation and within recirculation zones, which in turn affects the prediction of separation and reattachment. This behaviour is associated with the difficulty of linear eddy viscosity models in representing the non-equilibrium balance between production and dissipation under adverse pressure gradients and rapidly varying strain rates. Here, the expert
$E_2$
is formulated as a neural network derived from the field inversion and machine learning (FIML) method. This method is chosen because it has been widely applied and demonstrated success in separation corrections (Parish & Duraisamy Reference Parish and Duraisamy2016), produces spatially varying correction fields that are well suited for flows where modelling errors are localised, and is compatible with standard RANS solvers. It imposes a spatially varying correction term
$\beta$
into the production term of the transport equation:
\begin{equation} \mathcal{P}^* = \beta (\boldsymbol{x}) \boldsymbol{\cdot }C_{b1} \tilde {S} \tilde {\nu }. \end{equation}
This multiplicative correction allows the model to locally suppress or enhance turbulence production in recirculation zones without altering the model behaviour in the free stream.
Type III: constitutive relation modification for secondary flows. For corner flows where linear eddy viscosity models fail to predict secondary motions due to the isotropy assumption, which constrains Reynolds stress anisotropy and suppresses cross-stream momentum redistribution, the expert
$E_3$
modifies the stress–strain relationship itself. We adopt a data-driven calibration of the quadratic constitutive relation (QCR) (Spalart Reference Spalart2000):
\begin{equation} \tau _{ij,\textit{QCR}}=\tau _{\textit{ij}}^{linear} - C_{\textit{cr1}}(\boldsymbol{x})\left (O_{ik} \tau _{jk}+O_{jk} \tau _{ik}\right ), \end{equation}
where
$O_{ik}$
is the normalised rotation tensor. Unlike the standard QCR where
$C_{\textit{cr1}}$
is constant, the expert
$E_3$
predicts a spatially varying
$C_{\textit{cr1}}$
field to accurately capture corner vortices.
This modularity allows the PMoE framework to progressively integrate experts ranging from concise symbolic expressions to complex neural mappings without structural conflict, effectively decoupling the model architecture from the constraints of a single baseline formulation. This flexibility allows developers to select the most suitable expert design strategy for a given flow regime, and to update or expand the expert group in step with advances in the literature, thereby further enhancing the practicality, extensibility and efficiency of the framework.
While the formulation of individual experts reflects physical modelling considerations, the assignment of flow regimes during training and deployment is entirely feature-driven. No manual flow-type specification is required when processing a new case. The sampling strategy is likewise based on local kinematic indicators, ensuring minimal dependence on prior physical labelling.
3. Continual learning PMoE model with various cases introduced progressively
To demonstrate the sustainable generalisation capability of the PMoE framework, we employ a curriculum learning strategy. As shown in figure 4, each dataset is introduced sequentially on a case-by-case basis in order of increasing geometric complexity, starting from the aerofoil wake case (free-shear turbulence away from solid walls) and progressing to the channel flow cases (equilibrium wall-bounded flow with parallel flat wall), periodic hill cases (curved wall featuring separation), and the square duct case (featuring corner-induced secondary flow). It should be noted that the identification of whether a new dataset belongs to a known or novel regime is performed automatically by the router via the reconstruction-error-based confidence metric. However, the order in which datasets are introduced, and the design of the corresponding expert formulation at each stage, are prescribed by the user, analogous to a practitioner sequentially encountering new flow scenarios in an industrial workflow. This obviously relies on a modeller’s experience and physical insight. Importantly, because each router component and expert is trained independently with all previous modules frozen, the final model is invariant to the order of regime introduction.
Schematic diagrams of typical examples of various flow regimes: (a) 2DANW, (b) fully developed channel flow, (c) periodic hill flow, (d) square duct flow.
Figure 4. Long description
The image contains four schematic diagrams illustrating different flow regimes. Diagram (a) shows a 2DANW flow with a streamlined shape moving through a fluid. Diagram (b) depicts a fully developed channel flow with turbulent structures. Diagram (c) illustrates a periodic hill flow, highlighting the general flow direction and a recirculation zone. Diagram (d) presents a square duct flow with arrows indicating the direction of the inlet flow and the three-dimensional structure of the duct.
The evolution of the model is denoted by stages
$S_0$
to
$S_3$
, with the corresponding trained PMoE denoted as PMoE-S0, PMoE-S1, PMoE-S2 and PMoE-S3, respectively. At each stage, the router evaluates the novelty of the incoming data. We extract features from sampled points in the baseline RANS flow field, and construct input vector
$\boldsymbol{x}$
and matrix
${\unicode{x1D653}}\in \mathbb{R}^{N\times M}$
, where
$M=7,\ N=10\,000$
in our training process. The input data are then evaluated through the existing router modules. If the maximum confidence
$p_{\mathcal{K}}$
falls below the acceptance threshold
$T_{\textit{accept}} = 90\,\%$
, then the system triggers the expansion mechanism: a new autoencoder component is trained to recognise the regime, and a specialised expert is trained to close the RANS equations. The specific datasets and their role in the curriculum are detailed in table 2. It should be noted that each router (and expert model) is designed to handle a range of flow conditions over the same geometry. Thus for each flow regime, only some of the working conditions are used for the training, while the rest are left for model validation in § 4. To demonstrate the feasibility of the PMoE framework while maintaining minimal computational cost, the current training dataset is intentionally kept as small as possible. Further reducing the training set will inevitably degrade the model’s generalisability. For each flow regime, the router and the corresponding expert are trained using the same subset of cases. The numerical set-up of RANS calculation for each case is summarised in Appendix B.
Summary of flow cases used for training and validation.
Table 2. Long description
The table presents a summary of flow cases used for training and validation, organized into four stages. Each stage includes specific flow regimes, cases, abbreviations, stages, and references. Stage 0 focuses on the wake of a two-dimensional aerofoil, with the case 2DANW used for training, referenced from Nakayama 1985. Stage 1 involves a fully developed channel at different Reynolds numbers: C2000 for training, C5200 for training, and C8000 for validation, with references from Lee & Moser 2015 and Yamamoto & Tsuji 2018. Stage 2 examines periodic hills with varying alpha values: PH0p8 for validation, PH1p0 for training, PH1p2 for validation, and PH1p5 for validation, referenced from Xiao et al. 2020. Stage 3 covers square ducts at different Reynolds numbers: SD2500 for training, SD3500 for validation, and SD5693 for training, with references from Pinelli et al. 2010 and Vinuesa et al. 2018. The table provides a structured overview of the flow cases and their roles in the training and validation processes.
While completing the expansion of the router from PMoE-S0 to PMoE-S3, we also need to synchronously train the corresponding expert model for each newly identified flow regime. It should be particularly pointed out that the regimes are distinguished by the router automatically rather than pre-specified by humans during the training process.
3.1. Stage 0: baseline initialisation
The initial model, PMoE-S0, is established using the two-dimensional (2-D) aerofoil near wake (2DANW) case. This case is a verification of the aerofoil wake flow, with a uniform inflow, no-slip boundary conditions on the aerofoil surface, and a far-field free-stream condition. Detailed case configurations are referred to the NASA web page (data available online at https://tmbwg.github.io/turbmodels/airfoilwakeverif500c.html), and the grid used for the present RANS calculations is given in Appendix B. This regime represents simple free-shear turbulence where the standard SA model is known to perform adequately (Bin et al. Reference Bin, Huang and Yang2023). Consequently, we designate the SA model as baseline expert
$E_0$
, and its performance in the 2DANW case is presented in figure 5(a). Furthermore, the first router component
$C_0$
is trained on the wake features, establishing the baseline latent space distribution. As shown in table 3, the router identifies this regime with
$99.9\,\%$
confidence, serving as the anchor for future expansions.
Velocity profiles and contours predicted by the baseline SA model, the expert model and high-fidelity data: (a) expert for S0 trained by 2DANW; (b) expert for S1 trained by C5200; (c) expert for S2 trained by PH1p0; (d) expert for S3 trained by SD2500 and SD5693.
Figure 5. Long description
The image contains four subplots labeled (a), (b), (c), and (d), each showing different velocity profiles and contours. Subplot (a) compares experimental data with expert model predictions for S0 trained by 2DANW. Subplot (b) shows a graph comparing DNS, Baseline, and Expert models for S1 trained by C5200. Subplot (c) displays contour plots of the magnitude of velocity for DNS, Baseline, and Expert models for S2 trained by PH1p0. Subplot (d) presents contour plots of velocity components and magnitude for DNS, Baseline, and Expert models for S3 trained by SD2500 and SD5693. The graphs illustrate the differences in velocity predictions between the models and high-fidelity data.
Progressive construction of the PMoE model.
Table 3. Long description
The table presents the progressive construction of the PMoE model, detailing the introduction of new regimes and the confidence levels of new router components. It consists of four columns and four rows, including headers. The columns are labeled as New regime, MAX confidence with existing router, Introduced expert, and Confidence of new router component. The rows detail specific regimes such as 2DANW, C2000, C5200, PH1p0, SD2500, and SD5693, along with their respective confidence percentages. The table highlights the performance of different models and their confidence levels in identifying new regimes, with the SA baseline model serving as the initial expert. Notable trends include high confidence levels for the new router components across different regimes, indicating successful identification and integration of new regimes into the PMoE model.
3.2. Stage 1: adaptation to wall-bounded turbulence
In the second phase, the model is exposed to fully developed turbulent channel flows at
$Re_\tau = 2000$
and
$5200$
, which are driven by a volumetric body force, with cyclic boundary conditions at the inlet and outlet and no-slip walls. The grid sizes used for the present RANS calculations are provided in Appendix B. When these data are fed into PMoE-S0, the router yields confidence
$0.00\,\%$
as presented in table 3, correctly identifying that the physics of the wall turbulence differs fundamentally from the free-shear wake.
This triggers the creation of a new expert
$E_1$
. As discussed in § 2.2.3, we employ offline symbolic regression to derive a correction for the damping function
$f_{\nu 1}$
. The resulting expert formulation is
\begin{equation} f_{\nu 1}^* = 0.985 \tanh \left (\frac {\chi ^{1.3}}{31.1}\right ). \end{equation}
Applying this expert to RANS calculations, the corrected model is able to provide more accurate wall-normal profiles compared to the baseline SA model, agreeing well with the direct numerical simulations (DNS) result (Lee & Moser Reference Lee and Moser2015) as shown in figure 5(b). Simultaneously, a new router component
$C_1$
is trained. The updated model, PMoE-S1, subsequently recognises the channel data with
$99.5\,\%$
confidence, demonstrating its ability to capture wall-bounded flow characteristics (as represented here by channel flow), while retaining expert 0 for wakes.
3.3. Stage 2: capturing flow separation
Most industrial flows exhibit features such as strong pressure gradients, streamline curvature and separation, which violate the equilibrium assumption underlying conventional RANS models. The periodic hill flow is usually adopted as a canonical benchmark for separated turbulence, as it combines adverse pressure gradients, separation and reattachment in a simple, well-defined geometry, with the accurate prediction of the separation bubble remaining a key modelling challenge. The Reynolds number of the periodic hill cases is set as
$5600$
, consistent with DNS simulation (Xiao et al. Reference Xiao, Wu, Laizet and Duan2020), and the hill slope defined by
$\alpha = l/h$
varies across cases, where
$l$
and
$h$
are the hill width and height, respectively.
The curriculum thus proceeds to the periodic hill case with
$\alpha =1.0$
as shown in table 2. The PMoE-S1 router produces maximum confidence
$58.4\,\%$
as presented in table 3. Note that this relatively high level of confidence is likely due to shared wall-bounded features. Nevertheless, it falls below the
$90\,\%$
threshold, flagging the regime as novel.
To address the non-equilibrium physics in the recirculation bubble, we train a new expert
$E_2$
using the FIML method. This expert injects a source term
$\beta (\boldsymbol{x})$
into the production term as presented in (2.12). The details are elaborated in Appendix C. As shown in figure 5(c), the trained separation flow expert can more accurately capture the phenomena of separation and reattachment, while the estimation of the separation bubble size by the basic SA model is significantly larger than that of the expert model and DNS results. Furthermore, upon integrating
$E_2$
and its corresponding router component, the identification confidence for separated flows (as represented here by periodic hill flow) rises to
$93.5\,\%$
, yielding model PMoE-S2.
3.4. Stage 3: modelling corner-induced anisotropy
Finally, the model encounters the square duct flow at
$Re=2500, 5693$
. This case is driven by a volumetric body force, with cyclic inlet and outlet boundaries and no-slip walls on all duct surfaces. This case is physically distinct due to the secondary motions in the corner region, which are impossible for linear eddy viscosity models. The training data are obtained from published DNS simulations (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010; Vinuesa, Schlatter & Nagib Reference Vinuesa, Schlatter and Nagib2018). The Reynolds number is based on the edge half-length
$h=D/2$
of the square duct and the bulk velocity
$U_b$
as shown in figure 4(d).
The PMoE-S2 router yields confidence
$56.7\,\%$
as presented in table 3, again triggering expansion. A new expert
$E_3$
is then trained to predict the nonlinear stress coefficient
$C_{\textit{cr1}}(\boldsymbol{x})$
in the QCR formulation in (2.13). The details are elaborated in Appendix C. As shown in figure 5(d), the trained expert model successfully captures the generation of secondary flow, while the baseline SA model is unable to calculate the secondary flow. The final model, PMoE-S3, achieves recognition confidence
$97.0\,\%$
for this regime (as represented here by square duct flow).
4. Model validation
The validation of the PMoE framework focuses on three critical performance metrics: the interpretability of the unsupervised routing, the prevention of catastrophic forgetting on previously learned regimes, and the generalisation capability to cases with different operating conditions or configurations. All results presented herein utilise the final PMoE-S3 model to demonstrate its cumulative capabilities.
4.1. Validation for the PMoE router
A prerequisite for reliable expert selection is the router’s ability to distinguish flow regimes based on local RANS features. We analyse the four-dimensional latent space of the autoencoder components using the Mahalanobis distance (Mahalanobis Reference Mahalanobis1936)
$D_M$
, which measures the distance between a point
$\boldsymbol{z}$
and a distribution with mean
$\boldsymbol{\mu }$
and covariance
$\boldsymbol{\varSigma }$
:
\begin{equation} D_M^2(\boldsymbol{z}) = (\boldsymbol{z} - \boldsymbol{\mu })^{\rm T} \boldsymbol{\varSigma }^{-1} (\boldsymbol{z} - \boldsymbol{\mu }). \end{equation}
Unlike the Euclidean distance,
$D_M$
accounts for the correlation structure of the learned representations. Based on this metric, we compute both the intra-cluster dispersion, which quantifies the compactness of samples within a regime, and the inter-cluster dispersion measuring the separation between regime centroids, as illustrated in figure 6. The ratio between the averaged inter- and intra-cluster dispersions is adopted as a diagnostic indicator of the latent-space organisation. A high ratio implies that different flow regimes are topologically distinct within the component’s latent space, whereas a low ratio suggests possible overlap.
Schematic illustration of intra- and inter-cluster dispersion based on the Mahalanobis distance.
Figure 6. Long description
The diagram shows four clusters labeled Cluster 0, Cluster 1, Cluster 2, and Cluster 3. Intra-cluster dispersion is represented by arrows within each cluster, indicating the spread of data points within the same cluster. Inter-cluster dispersion is shown by arrows pointing between different clusters, indicating the distance between the centers of different clusters. The diagram also contrasts poorly dispersed and well-dispersed clusters, illustrating the difference in data point distribution.
This analysis is carried out for all four components. Each component is exposed to data from all four flow regimes to assess its discriminative behaviour, and table 4 summarises these statistics. The intra-cluster dispersion remains consistent as
$\approx 1.8$
across all components, indicating that the autoencoders compress different flow regimes with comparable efficiency. Furthermore, the ratios between the averaged inter- and intra-cluster dispersions are also presented, and all four components exhibit ratios exceeding
$3.4$
, indicating that the corresponding clusters have little overlap as illustrated by the well-separated configuration in figure 6. Although the router decision is ultimately based on the reconstruction error along the full encoder–decoder pathway, the latent-space analysis still provides an interpretable metric for its behaviour, as the observed clustering patterns reflect the feature representations that govern reconstruction quality and hence routing outcomes. This mechanism, combined with the competitive winner-takes-all strategy, ensures robust classification. As evidenced in table 5, the router achieves
$\gt \kern-2pt 90\,\%$
classification accuracy on the validation sets, confirming that the present router successfully distinguishes the tested regimes.
Dispersion analysis of the latent space clusters using Mahalanobis distance.
Table 4. Long description
The table presents dispersion analysis of latent space clusters using Mahalanobis distance for four flow regimes across four components: Wake, Channel, Periodic hill, and Square duct. It consists of four rows and four columns. The columns are labeled as Dispersion, Average intra-cluster, Average inter-cluster, and Inter-cluster/intra-cluster. The rows are labeled with the four components. The table shows the average intra-cluster dispersion values for each component: 1.8207 for Wake, 1.8383 for Channel, 1.8438 for Periodic hill, and 1.6947 for Square duct. The average inter-cluster dispersion values are: 6.6845 for Wake, 6.3698 for Channel, 7.6183 for Periodic hill, and 7.8923 for Square duct. The inter-cluster/intra-cluster dispersion ratios are: 3.6713 for Wake, 3.4650 for Channel, 4.1319 for Periodic hill, and 4.6569 for Square duct. The table indicates that the intra-cluster dispersion remains consistent across all components, suggesting that the autoencoders compress different flow regimes with comparable efficiency. The ratios between the averaged inter- and intra-cluster dispersions exceed one, indicating that the clusters have little overlap, as illustrated by the well-separated configuration in figure 6.
Confidence distribution of the PMoE-S3 router over different components for all flow cases.
Table 5. Long description
The table presents the confidence distribution of the PMoE-S3 router across various flow cases and components. It includes data for cases such as ANW, C2000, C5200, C8000, PH0p8, PH1p0, PH1p2, PH1p5, SD2500, SD3500, and SD5693. Each case is evaluated across several components: Unknown, Wake, Channel, Periodic hill, and Square duct. The table shows the percentage distribution for each component, highlighting the dominant component for each case. For instance, the ANW case has 99.7 percent confidence in the Wake component, while the SD2500 case shows 99.3 percent confidence in the Square duct component. The data indicates the router's ability to distinguish between different flow regimes effectively.
Boldface in the table highlights confidence values that exceed the threshold
$T_{\textit{accept}}$
of
$90\,\%$
.
4.2. Prevention of catastrophic forgetting via in-distribution tests
A primary failure mode in continual learning is catastrophic forgetting, where optimising for new tasks degrades performance on previous ones (French Reference French1999). To assess whether the current PMoE framework is able to avoid catastrophic forgetting, we apply the final PMoE-S3 model back to the initial training cases. The corresponding results are presented in figure 7. The number shown in the top-left corner of each panel denotes the relative root mean square error (RMSE) reduction for the PMoE model
$\textit{RMSE}_{\textit{PMoE}}$
with respect to the baseline
$\textit{RMSE}_{\textit{base}}$
, defined as
\begin{equation} \mathcal{R} = (\textit{RMSE}_{\textit{base}} - \textit{RMSE}_{\textit{PMoE}})/\textit{RMSE}_{\textit{base}}, \end{equation}
where
$\textit{RMSE} = \sqrt {({1}/{N})\sum _{i=1}^{N} (u_i - u_i^{\textit{ref}} )^2}$
. Here,
$u_i$
denotes the model prediction (baseline SA or PMoE, respectively) at the
$i$
th reference data location, the error is evaluated for the velocity component shown in each panel, and
${u_i^{\textit{ref}}},\ i=1,\ldots,N$
, are the corresponding reference values, with model predictions interpolated onto these locations when necessary. No value is shown in in figure 7(a) as expert
$E_0$
coincides with the baseline.
Mean velocity profiles for trained cases: (a) 2DANW (colours indicate streamwise locations), (b) C2000, (c) C5200, (d,e) PH1p0, (f–h) SD2500, (i–k) SD5693. Circles denote DNS/experimental data; green squares denote baseline SA; solid red lines denote PMoE-S3. The number shown in the top left corner of each panel is the relative error defined by (4.2).
Figure 7. Long description
The line graph displays mean velocity profiles for trained cases across different datasets and models. The subplots (a) through (k) represent different cases: (a) 2DANW with streamwise locations indicated by colors, (b) C2000, (c) C5200, (d) and (e) PH1p0, (f) to (h) SD2500, and (i) to (k) SD5693. Each subplot includes circles denoting DNS/experimental data, green squares representing baseline SA, and solid red lines indicating PMoE-S3. The relative error for each case is shown in the top left corner of each panel. The x-axis represents normalized velocity (u/U_inf or U/U_b), and the y-axis represents normalized height (y/c or y/h). All values are approximated.
Figure 7(a) presents the wake velocity profiles for the 2DANW case (stage 0), with each colour representing a different observation position: black for
$x/c=1.01$
, red for
$x/c=1.05$
, green for
$x/c=1.20$
, blue for
$x/c=1.40$
, yellow for
$x/c=1.80$
, and brown for
$x/c=2.19$
. Despite the model undergoing three subsequent expansions involving wall-bounded and secondary flows, the prediction for the free-shear wake remains virtually identical to the baseline SA model, and agrees well with experimental data. This suggests that in this near-wake region, the flow is dominated by the trailing-edge shear layer and vortex shedding, which are largely insensitive to wall effects. This confirms that the modular router successfully shields the baseline expert
$E_0$
from interference.
Similarly, figures 7(b) and 7(c) show the mean velocity profiles extracted from the channel flow cases at
$Re_\tau =2000$
and
$Re_\tau =5200$
(stage 1), while figures 7(d) and 7(e) present the periodic hill at
$\alpha =1.0$
(stage 2). The final PMoE-S3 model retains the specialised accuracy of experts
$E_1$
and
$E_2$
in according tests. Particularly in the periodic hill case, the reattachment locations are quantified in table 6. For the training PH1p0 case, the baseline model predicts delayed reattachments near
$x/h=7.31$
, whereas the PMoE model provides more accurate predictions for the separation bubble size with reattachment point
$x/h=5.87$
. This improvement reflects the model’s ability to capture the enhanced turbulence production and shear-layer dynamics responsible for earlier flow reattachment.
Reattachment locations of periodic hill cases.
Table 6. Long description
The table presents reattachment locations for periodic hill cases, comparing three models: DNS, Baseline, and PMoE-S3. It includes four cases: PH0p8, PH1p0, PH1p2, and PH1p5. The DNS column shows values 5.21, 5.04, 4.50, and 4.21. The Baseline column shows values 7.34, 7.31, 7.88, and 8.32. The PMoE-S3 column shows values 7.19, 5.87, 5.36, and 4.71. The table highlights the differences in reattachment locations predicted by each model.
The capability of the final model PMoE-S3 for a trained square duct (
$Re = 2500$
) is shown in figures 7(f–h). Figures 7(f) and 7(g) show the velocity profiles at
$z/h=0.3, h=D/2$
, while figure 7(h) represents
$y/h=0.25, 0.5, 0.75, 1.0$
from left to right just as the orange arrows shown in figure 4(d). The square duct case contains typical corner-induced secondary flow, which cannot be characterised by the SA model based on linear eddy viscosity. The results of PMoE-S3, however, have a significant improvement compared to the baseline model. This suggests that the trained PMoE-S3 model can effectively represent the secondary flow in the non-flow direction due to the introduction of nonlinear correction based on QCR. Moreover, figures 7(i–k) present the results of the SD5693 case. The performance of PMoE-S3 remains excellent, with significant improvements in both flow velocity and secondary flow representation compared to the baseline.
4.3. Generalisation to cases with different operating conditions
The robustness of the PMoE-S3 model is further tested against cases with different operating conditions never seen during training. The mean velocity profiles of all the cases are shown in figure 8. Compared with the cases in § 4.2, these new tests evaluate the model’s generalisation ability to different Reynolds numbers and periodic hill geometries.
Mean velocity profiles for unseen cases with different operating conditions: (a) C8000, (b) PH0p8, (d) PH1p2, (f) PH1p5, and (c,e,g,h) SD3500. Circles denote DNS/experimental data; green squares denote baseline SA; solid red lines denote PMoE-S3. The number shown in the top left corner of each panel is the relative error defined by (4.2).
Figure 8. Long description
The image contains eight subplots labeled (a) through (h), each showing mean velocity profiles for unseen cases with different operating conditions. Each subplot compares three datasets: DNS/experimental data represented by circles, baseline SA represented by green squares, and PMoE-S3 represented by solid red lines. Subplot (a) shows data for C8000, (b) for PH0p8, (d) for PH1p2, (f) for PH1p5, and (c), (e), (g), and (h) for SD3500. The x-axis represents normalized velocity (U+ or Ux/Ub), and the y-axis represents normalized height (y+ or y/h). Each subplot includes a relative error percentage in the top left corner, indicating the accuracy of the PMoE-S3 model compared to the baseline SA model. The trends show how well the PMoE-S3 model fits the DNS/experimental data across different operating conditions. All values are approximated.
The model is applied to a channel flow at
$Re_\tau =8000$
and a square duct at
$Re=3500$
. In both cases, the Reynolds numbers are different from the specific training conditions. The router successfully maps these flows to experts
$E_1$
and
$E_3$
, respectively, as shown in table 5. The velocity profiles show excellent agreement with DNS data, indicating that the learned corrections, including the damping function modification in
$E_1$
and the QCR coefficient field in
$E_3$
, capture the underlying physics rather than merely overfitting the training data.
We further test the periodic hill flows with varying slopes
$\alpha = 0.8, 1.2, 1.5$
. These geometric changes significantly alter the adverse pressure gradient and the separation bubble size. As shown in figures 8(b), 8(d) and 8(f), the PMoE-S3 model consistently outperforms the baseline, despite the model being trained only for the
$\alpha =1.0$
case. Similar to the conclusion drawn from § 4.2, the baseline SA model overestimates the size of the separation bubble in the lower hill area, while the PMoE-S3 model shows much better agreement with the result of DNS as shown in table 6. This suggests that the FIML-trained expert
$E_2$
has learned a generalised correction for separation control that is robust to geometric deformation.
Up to now, the PMoE-S3 model has shown satisfactory predictive accuracy in a series of testing cases with different operating conditions, successfully overcoming catastrophic forgetting. These results indicate that although each expert model is trained for one or two operating conditions, it can actually handle a relatively wide range of operating conditions within the same geometry. For extreme operating conditions that exceed the generalisation boundaries of the expert models, however, the PMoE framework can be extended by fine-tuning only the relevant expert and its router, a far more efficient approach than retraining a monolithic model, highlighting its practical advantage.
4.4. Computational efficiency
A key practical consideration for the PMoE framework is its computational cost during inference. The complete inference pipeline consists of three stages: (i) a baseline SA simulation, from which input features are extracted; (ii) autoencoder-based routing; and (iii) a corrected RANS simulation using the activated expert. In many industrial workflows, a baseline RANS solution is already available as part of the standard design or analysis process. In such cases, only stages (ii) and (iii) constitute additional cost. When the router assigns the case to expert
$E_0$
(the standard SA model), the baseline solution is itself the final prediction, and no additional simulation is required.
To quantify the computational overhead, we performed wall-clock timing tests using both the baseline SA model and the PMoE-S3 model on identical hardware. (All tests were conducted on a system equipped with an
$\textrm {Intel}^{\circledR }$
$\textrm {Xeon}^{\circledR }$
Gold 6530 CPU and an
$\textrm {NVIDIA}^{\circledR }$
AD102 (GeForce RTX 4090) GPU.) The results are summarised in table 7. Across all tested cases, the routing time remains below 10 s, confirming that the autoencoder-based router is sufficiently lightweight for its overhead to be negligible. The full-pipeline overhead, which accounts for the baseline SA run, ranges from approximately
$100\,\%$
to
$104\,\%$
, reflecting the cost of effectively running RANS twice. This represents the worst-case scenario in which no prior baseline solution exists.
Comparison of computational wall-clock time between the baseline SA and PMoE-S3 models. Here,
$T_{\textit{tot}} = T_{\textit{base}} + T_{\textit{route}} + T_{\textit{exp}}$
.
Table 7. Long description
The table presents a comparison of computational wall-clock time between the baseline SA and PMoE-S3 models for different cases. It includes columns for baseline simulation time, routing time, expert simulation time, total time, and the increase ratio. The cases listed are ANW, C5200, PHI1p5, and SD5693. For each case, the baseline simulation time, routing time, expert simulation time, total time, and increase ratio are provided. The routing time remains below 10 seconds for all cases, indicating a negligible overhead. The total time varies significantly across cases, reflecting the computational cost of running RANS twice in the worst-case scenario.
Importantly, owing to the sparse activation mechanism, model expansion through continual learning does not lead to a proportional increase in inference cost. Regardless of the total number of experts in the PMoE model, only one expert is activated per case, ensuring that the per-case simulation cost scales with the size of a single expert rather than the entire ensemble.
5. Discussion
The results in § 4 have established that the PMoE framework achieves accurate predictions across canonical flows with different operating conditions while effectively avoiding catastrophic forgetting through its modular design. Two key questions arise when considering practical deployment: (i) how the model generalises to flow configurations that differ from the training set not only in operating conditions but also in geometry and boundary condition type; and (ii) how the current routing strategy performs when multiple flow regimes coexist within a single computational domain, and what alternative strategies may be considered. To address these questions, § 5.1 presents additional test cases with geometrically distinct configurations, and § 5.2 investigates routing strategies for multi-regime flows using a representative three-dimensional (3-D) configuration.
5.1. Generalisation to geometrically distinct configurations
The PMoE framework’s generalisation capability arises from two complementary mechanisms: the router, which identifies the most relevant expert for a given flow based on learned feature-space representations, and the intrinsic generalisation of the expert models themselves, which encode physically meaningful corrections rather than case-specific data fitting. To probe the limits of both mechanisms, we introduce three additional test cases spanning distinct geometrical configurations, as shown in figure 9, including the 2-D NASA wall-mounted hump separated flow (2DWMH), the curved backward-facing step at
$Re = 13\,700$
(CBFS), and the rectangular duct flow at
$Re = 2500$
with aspect ratio 3 (RD2500). These cases are deliberately chosen to test the framework on geometries and boundary conditions that were not encountered during training, while retaining physical mechanisms, i.e. separation under adverse pressure gradients and corner-induced secondary flow, that overlap to varying degrees with the trained regimes. Notably, the 2DWMH case also incorporates an upstream zero-pressure-gradient boundary layer developing beneath a free stream, thereby simultaneously testing the coupling between free-stream, wall-attached and separated flow regions within a single configuration.
Schematic diagrams of additional test cases with different configurations: (a) 2DWMH, (b) CBFS, (c) RD2500.
Figure 9. Long description
The diagram presents three distinct test cases for computational fluid dynamics (CFD) simulations. The first configuration, labeled 2DWMH, shows a two-dimensional wall model with a highlighted recirculation zone and general flow direction indicated by an arrow. The second configuration, labeled CBFS, depicts a solid wall with a recirculation zone and general flow direction. The third configuration, labeled RD2500, illustrates a three-dimensional rectangular duct with inlet flow, dimensions labeled B and D, and internal flow patterns indicated by arrows. Each subfigure demonstrates different flow behaviors and configurations used in CFD testing.
5.1.1. Router discrimination on unseen geometries
Table 8 summarises the confidence distribution of the PMoE-S3 router for the three additional cases. The confidence levels span a considerable range, from a maximum of
$90.8\,\%$
for RD2500 to
$54.7\,\%$
for the most challenging case (2DWMH). These differences can be understood by examining how the geometric and physical characteristics of each new case align with the feature signatures learned during training, as characterised by the PFI analysis in Appendix A.
Confidence distribution of the PMoE-S3 router over different components for additional flow cases.
Table 8. Long description
The table presents the confidence distribution of the PMoE-S3 router for three additional cases, with percentages across different components. It includes five columns: Unknown, Wake, Channel, Periodic hill, and Square duct, each representing a different confidence level. The table has three rows, each corresponding to a different case: 2DWMH, CBFS, and RD2500. For the 2DWMH case, the percentages are 10.0 for Unknown, 0.0 for Wake, 35.2 for Channel, 54.7 for Periodic hill, and 0.0 for Square duct. For the CBFS case, the percentages are 0.9 for Unknown, 0.0 for Wake, 22.4 for Channel, 76.7 for Periodic hill, and 0.0 for Square duct. For the RD2500 case, the percentages are 0.0 for Unknown, 0.0 for Wake, 1.8 for Channel, 7.4 for Periodic hill, and 90.8 for Square duct. The data highlights significant variations in confidence levels across different cases and components.
Boldface in table indicates the highest router confidence of all components.
The RD2500 rectangular duct preserves the most similar features of the training square duct flow. It is a fully developed, pressure-driven wall-bounded flow with no-slip conditions on all walls and corner-induced secondary motions generated by Reynolds stress anisotropy. The only systematic deviation is the aspect ratio, which changes from 1 to 3, modifying the relative strength and spatial extent of the corner vortices. As a result, the flow retains the same wall-driven shear representation, and the router assigns the highest confidence (
$90.8\,\%$
) to the square duct class
$C_3$
.
The CBFS case shares the dominant physics of the periodic hill training case, an adverse-pressure-gradient-driving shear-layer separation followed by reattachment, but introduces a sharp-to-curved geometric expansion that alters the pressure-gradient distribution and shear-layer inception. Although the geometric expansion differs, the resulting flow structure remains governed by pressure–shear interaction. Consequently, the feature representation remains close to that of the periodic hill, and the router identifies class
$C_2$
with confidence
$76.7\,\%$
, reflecting a substantial overlap in the underlying flow dynamics.
The 2DWMH case likewise involves separation driven by an adverse pressure gradient, yet the separation occurs over a convex hump surface in an external-flow setting rather than within a confined channel, introducing greater geometric and boundary condition dissimilarity relative to the periodic hill. In particular, the presence of an upstream attached boundary layer and the absence of confinement lead to a more mixed flow character. This reduces the similarity in the encoded representation, resulting in a lower router confidence of
$54.7\,\%$
for
$C_2$
, with a secondary attribution of
$35.2\,\%$
to the channel class
$C_1$
.
The pattern emerging from table 8 is physically intuitive. Flows that retain stronger similarity to the training regimes yield higher router confidence, whereas configurations with more pronounced geometric and boundary condition deviations exhibit lower and more distributed confidence.
5.1.2. Expert prediction
Beyond the router’s discriminative capability, the expert models themselves exhibit strong generalisation when activated on the basis of the highest confidence assignment. For the CBFS case, the periodic hill expert
$E_2$
substantially reduces the overprediction of the separation bubble extent produced by the baseline SA model, yielding corrected velocity profiles that faithfully capture the interaction between the adverse pressure gradient and shear-layer growth, and accurately locate the reattachment point as shown in figure 10(a). The reference data are obtained from the LES of Bentaleb, Lardeau & Leschziner (Reference Bentaleb, Lardeau and Leschziner2012). The underlying correction, a recalibration of the non-equilibrium balance between turbulence production and dissipation under adverse pressure gradients, transfers effectively to this geometrically related but distinct separated-flow configuration.
Predictions obtained with the highest-confidence experts for additional test cases: (a) CBFS, (b) 2DWMH, (c–e) RD2500. Circles denote DNS/experimental data; green squares denote baseline SA; solid red lines denote expert model. The number shown in the top left corner of each panel is the relative error defined by (4.2).
Figure 10. Long description
The image contains five graphs showing predictions for additional test cases using different models. The first graph (a) displays the CBFS model with circles denoting DNS/experimental data, green squares for baseline SA, and a solid red line for the expert model. The second graph (b) shows the 2DWMH model with similar notations. The third and fourth graphs (c and d) present the RD2500 model with different variables on the x-axis. The fifth graph (e) combines data from different models. Each graph includes a relative error percentage in the top left corner. The circles, squares, and lines represent different data sources and models, with the red line indicating the expert model's predictions. All values are approximated.
The 2DWMH case provides an even more stringent test: although the router confidence is appreciably lower, activating
$E_2$
still yields improved pressure coefficient (figure 10
b) distributions, with particular gains in the separated region and the downstream recovery zone where the baseline model fails to predict the correct separation onset and pressure recovery. The reference data are obtained from the experiment of Naughton, Viken & Greenblatt (Reference Naughton, Viken and Greenblatt2006) and Greenblatt et al. (Reference Greenblatt, Paschal, Yao and Harris2006a
,
Reference Greenblatt, Paschal, Yao, Harris, Schaeffler and Washburnb
). This result demonstrates that the FIML-trained expert encodes transferable corrections for the fundamental mechanisms governing turbulent separation. Specifically, the suppression of excessive turbulence production in adverse-pressure-gradient regions that remain effective even when the external geometry and boundary layer state differ substantially from those of the training case.
Finally, for the RD2500 case, the square duct expert
$E_3$
, which predicts a spatially varying nonlinear stress coefficient
$C_{\textit{cr1}}(\boldsymbol{x})$
, accurately reproduces the corner-induced secondary vortices driven by anisotropic Reynolds stresses and cross-stream pressure gradients, and improves the streamwise velocity profile relative to the baseline SA prediction, despite the change in aspect ratio. This confirms that the expert has encoded the essential physics of corner-induced anisotropy rather than a geometry-specific correction. The reference data are obtained from the DNS of Vinuesa et al. (Reference Vinuesa, Schlatter and Nagib2018). The mean velocity profiles shown in figures 10(c) and 10 (d) are taken at
$y/h = 0.5,\ h=D/2$
, where
$D$
denotes the duct height, while those in figure 10(e) are extracted along the lines indicated in figure 9(c).
5.1.3. Synthesis
Collectively, these results demonstrate that the generalisation capability of the PMoE framework arises from the combined action of its two core components. The router provides a physically interpretable confidence measure that quantifies the relevance of each trained expert to an unseen flow, and serves as an indicator of similarity in the learned feature space. The expert models encode physically meaningful corrections rather than case-specific data fitting, and these corrections remain effective across variations in geometry and boundary conditions as long as the dominant physical mechanism is preserved. Even when the router confidence does not reach the acceptance threshold
$90\,\%$
, the expert with the highest confidence can still provide noticeable improvement over the baseline model, as demonstrated by the 2DWMH case. Conversely, when no expert is sufficiently relevant, the acceptance threshold
$T_{\textit{accept}}$
ensures that the framework defaults to the unmodified baseline model, avoiding the introduction of spurious corrections. This combination of a discriminative router and robust physics-based experts provides a degree of generalisation and fault tolerance that is not typically achieved by monolithic data-driven models. It also offers a practical strategy for extension. When the limit of an expert is reached for a substantially different flow configuration, only the relevant expert and its associated router component need to be refined, rather than retraining the entire model.
5.2. Routing strategies for more complex flows
5.2.1. Assumptions and limitations of the current routing strategy
The present PMoE router classifies flow types based on the baseline SA model’s predictions. The effectiveness generally rests on several key assumptions. Specifically, the flow under consideration can be adequately characterised by the features listed in table 1; and the baseline SA model provides a qualitatively correct flow field, from which meaningful features can be extracted. As demonstrated in the preceding cases, even when the SA model exhibits deviations, such as overestimating the separation bubble in the periodic hill case or failing to capture the secondary flows in the square duct, these inaccuracies have only a limited effect on the routing outcomes. This robustness stems from the design of the router, which considers a comprehensive set of normalised invariant features that characterise the local flow.
Conversely, the framework may become unreliable when these assumptions are violated. This can arise if the baseline model produces a qualitatively incorrect flow topology, such as failing to capture separation or transition, thereby corrupting the extracted features.
Furthermore, the results presented in §§ 4 and 5.1 employ a top-
$1$
(winner-takes-all) routing strategy, in which the expert with the highest router confidence is exclusively activated for the entire computational domain. This approach is effective when the flow field is dominated by a single physical regime, as demonstrated for the canonical test cases. For flows closely resembling the trained cases (e.g. RD2500), the corresponding expert can be directly activated for optimal corrections. For moderately similar flows (e.g. CBFS and 2DWMH), the router’s top recommendation provides effective corrections even at reduced confidence. For unknown cases identified by the PMoE router, the baseline SA model will be invoked, ensuring physically consistent results. Note that the top-1 routing is similar to industrial strategies such as GEKO (Menter & Matyushenko Reference Menter and Matyushenko2025). Applying a single expert globally ensures numerical stability. While a single expert is applied globally, the corrections that it produces are inherently locally determined by the pointwise input features, so the expert remains quiescent in regions where the local flow state does not trigger its correction (e.g.
$\beta (\boldsymbol{x}) \approx 1$
in equilibrium regions for expert
$E_2$
).
For complex industrial configurations, when multiple flow regimes coexist with comparable importance within a single domain, a single top-1 routing decision may no longer be sufficient to represent the flow physics. In such situations, additional reference data or more advanced routing strategies may be required to ensure reliable performance.
5.2.2. Multi-regime test: the 3-D diffuser
To examine the behaviour of the PMoE framework under conditions of genuinely coexisting flow regimes, we consider the 3-D diffuser case (Cherry, Elkins & Eaton Reference Cherry, Elkins and Eaton2008). As shown in figure 11(a), the geometry consists of a straight square duct inlet section that transitions into a rectangular expansion section, where one wall diverges at a prescribed angle while the opposite wall remains flat. The boundary conditions include a fully developed duct inflow, no-slip conditions on all walls, and a zero-gradient outflow at the diffuser exit. The upstream section features fully developed square duct flow with corner-induced secondary motions, while the expansion region introduces an adverse pressure gradient, streamline curvature, and ultimately 3-D flow separation from the diverging wall. The interplay between the corner vortices carried from the inlet section and the incipient separation renders this flow particularly challenging for RANS turbulence models (Abe & Ohtsuka Reference Abe and Ohtsuka2010).
Performance evaluation of the mixed model for the 3-D diffuser case. (a) Zero streamwise-velocity contours on the
$x$
–
$y$
plane at
$z=B/2$
(midplane), where black, green, and red denote DNS, baseline SA, and the PMoE experts, respectively. (b) Skin-friction coefficient (
$C_f$
) distribution along the bottom wall centreline (
$z=B/2$
), extracted along the yellow line (length
$L$
) indicated in (a). (c) Mean velocity profiles on the midplane. Circles denote DNS/experimental data; green squares denote baseline SA; solid red lines denote expert model.
Figure 11. Long description
The image contains three graphs analyzing the performance of the PMoE framework for a 3-D diffuser case. The first graph shows zero streamwise-velocity contours on the midplane, with black, green, and red lines representing DNS, baseline SA, and the PMoE experts, respectively. The second graph displays the skin-friction coefficient distribution along the bottom wall centerline, with circles denoting DNS/experimental data, green squares representing baseline SA, and solid red lines indicating the expert model. The third graph presents mean velocity profiles on the midplane, with circles for DNS/experimental data, green squares for baseline SA, and solid red lines for the expert model. The graphs illustrate the comparison of different models' performance in capturing flow characteristics.
This configuration is obviously beyond the scope of the current top-
$1$
strategy. The router identifies three competing regimes with comparable confidence: the periodic hill class
$C_2$
at
$40.7\,\%$
, the channel class
$C_1$
at
$36.5\,\%$
, and the square duct class
$C_3$
at
$21.5\,\%$
. No single component achieves a dominant confidence, reflecting the genuinely multi-physics nature of the case.
As a preliminary exploration of strategies beyond the top-
$1$
paradigm, we conducted a test in which all three major expert models (
$E_1$
,
$E_2$
and
$E_3$
) were simultaneously activated across the entire 3-D diffuser domain. This combined application exploits the inherently local nature of the expert corrections. Specifically,
$E_1$
primarily acts in near-wall regions through the damping function modification,
$E_2$
adjusts the production term via
$\beta (\boldsymbol{x})$
in zones of strong adverse pressure gradients and recirculation, and
$E_3$
introduces nonlinear stress corrections through
$C_{\textit{cr1}}(\boldsymbol{x})$
in corner regions. In flow regions that do not exhibit strong similarity to the training cases, i.e. where the local features do not trigger significant corrections, the predicted corrections remain small, and the model behaviour effectively reverts to that of the baseline SA model. As shown in figure 11(b), this simple combined strategy already yields encouraging results. In particular, the skin-friction coefficient
$C_f$
along the bottom wall shows clear improvement compared to the baseline SA prediction, and is notably closer to the DNS data of Ohlsson et al. (Reference Ohlsson, Schlatter, Fischer and Henningson2010). This demonstrates that the modular expert design, while developed for single-regime flows, has the potential to be meaningfully extended to multi-regime configurations.
It should be noted, however, that the 3-D diffuser involves a substantially higher degree of geometric and physical complexity than any individual training case, including 3-D corner-driven separation interacting with pressure-induced detachment. Consistently, the results reveal mixed performance: while the adverse-pressure-gradient effects near the bottom wall are better captured, as reflected by both the improved
$C_f$
distribution and the near-wall velocity profiles (figure 11
c). In contrast, the velocity profiles in the core region show larger deviations compared to the baseline SA prediction, indicating that the concurrent activation of multiple experts does not uniformly improve the solution throughout the flow field. Furthermore, the separation extent inferred from the
$u=0$
contour indicates that directly activating multiple experts shifts the separation location from the baseline’s overly downstream prediction to an overcorrected upstream position. The current expert library, designed and trained on canonical 2-D and simple 3-D configurations, is not expected to fully capture such interactions. The improvement observed here should therefore be viewed as a demonstration of the framework’s extensibility and of the transferability of the expert corrections, rather than as a claim that the current model is sufficient for this class of flows.
5.2.3. Toward region-aware routing
The 3-D diffuser results motivate several directions for extending the current routing strategy. A natural next step is region-wise application, in which the computational domain is partitioned into subregions and each subregion is assigned to the expert whose local confidence is highest. Conceptually, this approach exploits the spatial distribution of the router’s reconstruction error to create a patchwork of expert-corrected zones. This strategy has the potential to improve predictions in complex flows by locally matching the correction to the prevailing physics. However, it introduces practical challenges, including the treatment of interfaces between expert regions, where discontinuities in the correction field may arise and compromise solver convergence, and the sensitivity of the partitioning to the threshold and sampling density.
An alternative and potentially more robust approach is a top-
$k$
(
$k \gt 1$
) gating strategy with soft blending, in which the corrections from the
$k$
most confident experts are combined using continuous weights derived from the reconstruction-error-based confidence
$p_k$
. Because
$p_k$
varies spatially through its dependence on local features, such a strategy would enable smooth transitions between regimes without the hard boundaries inherent in region-wise partitioning. The modular autoencoder architecture of the PMoE router is inherently compatible with such extensions, as each component can independently evaluate reconstruction error at the point level.
It is important to note, however, that the current heterogeneous expert formulations, which range from symbolic expressions modifying
$f_{\nu 1}$
to neural-network-based corrections of different target quantities (
$\beta$
for production,
$C_{\textit{cr1}}$
for the constitutive relation), preclude straightforward soft blending, since the expert outputs are not directly commensurable. A future implementation targeting smooth region-wise routing would benefit from a unified expert format, e.g. all experts predicting the same correction type, such as a Reynolds stress anisotropy tensor or a common multiplicative field. This represents a fundamental design trade-off. The heterogeneous formulations adopted in the present work maximise physical fidelity by allowing each expert to employ the correction form best suited to its regime, whereas homogeneous formulations would sacrifice some of this flexibility in exchange for the numerical smoothness required for blending. Resolving this trade-off, and developing the associated interface-treatment and blending algorithms, will be a key focus of future work. We anticipate that such developments will enable improved performance not only for the 3-D diffuser case but also for other complex 3-D industrial flows in which multiple physical mechanisms interact.
6. Conclusions
In this work, we propose a progressive modular mixture-of-experts framework, termed PMoE, to address generalisation and continual learning in data-driven RANS turbulence modelling. The framework is based on a successively trainable autoencoder-based similarity router that enables the progressive introduction of specialised turbulence model experts as new flow regimes emerge. Using training data covering four representative regimes, namely aerofoil wakes, channel flows, periodic hill flows and square duct flows, we construct the PMoE-S3 model and demonstrate its ability to deliver accurate predictions for both seen and unseen test cases.
A key advantage of the proposed framework lies in its support for continual learning without catastrophic forgetting. New flow regimes can be incorporated by adding dedicated experts and router components, while previously trained components remain unchanged, preserving established predictive capabilities. Moreover, the modular design of PMoE allows experts of fundamentally different forms to coexist within a single framework, offering the flexibility to integrate turbulence models based on different formulations and assumptions as needed. Importantly, the training procedure relies on unsupervised feature extraction through the autoencoder-based router, making the framework particularly well suited to industrial workflows where labelled data and repeated full retraining are often impractical. Beyond its demonstrated performance, PMoE provides a scalable and extensible platform for RANS modelling. Owing to the sparse activation mechanism, model expansion does not lead to a proportional increase in computational cost during inference, which is essential for industrial CFD applications. Its plug-and-play structure enables modular updates and reduces both development and deployment costs, while its extensibility allows the turbulence modelling community to contribute new experts targeting previously unexplored flow regimes, establishing a living framework that can evolve alongside advances in modelling strategies.
Building on the modular and extensible nature of the framework, future developments may incorporate experts based on alternative turbulence closures, such as the
$k$
–
$\omega$
shear stress transport model (Menter Reference Menter1994) or Reynolds stress models (Launder, Reece & Rodi Reference Launder, Reece and Rodi1975), to broaden its applicability. For more complex flows, region-aware routing is expected to become increasingly important for applying different experts within different spatial zones. Although the current set of input features is sufficient for distinguishing the flow regimes examined here, further enrichment of the feature representation may be required as additional flow types are considered.
Funding
This work has been supported by the National Natural Science Foundation of China (grant nos 12588201, 12572247 and 12432010).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Feature importance analysis
Feature importance analysis is commonly employed to improve the interpretability of neural networks by quantifying the contribution of each input feature to the model prediction (Breiman Reference Breiman2001). In this work, to explore whether the importance distribution of different features varies in various flow regimes, and thereby demonstrate the feasibility of identifying flow types based on feature importance, it is necessary to conduct a feature importance analysis.
Among such approaches, PFI provides a model-agnostic assessment of the sensitivity of a fitted model to individual features using a prescribed tabular dataset (Fisher, Rudin & Dominici Reference Fisher, Rudin and Dominici2019; Mandler & Weigand Reference Mandler and Weigand2023). It is achieved by shuffling values of a single feature and subsequently assessing the consequent deterioration in the model performance. This approach can determine the extent to which the model depends on a specific feature by disrupting the intrinsic connection between the input features and the model prediction. The procedure of the PFI method is illustrated in figure 12. The inputs include a neural network used for reconstruction, an input feature matrix
${\unicode{x1D653}}\in \mathbb{R}^{N\times M}$
comprising
$M$
features and
$N$
sampling points, and the corresponding targets, which coincide with
$\unicode{x1D653}$
for this reconstruction task. For each feature
$q_i$
, a permuted matrix
${\unicode{x1D653}}_i^S$
is generated by shuffling the
$i$
th column of
$\unicode{x1D653}$
, and the reconstruction error
$e_i^S$
is recomputed as
$\|{\unicode{x1D653}}_{Ri}^S-{\unicode{x1D653}}\|$
. The importance associated with each feature
$q_i$
is defined as
$\varPsi _i$
, the change in error induced by permutation, providing a quantitative measure of the model’s sensitivity to that feature. All features are analysed in this manner, and the ranking is obtained from repeated trials with averaged results to minimise statistical variability.
Schematic illustration of the formulation and computational procedure of the PFI.
Figure 12. Long description
The diagram shows a matrix X with elements q, where columns are shuffled to create a new matrix X subscript S. This matrix is input into a neural network for reconstruction, producing output features X subscript R. The feature importance score is calculated using the differences between the original, shuffled, and reconstructed matrices. Arrows indicate the flow from input features through the network to output features and the calculation of error metrics.
To eliminate the interference caused by network capacity, we choose networks with sufficient parameters to test whether the four types of flowing data that we adopt have different orders of feature importance in the reconstruction task. To facilitate the comparison of relative feature importance both within and across heterogeneous datasets, the raw PFI scores were normalised by the sum of seven feature importance scores. In figure 13, the bar plots display the normalised mean PFI values, with error bars representing the standard deviation computed from
$20$
repeated permutation trials, indicating the stability of the importance estimates. Features with higher importance scores are considered more informative. The four flow regimes exhibit distinct importance profiles, indicating that different subsets of features govern the reconstruction. A more detailed inspection of the dominant feature combinations reveals that each canonical flow type admits a distinct low-dimensional physical representation in the SA-generated solution space. In the following, the four baseline flow categories are discussed separately. For the wake flows, the reconstruction is mainly governed by
$q_2$
,
$q_1$
and
$q_3$
, indicating that the flow is encoded by shear–rotation dynamics and their interaction with the pressure gradient, consistent with a free shear-layer structure. For the channel flow, the dominant features are
$q_5$
,
$q_4$
and
$q_7$
, corresponding to eddy viscosity, wall distance and shear anisotropy. This indicates that the flow is fully characterised by wall-bounded turbulence scaling, with negligible influence from rotational or pressure-driven effects. For the periodic hill flow, the importance is distributed across
$q_5$
,
$q_2$
,
$q_3$
and
$q_7$
, reflecting the coupled roles of turbulent transport, shear-layer instability, pressure gradient and recirculation. This confirms that separated flows require a multi-physics representation. For the square duct flow, the dominant features remain
$q_5$
,
$q_4$
and
$q_7$
, indicating a wall-driven anisotropic shear system. The weak contribution of
$q_2$
and
$q_3$
suggests that secondary-flow-related rotational and pressure effects are not explicitly represented in the SA solution. These class-dependent importance patterns lead to substantially different parameter distributions when training the modular experts of the router network. This explains why the approach of using reconstruction errors as the identification basis can effectively distinguish different flow types of data, and still maintain a high accuracy rate after modular expansion.
Normalised feature importance scores obtained from the PFI for the the training data and additional test cases: (a) wake, (b) channel, (c) periodic hill, (d) square duct, (e) 2DWMH, (f) CBFS, (g) RD2500.
Figure 13. Long description
The bar graph compares normalized feature importance scores across different test cases, including wake, channel, periodic hill, square duct, 2DWMH, CBFS, and RD2500. The x-axis represents different features labeled as q1, q2, q3, q4, q5, q6, and q7. The y-axis represents the normalized feature importance scores ranging from 0 to 0.6. Each subplot (a) to (g) corresponds to a specific test case: (a) wake, (b) channel, (c) periodic hill, (d) square duct, (e) 2DWMH, (f) CBFS, and (g) RD2500. The bars are color-coded to represent different features, with each color corresponding to a specific feature. The graph shows variations in feature importance across different test cases, highlighting which features are most significant in each scenario. All values are approximated.
For the CBFS case, the dominant contributions from
$q_2$
and
$q_3$
, together with secondary roles of
$q_1$
and
$q_7$
, place it within the same feature subspace as the periodic hill flow. Accordingly, the router identifies the periodic hill class
$C_2$
with confidence
$76.7\,\%$
, reflecting a strong overlap in the pressure–shear-dominated representation of separated flows. For the 2DWMH case, the PFI distribution remains dominated by
$q_2$
and
$q_3$
, accompanied by non-negligible contributions from
$q_1$
and
$q_5$
. This indicates that while the flow retains a separation mechanism governed by pressure–shear coupling similar to the periodic hill, the role of turbulent transport becomes more pronounced. At the same time, the presence of an upstream attached boundary layer introduces additional dependence on wall-related structures, leading to a more distributed feature representation. As a result, the router confidence for the periodic hill class
$C_2$
decreases to
$54.7\,\%$
, with a secondary attribution of
$35.2\,\%$
to the channel class
$C_1$
, reflecting the mixed character of external separation and partial wall-bounded behaviour. For the rectangular duct (RD2500), the dominant features remain
$q_5$
,
$q_4$
and
$q_7$
, closely matching the square duct flow. This consistency in feature space explains why the router assigns the highest confidence (
$90.8\,\%$
) to the square duct class
$C_3$
, indicating that the SA solution preserves the same wall-driven anisotropic shear representation.
Appendix B. Case settings
Numerical simulations in this paper are performed using the open-source computational framework OpenFOAM (open source field operation and manipulation, version 2012 ). The SA baseline model mentioned in this paper refers to the model and default parameter settings presented by Spalart & Allmaras (Reference Spalart and Allmaras1992). Grids and computation domain are shown in table 9. The numerical set-ups exhibit minor variations across cases due to differences in geometry and flow conditions, but follow consistent practices based on standard OpenFOAM configurations for canonical flows. A representative set-up employs second-order spatial discretisation, with Gauss linear schemes for gradients and diffusion terms, and a linear upwind scheme for the convective term of velocity. Turbulence transport equations are discretised using bounded upwind schemes. Pressure–velocity coupling is handled using the SIMPLE algorithm with one non-orthogonal corrector. Standard under-relaxation factors are applied. Convergence is assessed based on residual reduction and stabilisation of key flow quantities.
Computational domains and grid resolutions for all cases considered in this study.
Table 9. Long description
The table presents computational domains and grid resolutions for different cases in a study. It includes columns for Case, Box, and Grid, with specific values for each case. The Box column lists dimensions in various units, while the Grid column provides grid resolutions in pixel dimensions. The table contains 15 rows, each representing a different case with unique values for Box and Grid.
Appendix C. Training details
In this appendix, we provide the detailed training procedures for the three experts
$E_1$
,
$E_2$
and
$E_3$
, which are used in the PMoE framework. These procedures include the generation of training data, the optimisation objectives, and the implementation details of the symbolic regression and neural network models. The description is organised sequentially from
$E_1$
to
$E_3$
, highlighting the specific methodologies adopted for each expert.
C.1. Expert
$E_1$
: symbolic regression of the damping function
The damping function correction
$f_{\nu 1}$
for
$E_1$
is obtained via offline symbolic regression following the methodology of Weatheritt & Sandberg (Reference Weatheritt and Sandberg2016). Candidate expressions are generated and evolved using gene expression programming (GEP) under predefined function and terminal sets. The function set includes the basic operators
$\{+, -, \times , \div , \text{exponentiation}, \text{square}, \tanh \}$
, while the terminal set comprises the variable
$\chi = \tilde {\nu }/\nu$
, the constant
$1.0$
, and randomly generated constants sampled from the interval
$[-5, 5)$
. The GEP algorithm is initialised with a population of
$1000$
expressions, and evolved over
$10\,000$
generations through mutation and recombination.
Each candidate expression is evaluated on the training data using the objective
\begin{equation} \mathcal{J} = \frac {1}{N} \sum _{i=1}^{N} \left ( f_{\nu 1}^{\textit{pred}}(\chi _i) - f_{\nu 1}^{\textit{target}}(\chi _i) \right )^2, \end{equation}
where
$f_{\nu 1}^{\textit{pred}}$
denotes the symbolic prediction, and
$f_{\nu 1}^{\textit{target}}$
is the corresponding reference value. The training data are extracted from the corrected
$f_{\nu 1}(\chi )$
profiles reported by Bin et al. (Reference Bin, Huang and Yang2023), which were originally generated using a neural-network-based approach. These profiles capture the near-wall behaviour of the turbulence eddy viscosity, which is critical for accurately predicting the mean velocity distribution in wall-bounded flows. At the end of this offline regression stage, the optimal expression
$f_{\nu 1}(\chi )$
is selected according to the defined objective, and is subsequently hard-coded into the SA solver as
$E_1$
for RANS simulations, improving near-wall predictions without additional tuning.
C.2. Expert
$E_2$
: learning the correction of production term for flow separation
Expert
$E_2$
is trained to capture the non-equilibrium physics in the recirculation bubble using a two-step FIML procedure. First, an adjoint-based field inversion is performed on the periodic hill case PH1p0 using the DAFoam framework (He et al. Reference He, Mader, Martins and Maki2018, Reference He, Mader, Martins and Maki2020). The optimisation objective is to minimise the discrepancy between the RANS-predicted velocity field and DNS reference data (Xiao et al. Reference Xiao, Wu, Laizet and Duan2020) by adjusting a spatially varying production-term multiplier
$\beta (\boldsymbol{x})$
, yielding an optimal
$\beta$
field on the computational grid.
Second, a neural network is trained to learn the mapping from the input features
$\boldsymbol{x} = [q_1, \ldots , q_7]$
to the inverted
$\beta$
field. The dataset comprises grid-point-level pairs
$\{(\boldsymbol{x}_j, \beta _j)\}_{j=1}^{N_{\textit{grid}}}$
extracted from the PH1p0 case, and is randomly divided into
$80\,\%$
training and
$20\,\%$
testing subsets. Note that no extra sampling is applied at this stage. The network is implemented as a three-layer fully connected feed-forward network in PyTorch (Paszke et al. Reference Paszke2019). All input features are further normalised to the range
$ [0, 1 ]$
using min–max scaling, which is a standard pre-processing step to improve numerical conditioning during neural network optimisation. The network has seven input nodes, two hidden layers with
$64$
and
$32$
neurons, and a single scalar output. Rectified linear unit (Glorot, Bordes & Bengio Reference Glorot, Bordes and Bengio2011) activation functions are applied between layers to introduce nonlinearity. The network is trained using the Adam optimiser with a cosine-annealed learning-rate schedule and an adaptive reduction strategy triggered by validation-loss plateaus to enhance convergence robustness. The mean squared error loss function is used during training. This compact architecture provides sufficient representational power for nonlinear feature interactions while limiting the risk of overfitting. The trained network subsequently injects the predicted
$\beta (\boldsymbol{x})$
into the production term of the SA model, as defined in (2.12), for subsequent RANS simulations.
C.3. Expert
$E_3$
: learning the constitutive relation modification for secondary flows
Expert
$E_3$
is trained following the same two-step FIML procedure as
$E_2$
, but applied to the square duct cases SD2500 and SD5693. First, an adjoint-based field inversion optimises the spatially varying nonlinear stress coefficient
$C_{\textit{cr1}}(\boldsymbol{x})$
in the QCR formulation (2.13), with the objective of minimising the discrepancy between the RANS-predicted velocity field – including both streamwise and secondary-flow components – and DNS reference data (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010; Vinuesa et al. Reference Vinuesa, Schlatter and Nagib2018).
Second, a neural network with the same architecture and training procedure as
$E_2$
is employed to learn the mapping
$\boldsymbol{x} \to C_{\textit{cr1}}(\boldsymbol{x})$
. For brevity, we do not repeat the details of the network architecture and training procedure, which are identical to those described for
$E_2$
above. The trained network subsequently provides the
$C_{\textit{cr1}}(\boldsymbol{x})$
field for injection into the QCR-based RANS simulations.
Overall, the training of these three experts combines data-driven model extraction and physics-informed optimisation, enabling the PMoE framework to accurately capture both near-wall and non-equilibrium turbulent flow features across a variety of canonical cases.
𝙒
𝙎
ρ
ν
ν~
d
Lref
Ui
qi
q^i
qi∗
qi=q^i/(|q^i|+|qi∗|)
Ttot=Tbase+Troute+Texp
x
y
z=B/2
Cf
z=B/2
L
