On nonlinear coupled differential equations for corrugated backward facing step (CBFS) with circular obstacle: AI-neural networking

1.   Introduction

In fluid science, a mathematical formulation is really important. The study of fluid motion, which includes both liquid and gas motion, is known as fluid dynamics, and it primarily uses mathematical equations and models to describe and forecast fluid behavior. The foundation of fluid dynamics is a set of mathematical formulas referred to as the flow equations. The conservation of momentum, mass, and energy in a fluid is explained by these equations. The Navier-Stokes equations, which are most frequently used to govern fluid flow, are derived from Newton's laws of motion and conservation laws. These formulas serve as the foundation for researching a variety of fluid phenomena and offer a mathematical depiction of fluid dynamics like Demarco and Dvorkin ^[1] developed mathematical models for the process of metal forming. Badly conditioned matrices were produced when those mathematical models were discretized using velocity or u-p interpolated finite element formulations with the incompressibility enforced by penalty measures. For the solution of flow equations, the application of an iterative solver was covered in this study. Jan ^[2] provided a mathematical overview of the key processes involved in numerically simulating the inviscid flow in a sealing gap of a screw compressor. The Euler equations described the mathematical model of compressible inviscid flow. The cell-centered finite volume design of the MacCormack method was employed for solution. Due to the fact that it combines a variety of components, including gas, liquids, powders, and lumps of granular materials, the ironmaking blast furnace was once thought to be among the largest and most complicated industrial reactors. Nogami et al. ^[3] devised the mathematical simulator for this process. Since several flow mechanisms controlled the motions of these materials, they employed the multi-fluid treatment. The model's validity was established and it was able to accurately replicate the temperature, velocity, and reaction fields in the furnace. Sankar and Hemalatha ^[4] offered a mathematical formulation for blood pulsatile flow via a catheterized artery. The power index and yield stress are the two parameters of the Herschel-Bulkley fluid. The solution for equations was obtained by the perturbation method. There was a discussion of how the non-Newtonian character of blood and catheterization affected the flow field. The yield stress value increases, and the plug core region's width does as well. When yield stress increases and other parameters remain constant, wall shear stress rises and, in contrast, flow rate and velocity drop.

Siddiqui et al. ^[5] investigated pulsatility and the non-Newtonian nature of blood through a stenosed artery. The fact that the viscous flow region's thickness varies with axial distance was interesting to note. The critical yield stress value at which the behavior of the flow rate transitions between the two types was established. In the investigation, the velocity profiles and related physiological traits were identified. In this analysis, numerous standard results about the flow of Newtonian and Casson fluids were obtained as special examples. Pandey et al. ^[6] offered a theoretical investigation of power law peristaltic flow in three layers with varying viscosities. Low Reynolds number and long wavelength approximations were used during the investigation. Comparable results were obtained when the intermediate layer's viscosity was maintained constant while the outermost layer's viscosity was raised. Nonetheless, understanding how the intermediate layer's viscosity affects things makes it easier for them to reach the necessary flow rate without upsetting the farthest layer. Increased flow rates are also favored by raising the flow behavior index. In order to get beyond the asymptotic nature and inherent mathematical difficulties of direct kinetic theories, this characteristic suggested looking for potential other accurate methodologies. Tessarotto and Cremaschini ^[7] looked into the basic mathematical characteristics of the NS phase-space dynamical system that underpins INSE and is identified by IKT in this work. Specifically, for the NS dynamical system, it was shown that a global existence theorem and correspondence theorem with the INSE problem hold. An analytical method was presented by Sumbatyan and Tarasov ^[8] to investigate the elastic wing vibration in non-viscous fluid. They start by creating a simple integral equation in two dimensions. After that, a succession expansion was used alongside the short coordinate. As a result, the problem is reduced to an endless set of integral equations, which were asymptotically examined in relation to the parameter for the big aspect ratio. To evaluate the practicality of the suggested approach, a demonstration of optimizing wing thickness was given. A mathematical model based on modified memory that describes fluid flow in porous media was presented by Obembe et al. ^[9]. The suggested mathematical model worked well to explain the atypical diffusion behavior seen in extremely heterogeneous and disordered porous media, as well as in a medium with a fractal shape. An existing discretization method served as the basis for the numerical strategy. The findings indicate that a decrease in the anomalous diffusion exponent's magnitude leads to an increase in wellbore pressure and pore pressure. Aleem et al. ^[10] examined a magnetized viscous flow field. Two methods, Caputo fractional and Caputo–Fabrizio fractional time derivatives, were used to build the fractional model. After transforming the equations into a dimensionless system, analytical solutions were obtained using the Laplace transform method. Thet observed how the temperature, concentration, and velocity fields were affected by fractional and other flow characteristics. Consequently, they discovered that temperature and velocity had higher outputs in the model with the fractional derivative of Caputo-Fabrizio than in the Caputo one. Additionally, velocity was seen to have a dual nature for both large and small-time ranges.

