2.1. Governing equations and boundary conditions

Under the given assumptions, the continuity and momentum equations for this problem are as follows (Nield & Bejan Reference Nield and Bejan2006):

(2.2) \begin{equation} \boldsymbol{\nabla } \cdot \boldsymbol{\textit {v}}=0 , \end{equation}

(2.3) \begin{equation} 0=-\boldsymbol{\nabla }p+\frac {\mu _{e}}{\epsilon }{\nabla} ^2 \boldsymbol{\textit {v}}-\mu \boldsymbol{\textit {K}}^{-1}\boldsymbol{\textit {v}}-\frac {1.75\rho }{\sqrt {150}\epsilon ^{3/2}} \boldsymbol{\textit {K}}^{-1/2} \boldsymbol{\textit {v}}|\boldsymbol{\textit {v}}|. \end{equation}

In (2.2)–(2.3), $\boldsymbol{\textit {v}}$ is the velocity vector, $\rho$ is the density, $p$ is the pressure and $\mu$ is the viscosity of the fluid. Additionally, $\mu _{e}$ is the effective viscosity of the medium and $\epsilon$ is the porosity of the medium.

In the current study, the convective acceleration term $\boldsymbol{\textit {v}}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\textit {v}}$ is neglected because the porous medium is assumed to possess very low permeabilities ( $K_1$ and $K_2$ ). Under such conditions, the flow resistance is predominantly governed by the Darcy drag term $\mu \boldsymbol{\textit {K}}^{-1}\boldsymbol{\textit {v}}$ together with the Forchheimer inertial correction term $ ({1.75\rho }/{\sqrt {150}\,\epsilon ^{3/2}}) \boldsymbol{\textit {K}}^{-1/2}\boldsymbol{\textit {v}}|\boldsymbol{\textit {v}}|.$ The strong resistance induced by the low permeability significantly suppresses the fluid velocity, thereby rendering the macroscopic convective acceleration term negligible compared with the Darcy and Forchheimer resistance terms. This assumption is consistent with theoretical and experimental studies showing that, in low-permeability porous media, the pressure drop is primarily dominated by Darcy and Forchheimer effects (Vafai & Tien Reference Vafai and Tien1981).

It is evident that $\det (\boldsymbol{\textit {K}}) = K_1 K_2$ is non-zero, implying that the inverse $\boldsymbol{\textit {K}}^{-1}$ exists. The Forchheimer term involves the matrix expression $\boldsymbol{\textit {K}}^{-1/2}$ . To determine its existence, it is observed that $\boldsymbol{\textit {K}}$ is a symmetric positive semi-definite matrix with positive eigenvalues $K_1$ and $K_2$ . According to Horn & Johnson (Reference Horn and Johnson2012), there exists a unique Hermitian positive semidefinite matrix $\boldsymbol{\textit {L}}$ such that $\boldsymbol{\textit {L}}^m = \boldsymbol{\textit {K}}$ for a given $m$ . Specifically, for $m=2$ , this guarantees the existence of $\boldsymbol{\textit {K}}^{-1/2}$ . It can be computed explicitly using diagonalisation.

The computational domain of the channel is given by $ \varOmega =[0,x_{fd}] \times [-1,1],$ where $x_{fd}$ denotes the hydrodynamic development length of the channel. The boundary of the domain is represented by $ \varGamma =\varGamma _{\textit{in}}\cup \varGamma _{w}\cup \varGamma _{\textit{exit}},$ where $\varGamma _{\textit{in}}$ , $\varGamma _{w}$ and $\varGamma _{\textit{exit}}$ correspond to the inlet boundary, channel walls and outlet (exit) boundary, respectively, given by

(2.4) \begin{align} \varGamma _{\textit{in}}=\{ (x,y)\in \varGamma :x=0 \}, \end{align}

(2.5) \begin{align} \varGamma _{w}=\{ (x,y)\in \varGamma :y=-L, y=L \}, \end{align}

(2.6) \begin{align} \varGamma _{\textit{exit}}=\{ (x,y)\in \varGamma :x=x_{fd} \}. \end{align}

At the inlet ( $\varGamma _{\textit{in}}$ ), a Dirichlet-type boundary condition is prescribed by imposing a uniform axial inflow velocity,

