Hybrid Physics-Informed Neural Networks Integrating Multi-Relaxation-Time Lattice Boltzmann Method for Forward and Inverse Flow Problems

1. Introduction

Developing stable, accurate, and generalizable numerical simulation methods has long been one of the core challenges in computational fluid dynamics (CFD) [1,2,3]. Traditional numerical methods, such as the finite difference method (FDM), finite volume method (FVM), and spectral method (SM), have established the basic framework of modern CFD through discretization theory, mesh partitioning, and algebraic equation solving [1,4]. These methods have become benchmark tools for simulating physical phenomena such as fluid dynamics and heat transfer. However, when dealing with strong nonlinearity, discontinuous solutions, and complicated geometries, these traditional methods face significant limitations, including heavy reliance on refined meshes, numerical instability, and severe accuracy degradation or even failure [5,6,7].

In recent years, physics-informed neural networks (PINNs) have emerged as a novel computational method that combines the advantages of physical laws and data-driven learning [8]. Leveraging the mesh-free computation, PINNs demonstrate great potential in addressing diverse fluid problems, gradually becoming a key new paradigm to break through the bottlenecks of traditional numerical methods [9]. Specifically, PINNs embed physical laws directly into the loss function of neural networks and utilize parameter optimization to implicitly solve the governing equations, thus largely eliminating mesh dependence and enhancing adaptability to complex physical problems [10,11,12]. For instance, this approach has demonstrated excellent performance in shock capturing, turbulence modeling, and unsteady flow simulations [13]. Nonetheless, PINNs also suffer from challenges including the emergence of spurious solutions, convergence difficulties, and limited capability for dynamic system modeling [10,14,15]. As a result, introducing hybrid model architectures to improve robustness has become a promising research direction. Among these approaches, the integration of PINNs with the lattice Boltzmann method (LBM) represents a recently developed research direction.

The LBM, owing to its clear evolution mechanism, highly parallel computational structure, and flexibility in handling intricate boundaries and multi-physics coupling problems, has been widely adopted for numerical simulation in CFD [16]. Its fundamental governing equation, the lattice Boltzmann equation (LBE), originates from the underlying kinetic equation, i.e., the continuous Boltzmann equation. This intrinsic physical nature rooted in kinetic theory enables the LBE to systematically preserve the essential coupling between microscopic particle dynamics and macroscopic conservation law evolution in physical modeling, thus endowing it with inherent physical consistency, theoretical closure, and high generality for fluid flow solutions [17,18]. The LBE models the mesoscopic evolution of particle distributions through collision and streaming processes, and offers a solid theoretical basis for complex flow simulations.

Based on this, some researchers have proposed combining LBM with PINNs, replacing conventional partial differential equations, such as Navier-Stokes (N-S) equations, with the LBE as the underlying governing equations, thereby constructing hybrid models that integrate mesoscopic physical constraints [19]. For example, Lou et al. proposed a PINN architecture incorporating the Boltzmann-Bhatnagar-Gross-Krook (BGK) equation, substantially enhancing the model physical accuracy, and alleviating issues such as slow convergence and difficulties in modeling highly nonlinear scenarios [20]. However, due to coupling the relaxation of all hydrodynamic and non-hydrodynamic moments to a single parameter, the BGK model exhibits severe parameter coupling and numerical instability in flow simulations, thus limiting its applicability in complex scenarios [16,17].

To address these challenges, this paper proposes a PINN architecture incorporating the multi-relaxation-time lattice Boltzmann method (MRT-LBM). The MRT collision mechanism effectively decouples the relaxation processes of different order momentum components by assigning independent relaxation time to each physical mode [21,22]. This approach retains the advantages of mesoscale modeling while markedly improving the accuracy, numerical stability, and generalizability for diverse fluid structures. The remainder of this paper is organized as follows: Section 2 presents the theoretical foundation, along with the design of the proposed PINN-MRT architecture that integrates dual-network and physics-informed loss functions. Section 3 introduces the benchmark problem setup and model configurations. Section 4 demonstrates the applications of PINN-MRT to forward and inverse problems, including predictive accuracy and parameter identification at different Reynolds numbers (Re). Finally, Section 5 summarizes the main findings and outlines future research directions.

2. Method

2.1. MRT-LBM