Patel et al. ^[11] offered mathematical formulation for micropolar fluid flow with pertinent physical effects. The system of partial differential equations (PDEs) was reduced to system of ordinary differential equations (ODEs) with the aid of appropriate similarity variables and the homotopy analysis method (HAM) was then used to solve the ODEs. Both tabular and graphical forms were used to discuss the influence of several parameters on flow field. Micropolar lowers the fluid velocity whereas mixed convection, buoyancy force, and permeability parameters improve the fluid velocity. The temperature was enhanced towards heat generation parameter; however, the Prandtl number and unsteadiness parameter showed the opposite trend. Hanafi et al. ^[12] offer a mathematical analysis of a hybrid nanofluid heat transfer capability when employing a cylindrical single-jet impingement scheme to transfer heat from a target surface. The influence of coolant types on heat transfer performance was examined. Multiple nanofluids with 0–5% volume concentration, and pure water were the coolants utilized in the simulation. A thorough presentation and discussion of the aspects of the Reynolds number subject to heat transfer performance was given. When compared to single-particle nanofluids and pure water the results demonstrated that hybrid nanofluids performed the best in terms of heat transfer. In order to improve the heat transfer system, Thenmozhi et al. ^[13] investigated the fluid behavior under a number of particular circumstances. Both homogeneous and heterogeneous chemical reactions were examined in this study using the Darcy-Forchheimer model. PDEs were used as the basis for the mathematical model. The solutions were reported using the shooting method in conjunction with the finite difference scheme. Increased fluid velocity and a slight drop in fluid temperature are observed at higher values of the microrotation parameter. The fluid's temperature rises while its velocity decreases due to the Darcy-Forchheimer parameter. Gholampour et al. ^[14] used Cine PC-MRI to measure the pulsatile blood velocity which was taken into consideration as input data towards mathematical modeling. It was determined that the pulsatile deformation of brain tissue subject to time was the inlet velocity in cerebrospinal fluid (CSF). Concentration, Navier-Stokes, and continuity were the governing equations in each of the three fields. The brain material characteristics were defined using the Darcy law which has defined permeability and diffusivity values. The accuracy of the CSF pressure and velocity was confirmed using mathematical formulations. To assess the features of the intracranial fluid flow, they used Womersley, Reynolds, Peclet and Hartmann numbers. CSF velocity was at its highest and CSF pressure was at its lowest during the cardiac cycles mid-systole phase. For the hydrocephalus patients and the healthy subjects, CSF pressure maximum and amplitude as well as CSF stroke volume were computed and compared. Raza et al. ^[15] suggested a mathematical model to examine how interfacial nanolayers affect the processes of mass and heat transfer in nanofluid flow. Additionally, the current study included the uniform transverse magnetic flux. The nonlinear equations were reduced to a collection of ODEs through the use of a similarity transformation. For this problem, a phase simulation based on thermophoresis and Brownian factors was created. In order to achieve the required accuracy, the shooting method was supported by the numerical procedure. Using the Casson parameter while considering the impact of thermal conductivity in nanolayers produced terrible outcomes. The thermal performance was greatly improved by the increase in nanolayer thickness.

Recently, a non-Newtonian mathematical model called Phan-Thien-Tanner was proposed by Akbar et al. ^[16] to study the flow inside the human body that is governed by cilia. Cilia work together to create a rhythmic motion that travels along the channel wall as a metachronal wave. Exact solutions for the cilia-governed flow problem were found by non-dimensionalizing the complex problem using the required assumptions. The wavy movements that resulted from the cilias rhythmic actions were visible in pressure gradient graphs. Additionally, this ciliated channel flow problem also showed a parabolic flow pattern. Several pumping regions were observed in the interesting pressure rise graphical solutions and streamlines reveal the trapping phenomenon. Fluid circulation in the brain nutrient digestion ovum transit in the reproductive system pathogen removal in various tissues and mucus clearance in the respiratory tract all depended heavily on cilia-generated peristaltic flow. The micropolar fluid flow in a bounded vessel as a boundary value problem (BVP) was proposed by Baranovskii et al. ^[17]. They considered rotational degrees for a fluids elementary volume and had couple of shear stresses. The couple stress fluids governing equations were of order four. As a result, they encounter the difficult problem of defining additional velocity endpoint conditions that make sense from a mathematical and physical perspective. They presented the considered BVP in both weak and strong formulations and examined how they related to one another. By using a broader version of the Leray-Schauder fixed-point theorem, they demonstrated the existence of a weak solution and its uniqueness subject to extra assumptions about the model data. They specifically examined how the built solutions converge to the stationary Navier-Stokes system solutions.