(2.7) \begin{equation} \boldsymbol{\textit {v}}|_{\varGamma _{\textit{in}}}=(u_0,0)^T. \end{equation}

Along the channel walls ( $\varGamma _w$ ), the classical no-slip condition is enforced,

(2.8) \begin{equation} \boldsymbol{\textit {v}}|_{\varGamma _{w}}=(0,0)^T, \end{equation}

which is physically consistent with stationary impermeable solid boundaries. Although the porous medium is anisotropic, the directional permeability modifies the internal momentum resistance within the porous structure and does not alter the kinematic restriction imposed by rigid walls. Therefore, the no-slip assumption remains valid for the present DBF model.

At the outlet ( $\varGamma _{\textit{exit}}$ ), the fully developed boundary condition, which is of Neumann-type, is imposed,

(2.9) \begin{equation} \boldsymbol{\textit {v}}|_{{\varGamma }_{\textit{exit}}}=\{ (u,0):\frac {\partial u}{\partial x}=0 \}, \end{equation}

which assumes that the axial velocity gradient vanishes in the streamwise direction at the downstream boundary.

2.2. Non-dimensionalisation

The following dimensionless variables are considered for this problem:

(2.10) \begin{equation} \mathit{\boldsymbol {X}}=\frac {\mathit{\boldsymbol {x}}}{L},\boldsymbol{{V}}=\frac {\boldsymbol{\textit {v}}}{u_{0}},P=\frac {p}{\rho (u_{0})^2}. \end{equation}

After non-dimensionalisation, (2.2)–(2.3) become

(2.11) \begin{equation} \boldsymbol{\nabla } \cdot \boldsymbol{{V}}=0, \end{equation}

(2.12) \begin{equation} 0 = -\boldsymbol{\nabla }P + \frac {M}{\epsilon \textit{Re}}{\nabla} ^2 \boldsymbol{{V}} - \frac {1}{\textit{ReDa}}\varLambda \boldsymbol{{V}} - \frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} \varSigma \boldsymbol{{V}}|\boldsymbol{{V}}|. \end{equation}

In (2.11)–(2.12), $\varLambda$ and $\varSigma$ are second-order tensors, given by

(2.13) \begin{equation} \varLambda = \left [ {\begin{array}{cc} \sin ^2\alpha +K\cos ^2\alpha & \frac {1}{2}(K-1)\sin 2\alpha \\[3pt] \frac {1}{2}(K-1)\sin 2\alpha & K\sin ^2\alpha +\cos ^2\alpha \\ \end{array} } \right ]\!, \end{equation}

(2.14) \begin{equation} \varSigma = \left [ {\begin{array}{cc} \sin ^2 \alpha + \sqrt {K} \cos ^2 \alpha & \frac {1}{2}(\sqrt {K}-1)\sin 2\alpha \\[3pt] \frac {1}{2}(\sqrt {K}-1)\sin 2\alpha & \sqrt {K} \sin ^2 \alpha + \cos ^2 \alpha \\ \end{array} } \right ]\!. \end{equation}

Also, in (2.11)–(2.12), $M$ is the viscosity ratio, $\textit{Re}$ is the Reynolds number, $\textit{Da}$ is the Darcy number and $K$ is the permeability ratio, given by

(2.15) \begin{align} \begin{split} M=\frac {\mu _{e}}{\mu }, \quad \textit{Re}=\frac {\rho u_0 L}{\mu },\quad \textit{Da}=\frac {K_1}{L^2}, \quad K=\frac {K_1}{K_2}. \end{split} \end{align}

Dimensionless forms of the respective boundary conditions are

(2.16) \begin{equation} \boldsymbol{{V}}|_{\varGamma _{\textit{in}}}=(1,0)^T, \end{equation}

(2.17) \begin{equation} \boldsymbol{{V}}|_{\varGamma _{w}}=(0,0)^T, \end{equation}

(2.18) \begin{equation} \boldsymbol{{V}}|_{{\varGamma }_{\textit{exit}}}=\{ (V,0)^T:\frac {\partial V}{\partial X}=0 \}. \end{equation}

2.3. Weak formulation

The weak formulation provides a variational representation of the governing DBF equations by relaxing the differentiability requirements imposed on the classical solution. Instead of enforcing the governing equations point-wise, the weak formulation requires them to hold in an integral sense over the computational domain. This approach is particularly suitable for nonlinear porous medium flows and forms the theoretical foundation for establishing existence and uniqueness of weak solutions.