The LBM is an efficient numerical simulation approach at the mesoscopic scale, describing the statistical evolution of particle distribution functions through discretization in both space and time, as well as along specific discrete velocity directions. The fundamental dynamical model of LBM originates from the continuous Boltzmann equation [23]

$${\displaystyle \frac{\partial f}{\partial t}}+\xi ·\nabla f=\Omega \left(f\right),$$

(1)

where f denotes the particle distribution function in continuous velocity space, $\xi$ represents the microscopic particle velocity, and $\Omega \left(f\right)$ is the collision operator, characterizing the relaxation process toward local equilibrium. A classical simplification is the single-relaxation-time (SRT) or BGK approximation, which linearizes the collision operator to model the relaxation process toward a local equilibrium distribution function [24]

$${\displaystyle \frac{\partial f}{\partial t}}+\xi ·\nabla f=-{\displaystyle \frac{1}{\tau }}(f-f^{eq}),$$

(2)

$$f^{eq}={\displaystyle \frac{\rho }{{\left(2\pi RT\right)}^{D/2}}}exp\left[-{\displaystyle \frac{{|\xi -\mathbf{u}|}^2}{2RT}}\right].$$

(3)

In Eq.2, $\tau$ denotes the single relaxation factor, which directly determines the fluid kinematic viscosity, and $f^{eq}$ represents the equilibrium distribution function. However, the SRT model employs a unified relaxation parameter across all dynamical moments, which inherently leads to strong coupling among different physical processes and restricts the flexibility in independently controlling individual transport phenomena [24]. Fortunately, the MRT model introduces a transformation into moment space assigning independent relaxation parameters to distinct physical modes (mass, momentum, stress, etc.). The corresponding equation is shown as [25,26]

$${\displaystyle \frac{\partial f}{\partial t}}+\xi ·\nabla f=-\left({\mathbf{M}}^{-1}\mathrm{\Lambda }\mathbf{M}\right)(f-f^{eq}),$$

(4)

where $\mathbf{M}$ denotes an orthogonal transformation matrix that maps f from velocity space to moment space, and $\mathrm{\Lambda }$ is a diagonal matrix whose elements correspond to the relaxation times associated with individual moments. The D2Q9 discrete velocity model is adopted, and the resulting discretized evolution equation is given by

$${\displaystyle \frac{\partial f_i}{\partial t}}+{\mathbf{e}}_i·\nabla f_i=-{\left({\mathbf{M}}^{-1}\mathrm{\Lambda }\mathbf{M}\right)}_{ij}\left(f_j-f_j^{eq}\right),$$

(5)

where $f_i$ is the particle distribution function. $f_j^{eq}$ is the equilibrium discrete distribution function, and for the D2Q9 discrete velocity model can be expressed as

$$f_j^{eq}(\mathbf{x},t)=w_i\rho \left[1+{\displaystyle \frac{{\mathbf{e}}_i·\mathbf{u}}{c_s^2}}+{\displaystyle \frac{{({\mathbf{e}}_i·\mathbf{u})}^2}{2c_s^4}}-{\displaystyle \frac{\mathbf{u}·\mathbf{u}}{2c_s^2}}\right],$$

(6)

where $w_i$ is weight coefficients in the ith direction, which is given by