Deepalakshmi et al. ^[18] provided a mathematical analysis of the peristaltic pumping for Casson fluid through the debris-containing ureters inclined tapering wavy walls. Solid particles that were uniform and spherical were suspended in the Casson fluid which acts as the base liquid. The mathematical equations were utilized to model the motion of particles and fluids. In the transverse direction, an evenly distributed magnetic field was applied to the ureter flow. They used long wavelengths and a low Reynolds number to compute the flow-anchored equations. Using the Matlab software, the precise solutions were obtained in order to ascertain the fluids temperature contour pattern and particle velocity. The simulations show how electromagnetic radiation can be used to control urology pumping characteristics. A decrease peristaltic pumping rate is typically the result of an increase in particle volume fraction. Sanil et al. ^[19] conducted a new mathematical study on temperature dependent peristaltic Ree-Eying fluid flow. The study took into account a non-uniform channel with inclination and porous boundaries. Variable fluid properties that depend on temperature like changing viscosity and thermal conductivity were taken into account. The energy equations and the Navier-Stokes equation control the properties of fluid flow. The porous channel edge was subject to convective boundary conditions. The study took into account how chemical reactions affected the flow dynamics. A solution for velocity temperature concentration chemical reactions and stream functions was obtained by solving the non-linear coupled equations using the regular perturbation technique. Graphs were used to visually represent these variations after parametric analysis was applied. It was noticed that the velocity profiles decrease as the Ree-Eyring fluid parameter increases. Furthermore, there was an inverse relationship between the velocity profile and the velocity slip and porous parameters.

Owning to such mathematical formulation in fluid flow fields we offer mathematical analysis for Backward-facing step carrying viscous fluid. One should note that in CFD simulations and wind tunnel experiments, the BFS channels ^[20,21,22] are widely used to investigate the behavior of flow separation and its effect on aerodynamic performance. They aid in the optimization of the design of bluff bodies, such as cars, airplanes, and buildings, for increased aerodynamic efficiency by helping to understand the flow characteristics surrounding them ^[23,24]. Therefore, we consider viscous fluid flow with the parabolic pattern at the inlet and obstacle in the BFS channel. The flow is formulated mathematically in terms of PDEs and solved by using the finite element method. An AI-based neural networks model is constructed for the prediction of lift force faced by circular obstacles in the BFS channel. We are assured that the AI-based outcomes will be helpful for researchers to extend the idea for narrating the real-life applications of fluid flow around bodies such as, airplanes, buildings, and cars to mention just a few.

2.   Problem formulation

The backward-facing step (BFS) domain is considered in the present study. The inlet height is taken at 0.41m and the height of the outlet is taken at 0.61m. The expansion ratio is taken as 0.32m. The circular cylinder is taken as the obstacle with a diameter of 0.1m and it is centered in the BFS channel at (0.5, 0.4)m. The BFS is made corrugated by introducing triangular units at the bottom wall of the domain. The height of the triangular units is taken at 0.2m. The length of upstream and downstream is carried 0.5m and 1.6m respectively. The incompressible viscous fluid is introduced at an inlet. The said problem is two-dimensional and generally accepted fluid flow narrating differential systems are:

and the corresponding dimensionless form can be achieved by incorporating:

Therefore, Eq (2) turns into:

eliminating "*" results

In our case, the time independent viscous fluid flow is along the x-axis having velocity of $f(x, y)$ while the velocity $g(x, y)$ is along the y-axis. Therefore, Eq (5) reduces to the following forms:

The viscous fluid is taken at the inlet with a parabolic profile while it exits the BFS channel with Neumann conditions. The surface of the circular cylinder, upper and lower walls are taken with no-slip conditions. In this regard, the mathematical expressions as follows:

BFS inlet

BFS outlet

Obstacles/Side walls

The circular cylinder is installed in the BFS channel and striking of fluid flow will result in two hydrodynamic forces, namely drag and lift forces. The mathematical expressions as a drag coefficient (DC) and lift coefficient (LC) are:

It is important to note that the drag and lift coefficients offer a dimensionless, normalized way to compare and understand the forces being faced by natural and man-made obstacles across various flow situations and scales. Therefore, exploring coefficients instead of estimating forces results in more efficient modeling, flexibility in applications, and the ability to make predictions and general comparisons.

3.   Numerical methodology

In fluid science, various numerical techniques ^[25,26] are accessible to find the solutions of flow narrating differential systems, such as Hussain et al. ^[27], who recently used the Keller box method to report the solution of a system of flow equations for hybrid fluid flow with suspended nanoparticles.

Somiya and Kumar ^[28] utilized the bvp4c scheme to report the solution of a differential system constructed for the Maxwell flow field subject to a cylindrical surface. Salahuddin et al. ^[29] considered the Runge-Kutta scheme to find the solution of differential equations for the magnetized Williamson fluid flow towards a stretched cylinder. Most of the time, researchers used the reduction order (PDEs to ODEs) strategy to get the required solution, but for the best appropriate solution of Eqs (6)–(9), we adopted the finite element method.

The simulation is done with the assistance of Comsol Multiphysics ^[30,31]. For pressure and velocity, linear and quadratic polynomials are used. Further, Newton's method is carried out for the linearization of equations and PARDISO solved is adopted for inward approximation with convergence band of 10^-6.

3.1.   Meshing of BFS domain

The developed system of differential system is given as Eqs (6)–(9). For solution purposes used finite element method and the meshing of the field is the key step in this direction. The BFS is discretized by using hybrid meshing that includes triangular and rectangular elements. At level-1, the BFS domain contains 605 domain elements (DEs) and 91 boundary elements (BEs). Table 1 offers the complete description of the BFS domain subject to hybrid meshing. At level-6, the BFS domain carries 6996 DEs and 358 BEs. The corresponding meshing illustration is given in Figure 1. One can see that we have carried a fine meshing around the obstacle for better evaluation of hydrodynamic force.

3.2.   Grid independence test

The grid independence is tested by evaluating the lift coefficient via line integration around a circular obstacle. Table 2 is plotted in this regard. One can see that at a level for meshing level-7, 8, and 9, the LC seems independent of the choice of number of grid size. Therefore, for simulation purposes, we adopted the meshing level-9.

3.3.   Code validation

For validation of numerical scheme, we have considered as same configurations as carried out by Schafer et al. ^[32]. Both LC and DC are computed for circular obstacles at Re = 20, see Table 3. One can see that we found an excellent match with Schafer et al. ^[32], which yields the surety of present outcomes.

3.4.   Integration of FEM with multilayer perceptron (MLP) networks

In the present study, our interest is to construct the MLP networks for the prediction of lift coefficient subject to circular cylinder being placed in corrugated BFS. In the first step, the numerical values of the lift coefficient are obtained by using finite element-based solution. In the second step, such numerical values are used to develop MLP networks. In detail, the simulation is done at meshing level-9 by adopting the following particular values of flow parameter density $({r}_{d} = 1)$, viscosity $({a}_{1} = 0.001)$, characteristic length $({a}_{2} = 0.1)$, and mean inflow velocity and $({a}_{3} = 0.2)$, which is equivalent to say A = 20. The lift force is evaluated in terms of lift coefficient (LC) by doing line integration around a circular cylinder. Assisting mathematical relation in this regard is given in Eq (9). Tables 4–6 are evident in this direction. Particularly, Table 4 offers the LC values for a range of A = 60 to 91. Table 5 offers the statistics of LC for the range of A = 92 to 124. Table 6 gives the numerical values of LC towards A = 125 to 155. In an absolute sense, we have noticed that for advanced values of A, the LC shows inciting values. The lift force faced by the circular obstacle in the BFS domain is increasing the function of A.

4.   Results and discussion

The viscous fluid flow is considered in the BFS channel rooted with circular obstacle. The parabolic velocity profile is initiated at the inlet. The developed system of non-linear differential equations is solved by using FEM. The obtained results are provided as Figure 2. In detail, Figure 2(a) offers the velocity distribution of fluid flow in BFS domain rooted with circular cylinder.