Since the velocity field satisfies Dirichlet boundary conditions on the inlet and wall boundaries, the admissible velocity test space is defined as

(2.19) \begin{equation} W=\{\boldsymbol{{w}}\in H^1(\varOmega )\times H^1(\varOmega ):\boldsymbol{{w}}|_{\varGamma _{\textit{in}}\cup \varGamma _w}=0\}, \end{equation}

where $H^1(\varOmega )$ is the Sobolev space of square-integrable functions with square-integrable first-order weak derivatives. The pressure test space is chosen as

(2.20) \begin{equation} Q=L^2(\varOmega ), \end{equation}

which is the standard functional setting for incompressible flow problems. Since the pressure appears only through its first derivatives in the momentum equation, square integrability is sufficient. More details about the functional spaces are given in Appendix A.

The anisotropic permeability tensor $\boldsymbol{\textit {K}}$ is symmetric positive definite, and therefore the tensors $\varLambda$ and $\varSigma$ derived from $\boldsymbol{\textit {K}}^{-1}$ and $\boldsymbol{\textit {K}}^{-1/2}$ remain bounded and coercive. Consequently, the anisotropic permeability does not alter the admissible Sobolev framework and the variational formulation remains mathematically well posed.

Now, to derive the weak form, multiply (2.11)–(2.12) by the test functions $q \in Q$ and $\mathit{\boldsymbol {w}} \in W$ , and integrate over $\varOmega$ ,

(2.21) \begin{equation} \int _\varOmega (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{V}}) q \, {\rm d}\varOmega = 0, \end{equation}

(2.22) \begin{align} \begin{split} \int _\varOmega P \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega +\frac {M}{\epsilon \textit{Re}}\int _\varOmega ({\nabla} ^2 \boldsymbol{{V}}) \boldsymbol{\cdot }\boldsymbol{{w}} \, d\varOmega -\frac {1}{\textit{ReDa}}\int _\varOmega \left (\varLambda \boldsymbol{{V}}\right )\boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega \\-\frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} \int _\varOmega \left (\varSigma \boldsymbol{{V}}|\boldsymbol{{V}}|\right ) \boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega =0. \end{split} \end{align}

Since, the inlet condition is non-homogeneous, decompose the solution $\boldsymbol{{V}}$ into two parts, $\boldsymbol{{V}}=\boldsymbol{{V}}_{\!h}+\boldsymbol{{V}}_{\!P}$ , where $\boldsymbol{{V}}_{\!h}$ is the homogeneous solution space, which vanishes on Dirichlet boundaries, and $\boldsymbol{{V}}_{\!P}$ is a particular function that satisfies the non-homogeneous boundary conditions, i.e. $\boldsymbol{{V}}_{\!P}=(1,0)$ at the inlet.

Now, substitute the decomposed solution into (2.22) and integrating by parts the diffusion term, this equation can be written as

(2.23) \begin{equation} \chi (\boldsymbol{V}_{\!h},P,\boldsymbol{w})=\mathcal{L}(\boldsymbol{w}). \end{equation}

In (2.23), $\chi$ is a semi-linear operator and $\mathcal{L}$ is a linear operator, given by

(2.24) \begin{align} \begin{split} \chi (\boldsymbol{{V}}_{\!h},P,\boldsymbol{{w}}) &= -\int _\varOmega P \boldsymbol{\nabla }. \mathit{\boldsymbol {w}}\,{\rm d}\varOmega +\frac {M}{\epsilon \textit{Re}} \int _\varOmega \boldsymbol{\nabla }\boldsymbol{{V}}_{\!h} . \boldsymbol{\nabla }\boldsymbol{{w}}+\frac {1}{\textit{ReDa}}\int _\varOmega \left (\varLambda \boldsymbol{{V}}_{\!h}\right ) .\boldsymbol{{w}} \,{\rm d}\varOmega \\[3pt] &+\frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} \int _\varOmega \left (\varSigma (\boldsymbol{{V}}_{\!h}+\boldsymbol{{V}}_{\!P})|\boldsymbol{{V}}_{\!h}+\boldsymbol{{V}}_{\!P}|\right ).\boldsymbol{{w}}\, {\rm d}\varOmega , \end{split} \end{align}

