Particle Swarm Optimization for exploring Darcy–Forchheimer flow of Casson fluid between co-axial rotating disks with the Cattaneo–Christov model

In this paper, we carried out a numerical analysis of the fluid dynamics and heat transfer occurring between two parallel disks. The study accounts for the impact of temperature-dependent fluid viscosity and thermal conductivity. We systematically investigated various parameters, including viscosity, thermal conductivity, rotational behavior (rotation or counter-rotation), and the presence of stretching, aiming to comprehend their effects on fluid velocity, temperature profiles, and pressure distributions. Our research constructs a mathematical model that intricately couples fluid heat transfer and pressure distribution within the rotating system. To solve this model, we employed the 'Particle Swarm Optimization' method in tandem with the finite difference approach. The results are presented through visual representations of fluid flow profiles, temperature, and pressure distributions along the rotational axis. The findings revealed that the change in Casson factor from 2.5 to 1.5 resulted in a reduction of skin friction by up to 65%, while the change in local Nusselt number was minimal. Furthermore, both the viscosity variation parameter and thermal conductivity parameters were found to play significant roles in regulating both skin friction and local Nusselt number. These findings will have practical relevance to scientists and engineers working in fields related to heat management, such as those involved in rotating gas turbines, computer storage devices, medical equipment, space vehicles, and various other applications.


Problem statement and formulation
Consider two disks each with radii r occupy the planes z = 0 and z = h and are rotating with angular velocity 1 and 2 about z-axis in the same sense.A steady, incompressible and Darcy-Forchheimer flow of Casson fluid confined between the parallel disks.Both the disks are subject to radial stretching with stretching rates A 1 and A 2 and are assumed to be heating convectively.In this study, we considered that the flow model is subject to non-linear Boussinesq approximation (see 5th term on RHS of Eq. ( 3)).Owing to the practical utilization of non-linear Boussinesq approximation, the heating effect due to non-linear thermal radiation is accounted while heating due to porous medium is neglected.Here, the dynamic viscosity and thermal conductivity are assumed to be temperature dependent and other properties are independent of temperature.
The geometry illustrating the problem is shown in Fig. 1.The rheological equation which describes the relation of shear stress and strain rate for the incompressible flow of Casson fluid can be written as Upreti et al. 20 : In case of Casson fluid (non-Newtonian fluid), where π > π c , it is possible to say that µ = µ B + (2π) −1/2 P z and ϑ = µ B ρ 1 + 1 δ where δ = µ B √ 2π(P z ) −1 termed as Casson factor.Here, symbols Ŵ ij , µ , µ B , ϑ, P z , π = ς ij ς ij , π c and ς ij denotes Casson yield stress, absolute viscosity, plastic viscosity of non-Newtonian fluid, kinematic viscosity of Casson fluid, yield stress of the fluid, product of deformation rate, critical value of the product and component of deformation rate, respectively.
Under the above-stated assumptions along with the non-linear Boussinesq approximation 54 and the Cattaneo-Christov heat flux model 49,53 , the current problem is governed as 11,25 : (1) The flow is induced by the radial stretching or shrinking movement of the disks; also, the disks are supposed to be rotating along their axis of rotation.These result in centripetal and centrifugal forces, which in turn generate fluid flow.Additionally the disks undergo convective heating, therefore the respective boundary conditions are 44,49 : In above expressions, the symbols u, v and w denotes velocity components (m/s) in increasing direction of r, θ and z,ρ density kgm −3 of working fluid , C P specific heat J/kgK , k fH Forchheimer permeability, F o inertial coefficient of porous medium, g gravitational acceleration (m s −2 ), β a and β b linear and non-linear thermal expansion coefficients K −1 , p pressure (Pa),τ heat flux relaxation time, σ * Stefan Boltzmann constant Wm 2 K −4 , K * mean absorption coefficient m −1 , a and b coefficients of convective heat transfer Wm 2 K −1 , T fluid temperature (K) in the boundary layer.
Moreover, the dynamic viscosity (µ B (T)) and thermal conductivity (κ(T)) of the working fluid are considered being depending on temperature 10,55 : Here, µ * B , γ , κ ∞ and ε refers to ambient fluid viscosity, non-negative constant, ambient fluid thermal con- ductivity and thermal conductivity parameter, respectively.Now defining the following similarity transformation variables [10][11][12] Here, ς is the pressure parameter and ξ is the non-dimensional distance along the axis of rotation.
(3)  8) and ( 9) and similarity variables (11), the continuity equations i.e., Eq. ( 2) is satisfied identically and the Eqs.( 3)-( 6) are transformed to following dimensionless form and boundary conditions (7) transformed to following forms: Here, the symbols ̟ , Fr, Re, E, 1 , 2 , Ra, Pr, τ T , θ r , (A a , A b ), � , and (Bi a , Bi b ) present in the above expressions are refers to porosity parameter, Forchheimer parameter, Reynolds number, viscosity variation parameter, thermal buoyancy parameter, non-linear convection parameter, radiation parameter, Prandtl number, thermal relaxation parameter, temperature difference parameter, stretching parameter for lower and upper disks, rotation parameter, Biot numbers for lower and upper disks, respectively and are defined as Now, eliminating the pressure parameter ς from Eq. ( 11), and for this differentiating Eq. ( 11) with respect to ξ , we get the following equation Now, making use of (11) and (15), the following form of pressure parameter is obtained: ) To calculate the pressure term, Eq. ( 13) undergo integration with respect to ξ from the limits 0 to ξ , under the assumption of P(ξ = 0) = 0 we have The skin friction coefficient Cf x and local Nusselt number (Nu x ) at lower (ξ = 0) and upper (ξ = 1) disks are given as Khan et al. 44 where Re r = ϑ r� 1 h is local Reynolds number.