$$w_i=\left\{\begin{array}{cc}{\displaystyle \frac{4}{9}},\hfill & i=0\hfill \\ {\displaystyle \frac{1}{9}},\hfill & i=1,2,3,4\hfill \\ {\displaystyle \frac{1}{36}}.\hfill & i=5,6,7,8\hfill \end{array}\right.$$

(7)

The macroscopic fluid density $\rho$ and macroscopic fluid velocity $\mathbf{u}$ are shown as [27]

$$\begin{array}{ccc}\hfill & \rho =\sum _i{f}_i,\hfill & \hfill \rho \mathbf{u}=\sum _i{\mathbf{e}}_i{f}_i,\end{array}$$

(8)

where ${\mathbf{e}}_i$ discrete velocity along the ith direction. In this study, ${\mathbf{e}}_i$ is given by

$${\mathbf{e}}_i=\left\{\begin{array}{cc}(0,0),\hfill & i=0\hfill \\ \left(cos\left[(i-1){\displaystyle \frac{\pi }{2}}\right],sin\left[(i-1){\displaystyle \frac{\pi }{2}}\right]\right),\hfill & i=1,2,3,4\hfill \\ \sqrt{2}\left(cos\left[(2i-1){\displaystyle \frac{\pi }{4}}\right],sin\left[(2i-1){\displaystyle \frac{\pi }{4}}\right]\right).\hfill & i=5,6,7,8\hfill \end{array}\right.$$

(9)

The orthogonal transformation matrix is expressed as [27]

$$\mathbf{M}=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2\\ 4& -2& -2& -2& -2& 1& 1& 1\\ 0& 1& 0& -1& 0& 1& -1& -1\\ 0& -2& 0& 2& 0& 1& -1& -1\\ 0& 0& 1& 0& -1& 1& 1& -1\\ 0& 0& -2& 0& 2& 1& 1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1\end{array}\right].$$

(10)

The diagonal relaxation matrix $\mathrm{\Lambda }$ is shown as

$$\mathrm{\Lambda }=\mathrm{diag}\left({\tau }_{\rho }^{-1},{\tau }_e^{-1},{\tau }_f^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_{\nu }^{-1}\right),$$

(11)

where ${\tau }_{\rho }^{-1}$, ${\tau }_e^{-1}$, ${\tau }_f^{-1}$, ${\tau }_j^{-1}$, ${\tau }_q^{-1}$, and ${\tau }_{\nu }^{-1}$ represent the relaxation factors associated with density, energy, square root of energy, momentum flux component, energy flux component, and pressure tensor, respectively. In this study, the relaxation factors are set as follows, ${\tau }_{\rho }^{-1}={\tau }_j^{-1}=1.0$, ${\tau }_e^{-1}={\tau }_f^{-1}=0.8$, and ${\tau }_q^{-1}=1.1$. In the MRT model, the kinematic viscosity $\nu$ is related to the relaxation factor ${\tau }_{\nu }$ by

$$\nu =c_s^2\left({\tau }_{\nu }-{\displaystyle \frac{1}{2}}\right)\delta t,$$

(12)

where $c_s$ denotes the lattice sound speed, with $c_s=\sqrt{1/3}$ in the D2Q9 model. $\delta t$ is the time step.

2.2. PINN-MRT

A novel architecture (PINN-MRT) is developed in this study by integrating MRT-LBM into the PINNs architecture. As shown in Figure 1, the PINN-MRT architecture comprises two main components, i.e. the neural network architecture and the loss function design.

2.2.1. Neural Network Architecture

The neural network architecture has two sub-networks. The first sub-network, denoted as $N{N}_{eq}$, takes the spatial coordinates $(x,t)$ as input and outputs the macroscopic density $\rho$ and velocity components $(u,v)$, then calculates the equilibrium distribution function $f_i^{eq}$ using Eq.6. The second sub-network, $N{N}_{neq}$, shares the same input as $N{N}_{eq}$ but predicts the nine components of the non-equilibrium distribution function $f_i^{neq}$. This physically motivated separation leverages the smoothness of $f_i^{eq}$ and the complexity of $f_i^{neq}$, allowing each network to focus on distinct features and thereby improving stability and accuracy [28].

Based on the predicted outputs, the total distribution function at each spatial location is composed as

$$f_i=f_i^{eq}+f_i^{neq}.$$

(13)

2.2.2. Loss Function Design

For effective training of the PINN-MRT model, a loss function is formulated incorporating the physical residuals, boundary conditions, initial conditions, and data-driven terms, as shown in Eq.14 [29].

$$\mathcal{L}={\beta }_1{L}_{\mathrm{PDE}}+{\beta }_2{L}_{\mathrm{BC}}+{\beta }_3{L}_{\mathrm{IC}}+{\beta }_4{L}_{\mathrm{DATA}},$$

(14)

where ${\beta }_i(i=1,2,3,4)$ denote the weight coefficients corresponding to the respective loss components. To ensure consistency in comparisons with other models, this study does not implement weight optimization and instead adopts unit weight coefficients (i.e., ${\beta }_i=1$) for active loss terms, while setting ${\beta }_i=0$ for unused components.

The physics residual loss $L_{\mathrm{PDE}}$ ensures that the neural network predictions satisfy the governing dynamics derived from the MRT-LBM equations. The specific equation can be expressed as [30]

$$L_{\mathrm{PDE}}={\displaystyle \frac{1}{N_e}}\sum _{i=0}^{Q-1}\sum _{j=0}^{N_t-1}{\left|R_i({\mathbf{x}}_j,t_j)\right|}^2,$$

(15)

$$R_i={\displaystyle \frac{\partial f_i}{\partial t}}+{\mathbf{e}}_i·\nabla f_i+{\left({\mathbf{M}}^{-1}\mathrm{\Lambda }\mathbf{M}\right)}_{ij}\left(f_j-f_j^{eq}\right),$$

(16)

where $R_i$ denotes the residual of the MRT-LBM evolution in the i-th discrete velocity direction, Q refers to the total number of discrete velocity directions, and $N_t$ and $N_e$ specify the numbers of temporal and internal sampling points, respectively. The spatial coordinate for the j-th sampling point is given by ${\mathbf{x}}_j$, with $t_j$ representing its corresponding temporal coordinate.

The boundary condition loss $L_{\mathrm{BC}}$ is introduced to enforce that the predicted solutions strictly satisfy the prescribed boundary conditions, thereby ensuring physical consistency at the domain boundaries. In the LBM kinetic framework, it is further divided into three components

• Macroscopic boundary loss $L_{m_{BC}}$: Constrains the predicted macroscopic variables on the boundaries to match the prescribed physical values.
• Distribution function consistency loss $L_{f_{BC}}$: Constrains the predicted equilibrium and non-equilibrium distribution functions to match the theoretically defined boundary distributions.
• Boundary PDE residual loss $L_{e_{BC}}$: Imposes constraints on the neural network outputs by penalizing the residuals of the governing equations at the boundaries.

The equations for the three components are given by

$$L_{m_{BC}}={\displaystyle \frac{1}{N_b}}\sum _{j=0}^{N_b-1}{\left|{\varphi }_{N{N}_{eq}}({\mathbf{x}}_{b,j},t_j)-{\varphi }^{*}({\mathbf{x}}_{b,j},t_j)\right|}^2,$$

(17)

$$L_{f_{BC}}={\displaystyle \frac{1}{N_b}}\sum _{j=0}^{N_b-1}\sum _{{\xi }_i·\mathbf{n}>0}{\left|f_{i,N{N}_{neq}}({\mathbf{x}}_{b,j},t_j)-f_{i,neq}^{*}({\mathbf{x}}_{b,j},t_j)\right|}^2,$$

(18)

$$L_{e_{BC}}={\displaystyle \frac{1}{N_b}}\sum _{j=0}^{N_b-1}\sum _{{\xi }_i·\mathbf{n}<0}{\left|R_i({\mathbf{x}}_{b,j},t_j)\right|}^2.$$

(19)

Thus, the formula for $L_{\mathrm{BC}}$ is

$$L_{\mathrm{BC}}=L_{m_{BC}}+L_{f_{BC}}+L_{e_{BC}}.$$

(20)

In Eq.17, where ${\varphi }_{N{N}_{eq}}$ denotes the predicted macroscopic physical quantities (e.g., velocity or density), ${\varphi }^{*}$ represents the ground truth boundary condition, $N_b$ is the number of boundary sampling points, and ${\mathbf{x}}_{b,j}$ and $t_j$ denote the spatial and temporal coordinates, respectively, of the j-th boundary point. In Eq.18, $f_{i,N{N}_{neq}}^{neq}$ represents the non-equilibrium distribution function in the i-th discrete velocity direction predicted by the neural network, $f_i^{neq,*}$ is the theoretical reference value at the boundary, and ${\xi }_i·\mathbf{n}>0$ indicates that only outgoing velocity directions are considered (from the interior to the boundary). Eq.19 indicates that the evolution equation in the incident direction remains enforced at the boundary.

The initial condition loss $L_{\mathrm{IC}}$ ensures that the predicted initial state strictly conforms to the known reference solution. The equation is given by

$$L_{\mathrm{IC}}={\displaystyle \frac{1}{N_i}}\sum _{j=0}^{N_i-1}{\left|{\varphi }_{N{N}_{eq}}({\mathbf{x}}_j,0)-{\varphi }^{*}({\mathbf{x}}_j,0)\right|}^2,$$

(21)

with $N_i$ representing the number of sampling points.

The data-driven loss $L_{\mathrm{DATA}}$ minimizes the discrepancy between the predicted values and the experimentally measured or observed data. It is defined as

$$L_{\mathrm{DATA}}={\displaystyle \frac{1}{N_d}}\sum _{j=0}^{N_d-1}{\left|u({\mathbf{x}}_{u_j},t_{u_j})-u^{*}({\mathbf{x}}_{u_j},t_{u_j})\right|}^2,$$

(22)

where $N_d$ denotes the total number of data points, $u^{*}$ represents the observed velocity data, thereby guiding the neural network to better approximate the actual physical data.

3. Benchmark Modeling

The classical lid-driven cavity flow is selected as the benchmark case to evaluate the model performance. As shown in Figure 2, the computational domain is a unit square with spatial coordinates $(x,y)\in {[0,1]}^2$. A constant velocity of $U=1$, $V=0$ is prescribed on the top boundary, and the remaining walls are subjected to no-slip conditions (i.e., stationary). This setup induces a characteristic steady-state vortex structure that varies with Re, making it a widely adopted benchmark for assessing the predictive capability of both forward and inverse modeling approaches. The specific governing equations are written as [31]

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=0,\hfill \\ \hfill & u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+\nu \left({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right),\hfill \\ \hfill & u{\displaystyle \frac{\partial v}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial y}}+\nu \left({\displaystyle \frac{{\partial }^{2}v}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}v}{\partial y^2}}\right),\hfill \end{array}$$