(2.25) \begin{align} \begin{split} \mathcal{L}(\boldsymbol{{w}})&=- \frac {M}{\epsilon \textit{Re}}\int _\varOmega \boldsymbol{\nabla }\boldsymbol{{V}}_{\!P} . \boldsymbol{\nabla }\boldsymbol{{w}}-\frac {1}{\textit{ReDa}}\int _\varOmega \varLambda \boldsymbol{{V}}_{\!P}.\mathit{\boldsymbol {w}}\, {\rm d}\varOmega . \end{split} \end{align}

3.1. Step 3. Pressure recovery

Given the solution $\boldsymbol{{V}}_{\!h} \in W$ obtained from the weak formulation, the pressure field can be recovered from the momentum equation in weak form. In particular, the weak formulation $\chi$ in (2.23) contains the pressure contribution through the term $\int _\varOmega P \, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega ,$ where $\boldsymbol{{w}} \in W$ is an admissible test function. Rearranging the weak formulation yields the following variational problem for the pressure:

(3.31) \begin{align} \begin{split} \int _\varOmega P \, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega &= \frac {M}{\epsilon \textit{Re}} \int _\varOmega \boldsymbol{\nabla }\boldsymbol{{V}}_{\!h} . \boldsymbol{\nabla }\boldsymbol{{w}} \, {\rm d}\varOmega + \frac {1}{\textit{ReDa}} \int _\varOmega \varLambda \boldsymbol{{V}}_{\!h} \boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega \\[3pt] &\quad + \frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} \int _\varOmega \varSigma (\boldsymbol{{V}}_{\!h}+\boldsymbol{{V}}_{\!P}) |\boldsymbol{{V}}_{\!h}+\boldsymbol{{V}}_{\!P}| \boldsymbol{\cdot }\boldsymbol{{w}} \, {\rm d}\varOmega - \mathcal{L}(\boldsymbol{{w}}) \end{split} \end{align}

for all $\boldsymbol{{w}} \in W$ .

The existence and uniqueness of the pressure field are ensured by the Ladyzhenskaya–Babuška–Brezzi (LBB) or inf-sup condition (Brezzi Reference Brezzi1974), given by

(3.32) \begin{equation} \inf _{q \in Q} \, \sup _{\boldsymbol{{w}} \in W} \frac { \int _\varOmega q \, (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{w}}) \, {\rm d}\varOmega }{ \|\boldsymbol{{w}}\|_{H^1(\varOmega )} \, \|q\|_{L^2(\varOmega )} } \geqslant \beta \gt 0, \end{equation}

where $Q$ denotes the pressure space and $\beta$ is a positive constant independent of the discretisation parameter. This condition guarantees the stability of the mixed formulation and ensures that the pressure solution is uniquely determined up to an arbitrary additive constant.

Since pressure in incompressible flow problems is defined only up to a constant, an additional normalisation condition is imposed to obtain a unique solution. In the present work, the pressure is normalised by enforcing the zero-mean condition $ \int _\varOmega P \,{\rm d}\varOmega = 0.$

4.1. Neural network architecture

The computational framework is constructed using a multi-branch neural network architecture specifically designed to effectively handle the coupled flow variables and their complex interactions within the governing system (Biswas & Anand Reference Biswas and Anand2023; Aich & Bhargavi Reference Aich and Bhargavi2026b ). The architecture begins with a feature-learning module, which is herein referred to as the trunk, and consists of three fully connected layers with dimensionalities 128, 64 and 64 with hyperbolic tangent activation function. This is due to the fact that the first and relatively broad layer enables the network to capture large spatial characteristics of the flow, but the successive smaller layers are used to finesse these representations without generating unnecessary computational costs.

From this common latent representation, three separate output branches are constructed to independently predict the axial velocity $u$ , transverse velocity $v$ and pressure field $p$ . Each output branch consists of two hidden layers containing 20 neurons to effectively capture the coupled nonlinear flow behaviour; this is enough to provide the network with adequate expressiveness to represent the nonlinear variations without presenting a huge number of parameters to manage overall. The regular use of the tanh activation in the entire network guarantees smooth gradients; thus, making it easy to enforce governing equations using automatic differentiation. A schematic of the architecture is shown in figure 2.