The fluid enters with a parabolic pattern, and after striking with the obstacle, it gets bifurcated and travels downstream. The equivalent pressure outcome is given in Figure 2(b). The pressure is noticed maximum at the left face of the cylinder. The line graph study of the velocity of the fluid at different positions of the channel is given in Figure 2(c). We can see that at x = 0.0, the fluid enters with a parabolic profile while at x = 0.4, the fluid starts bifurcated. At x = 0.5, the fluid is completely bifurcated around the installed obstacle. At x = 0.9, the fluid is reforming the initial pattern. It is important to note that obtaining solutions for PDEs by using neural networks, such as physics-informed neural networks ^[33,34,35] and multilayer perceptron neural networks ^[36,37,38], gained significant attention. We constructed MLP neural networks to predict LC. It is well known that the Reynolds number (A) depends on viscosity, inflow velocity, density, and characteristics length. Therefore, for better prediction, we set all these flow parameters as inputs, namely ${r}_{d}, {a}_{1}, {a}_{2}, \text{and}\text{}{a}_{3}$, in the first layer of the ANN model. The LC is considered as an output in the last layer while 10 number of neurons are taken in the hidden layer.

The structure of the neural network is given in Figure 3(a). A total of 95 inputs are taken as input and 95 sample values are collected for LC as an output. The sample data is randomly distributed into three regimes: training, validation, and testing. Here, 67 (70%) values are used in training and 14 (15%) used for testing and validation. The Levenberg-Marquardt algorithm ^[39,40] utilizes gradient descent and the Gauss-Newton method. When the weights of models are near optimal values, it converges more quickly than traditional gradient-based methods. This can significantly reduce practice training time. Therefore, the ANN is trained by the Levenberg-Marquardt algorithm and the ultimate outcomes are offered in Figure 3(b) and (c). The activation functions used in hidden and output layers are as follows:

The performance of networks is tested by using mean square error (MSE) and coefficient of determination (R). The mathematical relation for MSE and R, are given as follows:

To be specific, Figure 3(b) gives the performance outcome for constructed ANN. We can observe that up to epoch 13, the graphical trends for training, validation, and testing approaches the lower error, and after epoch 13, all curves represent the asymptotic convergences to 10^-8. The best validation performance is noticed 7.6924e^-09 at epoch 25. The corresponding error histogram is given in Figure 3(c), which represents the lower error values. From both figures, we can conclude that ANN training is finished successfully.

The regression for all three stages is reported in Figure 4. For training, validation, and testing, the regression values are noticed 9.99996e^-01, 9.99997e^-01, and 9.99996e_,^-01 respectively. All values are approaches to R = 1. The regression value for all stages is also R = 1. The values show that the forecasted values and original values are strongly correlated. The corresponding performance parameters for the constructed ANN model are given in Table 7. The MSE values for each stage are up to mark, and hence the created ANN is the best to forecast the values of LC. The graphical illustration of EV for each data set of LC is presented in Figure 5(a). It can be observed that the error values for each prediction are very low, and the hence prediction made by ANN is trustful.

The graphical comparative study of predicted values of LC with targeted values of LC is given in Figure 5(b). One can see that predicted and targeted values are in great agreement. Therefore, the raised neural model is the best AI-based model that holds the capacity to estimate the lift force faced by the circular cylinder being installed in the BFS channel. The prediction by ANN can be extended to a wide range of Reynolds numbers.

5.   Conclusions

The corrugated backward-facing step (BFS) field is theoretically modeled in the presence of a viscous fluid. The circular cylinder is used as an obstacle in BFS. The flow field is analyzed by combining the finite element method and artificial intelligence-based MLP networks. The neural network model is used to predict the lift coefficient subject to circular obstacle. The main conclusions are as follows:

● The best validation performance for the built ANN model is 7.6924e^-09 at epoch 25. Particularly for training, validation, and testing, MSE values are 6.88504e^-09, 7.69240e^-09, and 1.123604e^-08, respectively.

● For training, validation, and testing, the regression values are 9.99996e^-01, 9.99997e^-01, and 9.99996e^-01, respectively. As a result, there is a significant correlation between the anticipated and target LC values.

● According to the MSE and R statistics, the constructed neural networks model is the best model to forecast LCs across a wide range of Reynolds numbers.

● In an absolute sense, the lift force experienced by a circular obstacle rooted in a backward-facing step field is found to be an increasing function of a higher Reynolds number.

Author contributions

K. U. Rehman: Writing – original draft, validation, investigation, formal analysis; W. Shatanawi: Writing – review & editing, visualization, investigation; supervision; Weam G Alharbi: Writing – review & editing, numerical solution, visualization, validation, conceptualization. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The authors would like to thank Prince Sultan University, Saudi Arabia, for the technical support through the TAS research lab.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.