(23)

where $(u,v)$ denote the velocity components in the x and y directions, respectively, p is the pressure, and $\nu$ represents the kinematic viscosity, which can be expressed in terms of the Re as

$$\nu =UL/Re,$$

(24)

here, U represents the lid velocity, and L denotes the characteristic length (i.e., the cavity length).

These governing equations are formulated from the general N-S equations under steady-state, incompressible, and two-dimensional assumptions.

To rigorously evaluate the performance of the proposed PINN-MRT architecture, two benchmark models are established for comparative analysis. All models share the same network architecture and training configuration, differing solely in the formulation of the residual term $L_{\mathrm{PDE}}$,which encodes the respective physical constraints [32]. The standard PINN adopts Eq.23 as its macroscopic residual constraints. Meanwhile, PINN-SRT introduces mesoscopic physical constraints based on SRT equation, which can be regarded as a special case of the MRT model when all the diagonal elements of the relaxation matrix are set to unity [32].

The detailed configurations for each model in inverse problem are listed in Table 1. The training data are obtained by solving Eq.23 using the traditional FDM, which also provides the reference numerical solutions for comparison in Section.Section 4.

Table 2 summarizes the corresponding configurations for forward problem. In this context, all physical parameters are known, and no additional data are introduced. A deeper network architecture (7-layers) is adopted to handle this more challenging data-free prediction task. Meanwhile, the sub-network for the non-equilibrium part uses more neurons (100 vs. 60) to better capture complex local features [10,14,33,34].