Figure 2.

Architecture of the PINN.

4.2. Computational domain and symmetry

The hydrodynamic entry length in internal flows strongly depends on the Reynolds number. Therefore, the axial coordinate is normalised as $X^+=X/\textit{Re}$ to preserve physical consistency across different flow regimes. Owing to the geometric and flow symmetry about the channel centreline $(Y=0)$ , only the upper half of the computational domain is considered, thereby reducing the computational cost. The dimensionless hydrodynamic development length is taken as $X^+_{fd}=x_{fd}/\textit{LRe}=0.4$ , resulting in a computational domain defined over $X^+\in [0,0.4]$ in the streamwise direction and $Y\in [0,1]$ across the half-channel height.

The full velocity and pressure fields are reconstructed after the simulation by applying symmetric and antisymmetric conditions:

(4.1) \begin{equation} \begin{aligned} U(X^+,-Y) = U(X^+,Y),\quad V(X^+,-Y) = -V(X^+,Y),\quad P(X^+,-Y) = P(X^+,Y). \end{aligned} \end{equation}

This approach reduces the number of collocation points required for training without compromising accuracy.

4.3. Training strategy

The network is trained using a two-stage optimisation procedure designed to manage the nonlinearity of the governing equations:

1. (i) Adam phase (20 000 iterations). Training begins with the Adam optimiser, initialised with a learning rate of $1 \times 10^{-3}$ . During the first 10 000 iterations, a ramping parameter $\lambda \in [0,1]$ is gradually increased to introduce the Forchheimer nonlinear term into the momentum equation. This gradual activation prevents abrupt changes in the loss landscape and stabilises convergence as the model transitions from a linear to a fully nonlinear regime.

2. (ii) L-BFGS phase (250 iterations). After the initial optimisation, the solution is refined using the limited-memory Broyden–Fletcher–Goldfarb–Shanno method (L-BFGS). L-BFGS is well suited for PINNs because it converges rapidly when the solution is already close to a minimum and it does so with minimal memory usage. This final stage reduces the residuals more precisely and yields a smoother, physically consistent solution.

4.4. Loss functions

The composite loss function carefully balances multiple physical constraints,

(4.2) \begin{equation} \textit{MSE}= \textit{MSE}_{\textit{PDE}} + 10*\textit{MSE}_{\textit{BC}}, \end{equation}

where $\textit{MSE}_{\textit{PDE}}$ and $\textit{MSE}_{\textit{BC}}$ are the mean squared errors with respect to the partial differential equations of the system and the boundary conditions, respectively. Let us now consider a feed-forward deep neural network (DNN), $N(\boldsymbol{x}; h)$ , parametrised by $h$ , where $\boldsymbol{x}$ represents the input vector $(X^+, Y)$ , and $h$ represents the set of weights and biases of the neural network. This neural network takes the inputs $X^+$ and $Y$ , and produces the outputs $U$ , $V$ and $P$ , represented by $N_U$ , $N_V$ and $N_P$ , respectively. Then, $\textit{MSE}_{\textit{PDE}}$ and $\textit{MSE}_{\textit{BC}}$ are given by

(4.3) \begin{equation} \textit{MSE}_{\textit{PDE}}=\textit{MSE}_{\textit{cont}}+\textit{MSE}_{x\text{-}mom}+\textit{MSE}_{y\text{-}mom}, \end{equation}

(4.4) \begin{equation} \textit{MSE}_{\textit{cont}} = \frac {1}{N_{\textit{PDE}}} \sum _{i=1}^{N_{\textit{PDE}}} \left ( \frac {1}{\textit{Re}}\frac {\partial N_U}{\partial X^+} + \frac {\partial N_V}{\partial Y} \right )^2, \end{equation}

(4.5) \begin{align} \begin{split} \textit{MSE}_{x\text{-}mom} &= \frac {1}{N_{\textit{PDE}}} \sum _{i=1}^{N_{\textit{PDE}}} \Bigg ( \frac {1}{\textit{Re}} \frac {\partial N_P}{\partial X^+} - \frac {M}{\epsilon \textit{Re}} \left ( \frac {1}{\textit{Re}^2} \frac {\partial ^2 N_U}{\partial (X^+)^2} + \frac {\partial ^2 N_U}{\partial Y^2} \right ) \\ & - \frac {1}{\textit{ReDa}} (aN_U+bN_V) - \frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} ({\rm d}N_U+eN_V)\sqrt {(N_U)^2+(N_V)^2} \Bigg )^2, \end{split} \end{align}