Numerical methodology
The non-dimensional mathematical model, detailed in Eqs. ( 12), ( 14) and ( 17) along with the boundary conditions (15), has been effectively addressed through numerical methods.Initially, we restructured all equations into a system of eight first-order differential equations, adhering to the specified boundary conditions (Eq.15).However, for solving this system using a finite difference scheme, all eight initial conditions are required.Unfortunately, the formulated problem only provides four initial conditions (Eq.15).Hence, to resolve this system, we must compute the remaining four missing initial conditions using the four known boundary conditions, also delineated in Eq. ( 15).To achieve this, we employ the "Particle Swarm Optimization" method in conjunction with the Runge-Kutta-Fehlberg scheme.We establish an error function, termed the root of residuals (RR), which is minimized to obtain the missing boundary condition and hence the final solution of the problem.A comprehensive description of this methodology can be found in 20 .
In PSO the swarms are updated by using the following equations.
The overall objective of the proposed method is to determine the initial such that the following error function is minimized.
The process is implemented with the help of computer code written in MATLAB.The minimization of this "RR" is achieved with the help of PSO under the convergence criteria of RR ≤ 10 −08 .The values of PSO conver- gence controlling parameters are taken as φ p = φ p = 1.49618 and ψ = 0.7298 .
The flow chart for the computational procedure is given below. (

Results and discussion
The current section deals with the visualization of the influence of parameters like A a , A b , Re, Bi a and on the fluid's velocity (axial F(ξ ) , radial F ′ (ξ ) and tangential G(ξ ) ), temperature θ(ξ ) and pressure P(ξ ) ; all the flow attribute profiles are plotted against distance along the axis of rotation (ξ ) .In the computations, the acceptable range assigned to the parameters for better graphical results is as follows  (15).And, the quantities like SFC and LNN are computed for the lower and upper discs for the aforesaid parameters and presented through Tables 2 and 3.
Figure 2 shows that increasing the lower disk stretching parameter A a , leads to an increase in the axial velocity F(ξ ) for each .This could be attributed to the fact that as the stretching parameter increases, the lower disk's deformation causes changes in the flow pattern.This alteration may result in an increased F(ξ ) , possibly due to changes in the effective area through which the fluid flows.Moreover, it is noticed that axial velocity is higher for the case = 0 followed by = −1 and 1 .Figure 3 that the rise in A a results enhancement of the radial velocity F ′ (ξ ) at the lower disk, but near the upper disk the velocity gets declined.This can be explained by the stretching effect causing changes in the radial distribution of velocities.The enhancement at the lower disk could be due to increased surface effects influencing the radial flow, while near the upper disk, other factors might be dominating and causing a decline.Figure 4 implies that, for = −1 and 0 tangential velocity G(ξ ) decrease with growing A a parameter while a dual behavior is recorded for = 1 , here for ξ < 0.3 velocity decreases while it increases for ξ > 0.3 .Here, negative values of velocity indicate that movement of upper disk is faster than the lower disk.Figures 5 and 6 specify that a reduction is temperature profile while an escalation in pressure distribution is Table 1.Comparison of results ω= 0, Fr= 0, Re = 1, E = 0, ε = 0, 1 = 0 = 2 ,Pr = 6.2, τ = 1, θ r = 0, and A a = 0 = A b .
Present results Turkyilmazoglu 56 Present results Turkyilmazoglu 56  recorded for growing values of A a .This is because of the changes in flow patterns and energy dissipation associ- ated with the increased stretching effect.Figures 7, 8, 9, 10, 11 are plotted to highlight the impacts of altering A b on flow properties (velocities, tempera- ture and pressure distribution) for co-rotation ( � = 1) , rotor-stator ( � = 0) , and counter rotation ( � = −1) situations, respectively.Figures show that an increase in the values of A b from −1 to 1 , results in a decrease in the axial velocity (see Fig. 7).Physically, increasing A b causes changes in the geometry of upper disk and surface effects, also the disks are assumed to be rotating about their axis, thus results in a decrease in F(ξ ) .Radial veloc- ity declines when ξ < 0.6 while the opposite pattern is recorded for ξ > 0.6 (see Fig. 8), and tangential velocity G(ξ ) increases for = −1 and 0 cases, while dual behavior is observed for = 1 , for ξ > 0.8 which velocity decreases whereas it increases for ξ < 0.8 (see Fig. 9). Figure 10 denotes that the temperature profile escalates with growth in; here, for a particular fluid temperature, it doesn't show any significant variation for = 1, 0, −1 .Physically, increasing A b leads to enhanced mixing and dissipation of energy, resulting in an overall increase in the temperature profile.and Fig. 11 shows that for ξ < 0.7 the pressure is in a positive direction and increases with an increase in the values of A b , whereas it falls in the remaining region.The positive pressure and its increase Table 3. Computed values of SFC and LNN when δ = 0.8, ω = 0.5,                 in a specific region might be attributed to changes in the flow patterns and dynamic forces caused by alterations in A b .The falling pressure in other regions may result from the counteracting effects of the stretching parameter on the flow dynamics.Figures 12, 13, 14, 15, 16 are shown in order to visually depict the implications of different Re on flow charac- teristics, such as velocities, temperature, and pressure distribution, under the scenarios of co-rotation ( � = 1) , rotor-stator ( � = 0) , and counter-rotation ( � = −1) .Axial velocity F(ξ ) of the fluid increases with increase in values of Re for each value of but velocity correspond to ( � = 1) is superior than other cases, further near the lower and upper disks the variation in fluid velocity is negligible, as depicted in Fig. 12. Physically, higher Re typically indicate an increase in fluid inertia compared to viscous forces.This leads to more streamlined and faster flow, resulting in an increase in F(ξ ) .Moreover, the negligible variation in F(ξ ) near lower and upper disks may be due to the influence of boundary conditions, geometric constraints, or other factors that limit the impact of Reynolds number on the flow near these surfaces.Figure 13 illustrates that F ′ (ξ ) of the fluid is a monotonically increasing function of Re for ξ < 0.6 , and behave as decreasing function in the remaining region.Variation in tangential velocity due to distinct Re is shown in Fig. 14, figure shows that for = 0 and 1 , velocity reduces but for = −1 it reduces for ξ < 0.6 and thereafter it increases.Figure 15 shows that fluid is hotter at the surface of lower disk compared to upper disk, fluid temperature reduces with growth in Reynolds number from 2 to 10 for each case of while for ξ > 0.6 a contrasting trend is observed is noticed for temperature profile correspond to = −1 and 1 .And, pressure reduces with growing Re for ξ < 0.6 and in the remaining region it elevates.The reduction in pressure might be associated with enhanced kinetic energy and turbulence at higher Reynolds numbers, while the increase in pressure in other cases could be influenced by the specific rotational configurations and their effects on flow patterns and dynamic forces.Figures 17, 18, 19, 20, 21 show the effects of changing Bi a on flow properties (velocities, temperature, and pressure distribution) for co-rotation ( � = 1) , rotor-stator ( � = 0) , and counter-rotation scenarios.Figure 17 illustrates that axial velocity near the surface of the upper and lower discs has insignificant variation with distinct values of Bi a and , but for ξ > 0.2 velocity field increases with increasing Bi a , and velocity becomes negative for ξ ≥ 0.5 .The stability of axial velocity for lower values of ξ suggests that at these regions, the heat transfer effects do not significantly influence the axial flow.And, increase in velocity for ξ > 0.2 may be indicative of stronger interactions between heat transfer and fluid dynamics, leading to changes in axial flow patterns.The negative velocity for ξ ≥ 0.5 indicate a reverse flow or backflow in certain regions.Radial velocity F ′ (ξ ) and tangential velocity G(ξ ) reduce with growth in Bi a , see Figs. 18 and 19, respectively.The reduction in F ′ (ξ ) and G(ξ ) with increasing Bi a is likely a consequence of enhanced heat transfer effects, higher Bi a indicate a stronger coupling between the fluid and the solid surfaces, impacting the radial and tangential components of the flow.Figure 20 illustrate that near the lower disc, fluid is hotter compared to fluid at upper disc and fluid temperature upsurges with growth in Bi a ; similar behavior is recorded for pressure profiles see Fig. 21.Higher Bi a likely lead to varia- tions in pressure due to the intensified heat transfer effects on the fluid.Figures 22, 23, 24, 25, 26 plotted to illustrate the consequences of varying Re on flow features (velocities, temperature and pressure distribution) for distinct values of Bi a .Figure 22 shows that near the lower and upper discs the difference between axial velocities is insignificant for distinct Re , but for 0.2 < ξ < 0.4 it decreases with growth in Re while for 0.4 < ξ < 0.8 it rises.Figure 23 illustrate that, for a particular Re radial velocity decreases for ξ < 0.5 while it rises thereafter, however the distribution of radial velocity does not show any sig- nificant variation with distinct values of Re for ξ < 0.3 whereas for 0.3 < ξ < 0.6 it increases and for ξ > 0.75 it reduces.Figure 24 shows that, at both lower and upper disk, tangential velocity is higher, and velocity distribution decreases with increase in Reynolds number.Usually, higher tangential velocity near the discs is likely due to the rotational motion of the system, which is more pronounced near the rotating surfaces; and the decrease in velocity distribution with an increase in Reynolds number is because of the changes in the flow dynamics associated with different Reynolds numbers.Figure 25 shows that, near the lower disk, fluid have higher temperature, and as temperature falls as fluid moves away from lower disk, as it reaches the upper disk fluid temperature elevates, thus higher temperature near the lower disc indicates more efficient heat transfer from the solid surfaces to the fluid in that region.Moreover, fluid temperature corresponds to Bi a = 0.5 is higher compared to Bi a = 0.3 and 0.1 ; this suggest that the heat transfer characteristics are influenced by the thermal properties of the system.Pressure distribution profiles elucidate that it is minimum at the extreme ends i.e., lower and upper discs while it is maximum at the mid-way of the discs; moreover, it reduces with escalation in Re and distribution profile correspond to Bi a = 0.5 is higher followed by Bi a = 0.3 and 0.1 , respectively, see Fig. 26.

Table 4. Computed values of SFC and LNN when
It is noted from Table 2 that as the lower disc undergoes the transition from radial shrinking to stretching, the lower disc has a higher SFC, whereas the LNN is higher at the upper disc; the opposite pattern is recorded when the upper disc undergoes a similar transition while the lower disc is radially stretching.It is anticipated from Table 3 that, when discs are in co-rotation, the LNN at the lower disc is higher compared to the upper disc for a higher Reynolds number, whereas the upper disc has a higher SFC than the lower disc; on the other hand, for counter-rotation and no rotation, the lower disc has a higher SFC.For the case where the lower disc is exposed to convective heating, it is noted from Table 3 that for each case of , the upper disc has a higher local Nusselt number; however, the variation between LNN for these discs is less than 1%.A similar pattern is recorded for LNN when discs are in co-rotation and Re and Bi a are varying; here, both LNN and SFC increase with increasing magnitudes of Re and Bi a .
In the case of counter rotation of discs, Table 4 presents the variations of LNN (local Nusselt number) and SFC (skin friction coefficient) concerning the Casson factor, thermal conductivity parameter, and viscosity variation parameter.From the results, it is evident that increasing the Casson parameter from 0.25 to 1.5 leads to a significant reduction in SNF (skin friction) on both lower and upper discs, by up to 65%.However, the increase in LNN is minimal, being less than 0.06%.Similarly, augmenting the viscosity variation parameter from 0 to 0.2 results in a decrease in SFC of up to 13%.Conversely, an increase in the thermal conductivity parameter exhibits a notable impact on LNN, affecting both lower and upper discs, with a change of up to 17% in LNN with a variation in the thermal conductivity parameter from 0 to 0.6.

Summary and conclusion
In this investigation, we have conducted a comprehensive examination of how various parameters influence fluid dynamics and heat transfer between two parallel disks.Our study provides a visual representation of how these parameters impact fluid velocity, encompassing axial, radial, and tangential components, as well as temperature distribution and pressure profiles along the axis of rotation.The quantitative dimensions of our research is encapsulated in the analysis of key metrics, namely the Skin Friction Coefficient (SFC) and the Local Nusselt Number (LNN), for both the lower and upper disks.We have also explored the behavior of SFC and LNN as the lower and upper disks transition between radial shrinking and stretching.
From these investigations, several salient conclusions emerge: • The Casson parameter plays a crucial role in reducing the skin friction coefficients on the rotating disks to a very high extent, while the thermal conductivity parameter primarily impacts the heat transfer rate (LNN) with minimal effect on the skin friction coefficient.• During the transition from radial shrinking to stretching of the lower disc, we observed a higher Skin Friction Coefficient (SFC) for the lower disc, while the Local Nusselt Number (LNN) is elevated at the upper disc.Conversely, when the upper disc undergoes a similar transition while the lower disc is radially stretching, the roles are reversed, with the lower disc exhibiting a higher SFC and the upper disc displaying a higher LNN.
• In scenarios involving co-rotation, the LNN at the lower disc exceeds that at the upper disc for higher Reyn- olds numbers.However, it is consistently noted that the upper disc maintains a higher SFC than the lower disc under identical conditions.• In cases of counter-rotation and no rotation, the lower disc consistently demonstrates a higher SFC compared to the upper disc.• When convective heating is applied to the lower disc, the upper disc registers a higher Local Nusselt Number.
Nevertheless, the variation in LNN between these two discs remains minimal, at less than 1%.
These findings collectively contribute to a deeper understanding of the dynamic behavior of parallel rotating disks under diverse circumstances.They are of significant relevance to engineers and researchers working on systems involving rotating components, providing essential insights into optimizing performance and heat transfer management.
Finally, fluid flow between coaxial disks may have industrial uses, such as the construction of mixing devices, lubricating systems, and cooling systems.Future study might concentrate on improving these applications' efficiency, sustainability, and cost-effectiveness.

Figure 1 .
Figure 1.Schematic representation of the flow model.

Figure 2 .
Figure 2. Consequences of varying A a on axial velocity F(ξ ) for distinct .

Figure 3 .
Figure 3. Consequences of varying A a on radial velocity F ′ (ξ ) for distinct .

Figure 4 .Figure 5 .
Figure 4. Consequences of varying A a on tangential velocity G(ξ ) for distinct .

Figure 6 .
Figure 6.Consequences of varying A a on pressure P(ξ ) for distinct .

Figure 7 .Figure 8 .
Figure 7. Consequences of varying A b on axial velocity F(ξ ) for distinct .

Figure 9 .
Figure 9. Consequences of varying A b on tangential velocity G(ξ ) for distinct .

Figure 12 .
Figure 12.Consequences of varying Re on axial velocity F(ξ ) for distinct .

Figure 18 .
Figure 18.Consequences of varying Bi a on radial velocity F ′ (ξ ) for distinct .

Figure 21 .
Figure 21.Consequences of varying Bi a on pressure P(ξ ) for distinct .

Figure 22 .Figure 23 .
Figure 22.Consequences of varying Re on axial velocity F(ξ ) for distinct Bi a .

Figure 24 .Figure 25 .
Figure 24.Consequences of varying Re on tangential velocity G(ξ ) for distinct Bi a .

Figure 26 .
Figure 26.Consequences of varying Re on pressure P(ξ ) for distinct Bi a .
To validate the methodology, we solved a specific case of the presented problem and compared the results with those from previously published work.The comparison presented in Table1below serves to illustrate the credibility and validity of the proposed methodology.

Generate the random population for unknown initial conditions Generate the finite difference solution (using Runge-Kutta Fehberg method) for system of coupled non-linear differential equations (12, 14, 17) with ICs/BCs (15) Assign the random positions and velocities to the particles of randomly generated population and Determine fitness function (RR (23)) Initialize, Pbest, and determine Gbest Update velocities and positions according to equations (21) & (22) Evaluate the fitness of each particle (RR) and update Pbest and Gbest Convergence check No Yes Optimum value of RR and unknown initial conditions achieved Terminate the process Iter>MaxIter Yes No obtained
5and other parameters are treated as constant quantities.From the graphs, it is clear that the solutions satisfy both the initial and boundary conditions

Table 2 .
Computed values of SFC and LNN when δ