4. Results

4.1. Inverse Problem

The subsequent analysis of the inverse problem is divided into two parts: Flow Field Analysis, focusing on the reconstruction of physical fields, and Parameter Inversion, addressing the identification of underlying model parameters.

4.1.1. Flow Field Analysis

The comprehensive performance comparison among the standard PINN, PINN-SRT and PINN-MRT models is conducted up to Re = 1000, since the standard PINN and PINN-SRT models become numerically unstable and fail to produce reliable results when the Re exceeds 1000. In this range, all models are systematically analyzed in terms of vortex structure reconstruction, error distribution, and velocity profile consistency. To further verify the robustness of the PINN-MRT model under more challenging flow conditions, additional evaluations at Re = 2000 and 5000 are conducted exclusively for the PINN-MRT model.

(A) Vortex Structure Reconstruction and Error Distribution

Figure 3 illustrates the predicted u and v velocity fields at Re = 100, showing small differences among three models. The absolute error (Abs) distributions in Figure 4 further confirm this observation. While the errors of all three models remain comparably low, PINN-MRT achieves the smallest and most localized errors. These results motivate further exploration at higher Re.

Figure 5 presents the flow field results at Re = 1000, where more complex structures appear compared to low Re. The primary vortex is mainly identified in the u velocity field as a large central rotational structure, whereas the secondary vortices are clearly observed in the v velocity field near the bottom corners. Among the models, the PINN-MRT accurately reconstructs both primary and secondary flow structures features, capturing their rotational directions and spatial distributions. By comparison, the PINN-SRT partially recovers the primary vortex but fails to resolve the secondary features. The standard PINN exhibits a complete breakdown in prediction. As shown in Figure 6, the PINN-MRT maintains low error levels, with peaks errors below 0.10 for u and 0.15 for v. Meanwhile, the PINN-SRT reaches peak errors of 0.30 with wider distributions, and the standard PINN yields errors exceeding 0.40, indicating clear failure.