(4.6) \begin{align} \begin{split} \textit{MSE}_{y\text{-}mom} &= \frac {1}{N_{\textit{PDE}}} \sum _{i=1}^{N_{\textit{PDE}}} \Bigg ( \frac {\partial N_P}{\partial Y} - \frac {M}{\epsilon \textit{Re}} \left ( \frac {1}{\textit{Re}^2} \frac {\partial ^2 N_V}{\partial (X^+)^2} + \frac {\partial ^2 N_V}{\partial Y^2} \right ) \\ & - \frac {1}{\textit{ReDa}} (bN_U+cN_V) - \frac {1.75}{\sqrt {150}\epsilon ^{3/2}\sqrt {\textit{Da}}} (eN_U+fN_V)\sqrt {(N_U)^2+(N_V)^2} \Bigg )^2, \end{split} \end{align}

(4.7) \begin{equation} \textit{MSE}_{\textit{BC}} = \frac {1}{N_{\textit{BC}}} \sum _{i=1}^{N_{\textit{BC}}} \left | N(\boldsymbol{\textit {x}}_{\textit{BC}}; h) - \mathcal{B} \right |^2. \end{equation}

Here, (4.4)–(4.6) indicate the losses associated with the continuity equation, $x$ -momentum equation and $y$ -momentum equation, respectively. Here, $N_{\textit{PDE}}$ and $N_{\textit{BC}}$ refer to the number of collocation points inside the domain and on the boundary, respectively, while $\boldsymbol{\textit {x}}_{\textit{BC}}$ represents the collocation points on the boundary. The term $\mathcal{B}$ denotes the prescribed boundary condition values corresponding to the neural network outputs at the boundary locations. Additionally, $ a = \sin ^{2}\alpha + K\cos ^{2}\alpha ,\ b = ({1}/{2})(K-1)\sin 2\alpha ,\ c = K\sin ^{2}\alpha + \cos ^{2}\alpha ,\ d = \sin ^{2}\alpha + \sqrt {K}\cos ^{2}\alpha ,\ e = ({1}/{2})(\sqrt {K}-1)\sin 2\alpha \text{ and } f = \sqrt {K}\sin ^{2}\alpha + \cos ^{2}\alpha$ are anisotropic constants.

The computational framework leverages PyTorch’s automatic differentiation capabilities for accurate gradient computation through the neural network and PDE residuals. The weights and biases of the network were initialised randomly. For training, 10 000 collocation points were used inside the domain, 2000 points were used near the inlet region, 1000 points were used along inlet and outlet boundaries, and 500 points were used along each wall.

The computational efficiency of the proposed framework was also assessed by comparing it with conventional numerical approaches previously employed for anisotropic DBF flows. In earlier FVM simulations reported by Aich & Bhargavi (Reference Aich and Bhargavi2026a ), the computational time required for solving similar anisotropic porous channel flow problems was approximately $10\,000$ s. However, the present PINN implementation was able to obtain converged solutions with NVIDIA T4 GPU acceleration in approximately 2000–2400 s. The low computational cost is mainly due to the fact that the PINN does not require any generation of the mesh, iterative pressure-velocity coupling processes and repeated discretisation steps that are involved in conventional methods. In addition, the governing equations and boundary conditions are directly included in the loss function via automatic differentiation, allowing for an efficient treatment of highly coupled nonlinear systems with anisotropic permeability and Forchheimer resistance effects. Another benefit of the proposed approach would be that the computational structure could be easily adapted for complicated geometries and parameter changes, without significant changes in the overall structure. In contrast to numerical solvers, PINNs can generate continuously differentiable solutions over the entire domain. For multiphysics porous medium simulations with nonlinear transport phenomena and anisotropic flow behaviour, these features prove to be very appealing for the present approach.

Figure 3.

Training loss histories for different Da and $\epsilon$ values.

Figure 4.

Training loss histories for different $K$ and $\alpha$ values.

Figure 5.

Training loss histories for different $M$ values.