(B) Velocity Profile

Centerline velocity profiles are analyzed to evaluate the accuracy of flow reconstruction along key directions. Table 3 summarizes the absolute error ranges of the horizontal and vertical velocity components for Re = 100 and 1000. The PINN-MRT model achieves the narrowest error bounds in both directions, confirming its superior predictive accuracy.

(C) High-Re Error Evaluation of PINN-MRT

The summarized relative $L^2$ errors of the PINN-MRT model across all tested Re are presented in Table 4. Here, ${\mathrm{err}}_{\mathrm{u}}$ and ${\mathrm{err}}_{\mathrm{u}}$ represent the horizontal and vertical velocity component errors, respectively, and ${\mathrm{err}}_{\mathrm{total}}$ denotes the total velocity magnitude error. As the Re increases from 100 to 5000, the errors naturally rise, with ${\mathrm{err}}_{\mathrm{u}}$ increasing from 9.067% to 15.281%, ${\mathrm{err}}_{\mathrm{v}}$ from 11.942% to 21.262%, and ${\mathrm{err}}_{\mathrm{total}}$ from 4.549% to 14.121%. Nevertheless, the error levels remain well controlled, demonstrating the robustness and practical applicability of the PINN-MRT model.

4.1.2. Parameter Inversion

Building on the successful flow field reconstruction, the proposed PINN-MRT model is further assessed for its capability in parameter inversion, focusing on both the kinematic viscosity coefficient $\nu$ and the lid velocity $U_{\mathrm{lid}}$. The inversion accuracy is systematically compared among the standard PINN, PINN-SRT, and PINN-MRT models across Reynolds numbers ranging from $\mathrm{Re}=100$ to 5000, as well as under noise-contaminated conditions.

As shown in Figure 7, the left panels present the relative $L^2$ errors of the inferred $\nu$ and $U_{\mathrm{lid}}$ for different models, while the right panels illustrate the corresponding robustness tests with noisy data, where the shaded regions denote a $\pm 5\%$ uncertainty band.

For the noiseless cases (left column), both PINN and PINN-SRT exhibit a rapid degradation of accuracy as Re increases, with relative errors exceeding $80\%$ and $50\%$, respectively, at $\mathrm{Re}=5000$. In contrast, the PINN-MRT maintains errors below $20\%$ for both parameters across the entire Reynolds number range, demonstrating remarkable stability and resistance to ill-conditioning at high Reynolds numbers.

Under noisy inputs (right column), the PINN-MRT results remain within the $\pm 5\%$ confidence interval, showing only marginal increases in error. This confirms that the MRT-based regularization effectively suppresses the amplification of noise during inverse inference, thereby ensuring both robustness and consistency of the reconstructed physical parameters.

4.2. Forward Problem

4.2.1. Flow Field Analysis

The forward problem focuses on evaluating the capability of each model to predict steady-state flow fields given known physical parameters. In this scenario, no additional data-driven constraints are introduced. Comparative evaluations are conducted at Re = 100 and 400 for all models (PINN, PINN-SRT, and PINN-MRT). For higher Re (Re > 400), only the PINN-MRT model is considered, as the other two models fail to produce reliable predictions under these more challenging conditions.

(A) Vortex Structure Reconstruction and Error Distribution

As shown in Figure 9, at Re = 100, the absolute error distributions demonstrate that the PINN-MRT model achieves the lowest error levels, highlighting its superior accuracy under low Re conditions.

At Re = 400 (Figure 11), the absolute error distributions clearly illustrate significant performance differences among the models. Compared to Re = 100, the errors of the PINN-MRT model show a slight increase. In contrast, both the standard PINN and PINN-SRT models have completely failed.

(B) Velocity Profile

As summarized in Table 5, at Re = 100, both the standard PINN and PINN-SRT models exhibit error ranges in $\Delta u$ and $\Delta v$ approximately twice as large as those of PINN-MRT. When Re increases to 400, the $\Delta u$ of PINN reaches nearly three times, and that of PINN-SRT about five times, compared to PINN-MRT. For $\Delta v$, the errors are even more pronounced, with PINN and PINN-SRT reaching approximately seven and six times larger than PINN-MRT, respectively.

(C) High Re Error Evaluation of PINN-MRT

As presented in Table 6, the PINN-MRT model performs robustly across all tested conditions. Notably, the increase in the vertical velocity component error is minimal, only 1.801% across the entire Re range. The horizontal velocity and total velocity component errors exhibit increases of approximately 5.000%. These restrained error increments underscore the model’s robust stability and accuracy, even under high Re conditions.

4.2.2. Viscosity Sensitivity

To further evaluate the stability of the model with respect to parameter perturbations, a sensitivity study was conducted by introducing controlled variations in the kinematic viscosity coefficient $\nu$. Figure 12 presents the relative $L^2$ errors of the reconstructed velocity fields (u, v, and total) under viscosity perturbations ranging from $0\%$ to $20\%$ for different Reynolds numbers. The shaded regions represent $\pm 3\%$ uncertainty bands.

Across all cases, the PINN-MRT model maintains stable error levels, with total relative errors remaining within the $\pm 3\%$ range, indicating low sensitivity to variations in $\nu$. The v-component shows relatively higher errors than u, particularly at higher Reynolds numbers. This behavior is attributed to the weak vertical velocity magnitude and the strong secondary vortex effects in the lid-driven cavity flow, which amplify relative deviations when minor physical perturbations occur. Overall, the results confirm that the proposed PINN-MRT framework achieves consistent flow reconstruction accuracy and robustness against viscosity perturbations in the forward modeling task.

5. Conclusions

In this study, the proposed PINN-MRT model is systematically applied to both forward and inverse problems, covering Re ranging from 100 to 5000. For inverse problems, it achieves high-fidelity parameter identification and flow field reconstruction, with remarkable robustness. For forward problems, the architecture exhibits outstanding predictive accuracy, precisely capturing intricate flow structures and velocity profiles. In addition, a sensitivity analysis on the kinematic viscosity $\nu$ demonstrates that the model maintains stable convergence under $\pm 3\%$ perturbations, confirming its insensitivity to small variations in physical parameters. These results highlight the robustness and predictive accuracy of the PINN-MRT framework across a broad range of flow conditions, including challenging high Re scenarios.

The primary contribution of this study is the introduction of the MRT collision mechanism into the PINN architecture. By assigning independent relaxation factors to different moment components, this approach effectively mitigates the critical trade-offs between stability, accuracy, and generalization that often constrains conventional PINN methods. This leads to a significant enhancement in the overall robustness and convergence of the training process.

It is worth noting that the current implementation employs basic fully connected networks, simple optimization strategies, and the fundamental MRT formulation.

In future work, this framework can be further improved by leveraging advanced neural architectures and more efficient training strategies to enhance convergence and scalability. For example, adaptive weighting mechanisms and dynamic optimization schemes will be introduced to better balance multiple loss components and improve numerical stability. Exploring advanced PINN variants, such as XPINNs [35], PDPINNs [36], and GPINNs [37], will also extend the applicability of the framework to high-dimensional and noisy flow systems. Furthermore, the PINN-MRT architecture can be extended through integrating specialized physics modules, such as pseudopotential models for multiphase and multicomponent flows, thermal formulations for heat transfer problems, and advanced boundary schemes for complex geometries. These extensions are expected to improve the capability of the framework to resolve complex interactions such as multiphase interfaces, thermal effects, and multiscale dynamics, while the pursuit of higher generalizability remains a key objective. In summary, the PINN-MRT model can be viewed as a bridge between LBM and deep learning, offering a promising foundation upon which more sophisticated and robust hybrid models can be built in future research.