Figures 3–5 display the training history for different values of the physical parameters. The optimisation was carried out in two stages. First, the models were trained using a version of the Adam optimiser, without mini-batching. There was a reduction in the loss initially, but in later epochs, it did not improve significantly. However, during the subsequent optimisation stage using the L-BFGS optimiser, the loss value decreased significantly. The second phase resulted in a substantial improvement in the accuracy of the solution, showing that L-BFGS is a very effective method for fine tuning the solution after the initial training using Adam. The observed behaviour shows that Adam offers a good starting point for convergence, whereas L-BFGS is good for final optimisation.

5.1. Validation with exact solution

For a fully developed flow, the velocity profile is given by $ \boldsymbol{{V}} = (U(Y),0)$ . By neglecting the nonlinear Forchheimer term, the governing equation simplifies to

(5.1) \begin{equation} 0 = -\frac {{\rm d}P}{{\rm d}X} + \frac {M}{\epsilon \textit{Re}} \frac {{\rm d}^2 U}{{\rm d}Y^2} - \frac {a}{\textit{ReDa}} U. \end{equation}

The general solution to this second-order ordinary differential equation (ODE) is

(5.2) \begin{equation} U(Y) = C_1 e^{\lambda Y} + C_2 e^{-\lambda Y} - \frac {\textit{Re}\,\textit{Da}}{a}\frac {{\rm d}P}{{\rm d}X}, \end{equation}

where

(5.3) \begin{equation} \lambda = \sqrt {\frac {a\epsilon }{M\,\textit{Da}}}. \end{equation}

Applying the boundary conditions $ U(\pm 1) = 0$ ,

(5.4) \begin{equation} C_1 e^{-\lambda } + C_2 e^{\lambda } = \frac {\textit{Re}\,\textit{Da}}{a}\frac {{\rm d}P}{{\rm d}X} \quad \text{and} \quad C_1 e^{\lambda } + C_2 e^{-\lambda } = \frac {\textit{Re}\,\textit{Da}}{a}\frac {{\rm d}P}{{\rm d}X}. \end{equation}

Solving for $ C_1$ and $ C_2$ ,

(5.5) \begin{equation} C_1 = C_2 = \frac {\textit{Re}\,\textit{Da}}{2a\cosh \lambda } \frac {{\rm d}P}{{\rm d}X}. \end{equation}

Substituting back, the exact solution for $ U(Y)$ becomes

(5.6) \begin{equation} U(Y) = -\frac {\textit{Re}\,\textit{Da}}{a}\frac {{\rm d}P}{{\rm d}X} \left ( 1-\frac {\cosh (\lambda Y)} {\cosh (\lambda )} \right )\!. \end{equation}

The average velocity is computed as

(5.7) \begin{equation} \overline {U} = \frac {1}{2}\int _{-1}^{1} U(Y)\,{\rm d}Y = -\frac {\textit{Re}\,\textit{Da}}{a}\frac {{\rm d}P}{{\rm d}X} \left ( 1-\frac {\tanh (\lambda )}{\lambda } \right )\!. \end{equation}

The normalised velocity $ U_a(Y)$ is then

(5.8) \begin{equation} U_a(Y) = \frac {U(Y)}{\overline {U}} = \left ( 1-\dfrac {\tanh (\lambda )}{\lambda } \right )^{-1} \left ( 1-\dfrac {\cosh (\lambda Y)} {\cosh (\lambda )} \right )\!. \end{equation}

Figure 6 compares the analytical solution $ U_a(Y)$ with the corresponding results from the PINN at $ X^+ = 0.4$ , demonstrating excellent agreement and validating the numerical approach.

5.2. Validation with published literature

The PINN code developed for this study has been rigorously validated against established benchmark solutions from the literature (Hooman & Gurgenci Reference Hooman and Gurgenci2007). As demonstrated in figure 7, excellent agreement is observed between our computed axial velocities at various X locations and the reference data, confirming the accuracy and reliability of the computational approach. This validation provides strong confidence in the PINN methodology for solving the DBF equations.

Figure 6.

Comparison of the PINN solution with the analytical solution at the exit.

Figure 7.

Comparison of the PINN solution with Hooman & Gurgenci (Reference Hooman and Gurgenci2007) for different $X^+$ values.