MHD Mixed Convection of Non-Newtonian Bingham Nanoﬂuid in a Wavy Enclosure with Temperature-Dependent Thermophysical Properties: A Sensitivity Analysis by Response Surface Methodology

: The numerical investigation of magneto-hydrodynamic (MHD) mixed convection ﬂow and entropy formation of non-Newtonian Bingham ﬂuid in a lid-driven wavy square cavity ﬁlled with nanoﬂuid was investigated by the ﬁnite volume method (FVM). The numerical data-based temperature and nanoparticle size-dependent correlations for the Al 2 O 3 -water nanoﬂuids are used here. The physical model is a two-dimensional wavy square cavity with thermally adiabatic horizontal boundaries, while the right and left vertical walls maintain a temperature of T C and T H , respectively. The top wall has a steady speed of u = u 0 . Pertinent non-dimensional parameters such as Reynolds number ( Re = 10,100, 200,400), Hartmann number ( Ha = 0,10, 20), Bingham number ( Bn = 0,2,5,10,50,100,200), nanoparticle volume fraction ( φ = 0,0.02,0.04), and Prandtl number ( Pr = 6.2) have been simulated numerically. The Richardson number Ri is calculated by combining the values of Re with a ﬁxed value of Gr , which is the governing factor for the mixed convective ﬂow. Using the Response Surface Methodology (RSM) method, the correlation equations are obtained using the input parameters for the average Nusselt number ( Nu ), total entropy generation ( E s ) t , and Bejan number ( Be avg ). The interactive effects of the pertinent parameters on the heat transfer rate are presented by plotting the response surfaces and the contours obtained from the RSM. The sensitivity of the output response to the input parameters is also tested. According to the ﬁndings, the mean Nusselt numbers ( Nu ) drop when Ha and Bn are increased and grow when Re and φ are augmented. It is found that ( E s ) t is reduced by raising Ha , but ( E s ) t rises with the augmentation of φ and Re . It is also found that the φ and Re numbers have a positive sensitivity to the Nu , while the sensitivity of the Ha and Bn numbers is negative.


Introduction
When buoyancy forces are added to a pressurized flow, mixed convection (which involves both free and forced convection) occurs [1]. In recent years, the mixed convection in fluid-filled lid-driven cavities has attracted much interest due to its significance in the field of heat and mass transfer and its multiple applications in engineering and sciences [2][3][4]. For example, applications include food processing, solar heat transfer, crystal growth, thermal storage system, cooling ponds, and glass production [5]. Two areas of the metal field are also distinguished: the first one is an unyielding zone where the fluid moves rigidly, and the second one is where the fluid moves like a viscous liquid. The second characteristic of the extra stress is greater than the yield stress in the yielded region, whereas it is equal to or lower than the yield stress in the unyielded zone. According to Migórski and Dudek [6], below a stress function's yield point, the medium behaves like a rigid body, and above the limit, it flows like an incompressible viscous fluid. As a result, determining the position and form of the yield surface(s), or the interaction between these two sets, is necessary for solving flow issues involving these fluids. Therefore, we employ a novel numerical method based on the finite volume method (FVM) to resolve the flow of a Bingham fluid in a lid-driven wavy square cavity with a heated left wall.
It is challenging to investigate the properties of a particular non-Newtonian fluid classified as viscoplastic in fluid mechanics. In our surroundings, several viscoplastic materials are used in geological and industrial processes. These compounds include suspensions, drilling mud, paste, and food items such as mayonnaise and ketchup. The fundamental characteristic that sets this fluid apart is the division of the flow domain into yielded and unyielded regions. At low pressures, it behaves like a solid object, but at considerable pressures, it flows like a viscous fluid [7]. Yield stress is the crucial stress value below which the material does not deform and over which it does. The properties of such non-Newtonian fluids can be modeled mathematically using the Bingham model [8,9], the Herschel-Bulkley model [10], and the Casson model [11]. The Bingham approach is the one that best captures the properties of viscoplastic material. Using numerical methods to solve the model might be challenging because it is discontinuous. The Papanastasiou regularization technique [12], a different strategy suggested by Papanastasiou, can be used to overcome these challenges.
Nanoparticles have higher thermal conductivity than base fluids [13] and thus are commonly used to boost the heat transfer phenomenon. This is a new kind of heat transfer fluid produced by dispersing non-metal or metal nanosize particles with less than 100 nm in size combined with a base fluid such as water, motor oil, and other coolants such as metallic, non-metallic, and carbon nanotubes (CNTs), all of which have reduced heat conduction [14]. Because of its heat conductivity, it also boosts the pressure gradient and pressure distribution [15]. Ferrofluids are also a form of nanofluids that work with magnetically field-dependent compounds [16]. Nanofluids have the potential to be employed in the production of vehicles, aircraft, microreactors, and other products, as shown by recent technical advances [17]. Nano and ferrofluids have revolutionized various industrial applications due to their ability to increase heat transmission [18][19][20]. As a result, many researchers have been attracted to studying nano and ferrofluids for their thermal properties in addition to their non-Newtonian behavior. To do so, Lajvardi et al. [21] performed an experiment to look at the ferrofluid flow within a heated copper tube. They used different magnetic field configurations and noticed a noticeable improvement in the ferrofluid's capacity to transfer heat. The next year, Hojjat et al. [22] conducted an experimental investigation of the rheological properties of the base fluid and three different types of non-Newtonian nanofluids with varied nanoparticle concentrations. After that, Nadooshan et al. [23] investigated a hybrid Fe 3 O 4 -MWCNTs/ethylene-based nanofluid for the change of volume fractions to examine flow behavior that was non-Newtonian. In order to analyze the non-Newtonian fluid flow in a cavity containing two cylinders, Barnoon et al. [24] created a two-phase mixing model.
Numerous studies on heat transmission and fluid flow in lid-driven cavities have been studied by different authors [25,26]. In addition, numerous research studies that investigate non-Newtonian Bingham liquid motion within various geometries of cavities or enclosures have been carried out both experimentally and numerically. A chronological account of these works was provided in the review study of Mitsoulis [27]. Bercovier and Engelman [28] made the first attempt at numerical modeling of the Bingham fluid in 1980. Using the finite element approach, they obtained a two-dimensional reference solution for lid-driven square cavity flow. Mitsoulis et al. published several works [29,30] that consider Bingham fluids. In the area of numerical analysis of Bingham fluid flow, Mitsoulis and Zisis' work [31] is remarkable since it employed the finite element approach with two separate grids; each had 400 and 1600 elements, and the numerical solutions were displayed graphically. The findings are very useful in determining the location and size of the yielded and unyielded zones for various Bingham number values. Another outstanding work was carried out by Vola et al. [32] employing the Galerkin technique to examine the unsteady laminar cavity flow of Bingham fluid in the condition of creeping motion, that is, Re = 0. Neofytou [33] used the third-order upwind finite-volume approach to conducting a numerical simulation of cavity (square) flow. The working substance, broad Newtonian fluid, was the focus of the experiment. In this study, various mathematical models, including the power-law, Bingham, Casson, and Quemada models, were taken into consideration when employing the multiple parameters of the corresponding mathematical models displayed graphically.
Several further studies of Bingham fluid flow take into account limited volumes in the literature. For instance, Syrakos et al. [34] used the Papanastasiou regularization strategy to study the viscoplastic (Bingham) flows and provide the numerical findings in both tabular and graphical form. They also looked at how well the finite-volume technique handled Bingham fluid flow in a square cavity driven by a lid. The same authors expanded on their study from a slow-moving river Bingham fluid, taking into account the same geometry and boundary conditions, to the investigation of inertia effects [35]. Pressure-driven flows of Bingham polymers through a square hole were studied by Mitsoulis et al. [31]. In a numerical analysis using the smoothed particle technique (SPH), Rafiee [36] represented the non-Newtonian characteristics of the fluid in a lid-driven cavity by using the generalized Newtonian fluid. Santos et al. [37] used a stabilized finite element approximation to numerically study the flow field's fundamental physics while considering the impact of inertia and rheology factors on the flow of viscoplastic fluids inside a lid-driven cavity. Through numerical simulation, Mahmood et al. [38] examined the cavity flow of Bingham plastic while considering single and double-lid-driven cavities. Due to the physical model's straightforward geometrical design, the Cartesian mesh was used throughout these studies. However, real-world issues have a more intricate geometry than a simple square cavity.
Bingham-type fluid mathematical models are used to investigate the behavior of materials such as paints, pastes, foams, suspensions, cements, and oils [39]. While the Bingham model is considered, two numerical approaches help study the flow of viscoplastic fluids. First, a suitable augmented Lagrangian functional is minimized; second, a variational inequality is resolved; for a description, see Huilgol et al. [40]. Sanchez et al. [41] used a first-order operator splitting technique to resolve the associated variational inequality for the flow of a Bingham fluid while analyzing lid-driven chamber flows. Dean and Glowinski [42] then talked about numerically simulating unsteady flows of a Bingham fluid based on temporal discretization and used it to the flow in a lid-driven cavity to illustrate their point. The operational splitting approach was recently used by Huilgol and Kefayati [43] to replicate spontaneous convection in a cavity filled with a Bingham fluid. A quick scan of the literature in the field reveals that the vast majority of these investigations focus on what are known as "interior flows", such as those in pipes with circular and noncircular cross-sections and slits, laminar and chaotic flow regimes [44], batch mixing containers [45,46], in porous media [47,48], etc. The fluid mechanical characteristics of these processes have also been researched far more thoroughly than their equivalent heat and mass transport. In contrast, a large portion of the literature on flows of the external surface layer type refers to a sphere [49,50] or two-dimensional geometries such as square [51], circular [52], or elliptical cylinders.
In the mixed convective flow, a wavy chamber can significantly affect fluid movement and heat transmission. One of a wavy cavity's main impacts is that it can create secondary movements, for example, vortices and recirculation zones, which can speed up the rate of heat transmission inside the cavity. These secondary flows usally are the result of buoyancy forces and can be particularly pronounced in mixed convection, where both natural and forced convection is present. A wavy cavity can improve heat transmission while also serving as a useful tool in a number of applications. For instance, wavy cavities have been used in heat transfers and cooling systems for electronics as well as microfluidic devices for mixing or dividing fluids. Based on the aforementioned literature, it is evident that further research in the wavy cavity using the Bingham model's effect on the mixed convection combined with the temperature-dependent non-Newtonian nanofluid is required. In this regard, the present study aims to analyze the effects of wavy cavities using finite volume simulations of MHD mixed convection and entropy production of temperature-dependent Bingham nanofluid flow. Body-fitted non-orthogonal grids are considered to derive governing equations in an irregular shape domain. For intricate or irregular geometry, the lid-driven square wavy cavity flow is inadequate. Benchmark results using non-orthogonal grids are likewise limited for comparing the numerical solutions for this kind of geometry. The body-fitted non-orthogonal grids and the Bingham model with the Papanastasiou regularization technique are employed. The correlation equations are determined using the Response Surface Methodology (RSM) approach using the input values for the average Nusselt number, total entropy production, and Bejan number. Plotting the response surfaces and contours generated from the RSM shows the interaction impacts of the relevant factors on the heat transfer rate. The output response's sensitivity to the input parameters is also examined. Given the relevance and variety of applications of the wavy cavity, as well as the influence of Bingham fluid on mixed convection, a sensitivity of the average Nusselt number to the effective parameters has been performed.

Physical Model Description
The schematic diagram with the coordinate system and mesh composition is shown in Figure 1, where Figure 1a shows the model geometry with coordinate system, and Figure 1b shows the mesh distribution of the system. A laminar fluid flow in a chamber is characterized as incompressible with the Boussinesq approximation. A horizontal magnetic field is applied to the flow. The Bingham non-Newtonian Al 2 O 3 -water nanofluid is assumed to be in the enclosure. Furthermore, it is assumed that the base fluid and solid nanoparticles are flowing at thermal equilibrium and at the same velocity. The vertically heated left wavy wall contributes to the heat transmission process in the square chamber. In Figure 1a, it is shown that the top (driven wall) wall maintains an initial velocity of u = u o . When it comes to the heated and cold walls, there is a difference between the left vertical wall's T H and the right wall's T C (T H > T C ). According to the Figure 1b, the mesh composition is densely distributed near the walls and sparsely distributed in the middle of the enclosure. The enclosure's wavy left and right walls can be expressed as follows: where a = 0.05 is the fixed amplitude of the waviness, and the number of periods N = 6.

Bingham Model
A Bingham fluid is a viscoplastic substance that acts like a solid or rigid body under low stress but flows like a viscous fluid under high stress [53]. When the shear stresses approach a threshold level, a large number of non-Newtonian fluids tend to slide across solid surfaces [54]. The Bingham model, which exhibits a linear stress-to-rate-of-strain ratio during flow, is the most basic mathematical model for a viscoplastic fluid. As a result, the constitutive equation for the Bingham model is [55]: where τ is the extra stress tensor,γ is the rate-of-strain tensor, and µ o is the plastic viscosity. There are numerous approaches for removing the Bingham model's discontinuity. Among these, the Papanastasiou [12] regularization technique is a widely accepted method. By using the Papanastasiou [12] regularization technique, the constitutive Equation (2) can be modified as below: where the constant viscosity is denoted as µ f , yield stress is denoted by τ y , and stress growth parameter is denoted by m.

Physical Properties of Non-Newtonian Bingham Nanofluid
The nanofluid's effective density (ρ n f ) can be computed using the mixture rule stated by [56,57]: where ρ f is the base fluid's density, ρ s is the solid particle's density, and φ is the nanoparticle volume fraction. The heat capacity (ρC p ) n f and the thermal expansion coefficient (ρβ) n f of the nanofluid are determined using the mass averaging technique for nanofluid [57] are as follows: where (ρC p ) f and (ρβ) f are the heat capacity and thermal expansion coefficient of the base fluid, respectively, and (ρC p ) s and (ρβ) s represent the heat capacity and thermal expansion coefficient of the solid particle, respectively. The electrical conductivity for the nanofluid is defined by σ n f from the Maxwell-Garnetts (MG) models [58,59] for nanofluid as: where σ f is the electrical conductivity for the base fluid, and σ s is the electrical conductivity for the solid particle. The extra-stress tensorτ is expressed in the dimensional form in terms of the shear rateγ as follows:τ =μ(|γ|)γ, (8) In this study, the Papanastasiou [12] regularization technique has been used to construct the stress-deformation behavior of the fluids with yield stress, which is defined by (dimensional):τ where constant viscosity is denoted as µ f , yield stress is denoted by τ y , and stress growth parameter is denoted by m. From the model of Corcione et al. [60,61], the nanofluid's viscosity (µ n f ) is defined as: In this study, the nanoparticle size d solid = 25 nm is considered, and the d f luid (base fluid's diameter) is given by: where B = 18.01528 × 10 −3 Kgmol −1 is denoted by the base fluid's molecular weight, and N A = 6.022 × 10 23 mol −1 is denoted by the Avogadro number. The temperaturedependent effective thermal conductivity is denoted by K n f , which is shown below, as defined in [60,61]: where Re = After combining Equations (9) and (10), we obtain the dimensional apparent viscosity as:μ where |γ| is the shear-rate's magnitude: The thermophysical properties for the solid ferroparticles and base fluid have been used from [62][63][64].

Governing Dimensional Equations
Under the Boussinesq approximation, the problem can be modeled as a system of coupled non-linear partial differential equations as below [57]: where the conservation of mass is enforced by (15), the conservation of the linear momentum is enforced by (16) and (17), and the conservation of energy is enforced by (18), together with the initial conditions and boundary conditions where g is the gravitational acceleration. B o is the applied magnetic field, and T is the temperature.

Non-Dimensional Governing Equations
Introducing the following non-dimensional variables and parameters, into the dimensional governing Equations (15)-(18), we present the following dimensionless form [65,66]: u-Momentum Equation: Energy Equation: The thermophysical properties for the solid nanoparticles and base fluid have been used from [57,67,68] where the non-dimensional form of apparent viscosity is given by where Bn = τ y H µ f u 0 is the Bingham number, M = mu o H = 1000 is the growth factor of dimensionless stress and |γ| is the dimensionless strain tensor rate, which is defined as: The initial condition and the dimensionless boundary can be expressed below,

Rate of Heat Transfer
The average rate of heat transmission along the wavy hot wall can be determined using the average Nusselt number (Nu) [57,65,69]: where Nu(y) = − k n f k f (T x ) x=0 , Q is the wavy wall's area.

Entropy Generation
Fluid friction, magnetic field, and heat transfer all worked together to create irreversibility in this study. Thus, the total entropy production is defined as the sum of the entropy produced by heat gradients, viscosity, and the magnetic field can be represented as follows [13,57,70]:Ē where entropy produced by fluid friction (Ē F ), heat transfer (Ē T ) and magnetic field (Ē M ) are computed as follows [13]: From Equation (13), putting the value of µ n f we can write, The Bejan number (Be) is the ratio of entropy produced by heat transfer to total entropy production, and it can be defined as follows [13,57,71]: The dimensionless equations can be written as follows after non-dimensionalizing: where D Bn is defined as Equation (31), and The non-dimensional local Bejan number is computed by The non-dimensional entropy production is numerically integrated throughout the cavity volume to determine the overall non-dimensional irreversibilities. The following is how the integral equations are stated:

Numerical Methodology
The governing equations are transformed into the Cartesian coordinates to handle the complex wavy enclosure to a simple square computational domain. The final equations are discretized by the finite volume method (FVM) with a collocated grid system. Grid generation is the crucial step of any numerical method such as FVM. It is a discrete implementation of the computational domain, and the problem is solved based on this domain [55]. This creates limited control volumes (CV) within the field. Each CV for the 2D problem represents a cell with four faces. Figure 2 represents a computational node of a central CV, whose center is at P. E, W, N, and S are the centers of its neighboring CVs. The central CV's four cells are denoted by e, w, n, and s. The variables' values at the cell faces are calculated via interpolation based on the node points. The surface and volume integrals are assessed using the suitable quadrature formula to produce an algebraic equation for each CV. The numerical methodology is described in excellent detail in [72]. In this study, the FORTRAN code is updated to solve the governing equations indicated above [73][74][75]. The code is double precision, entirely implicit, and second-order precise in both space and time. The SIMPLE algorithm has been used in conjunction with a colocated grid arrangement of the variables and a time-marching pressure-correction approach.

Grid Independence Test
The grid independence test (GIT) is an important examination after code validation before the beginning of the numerical simulations for a new study.

Code Validation for the Bingham Fluids
To validate the current program, the acquired results are compared to earlier results for Bingham fluid flow inside the square lid-driven chamber using several values of Bn and Re. Figure 3 shows the u-velocity comparison of the current study with Syrakos et al. [35]. The outcome of this comparison is depicted in Figure 3. The consistency with the results of Syrakos et al. [35] is good in the instance of Re = 1000 and Bn = 10.

Results and Discussion
The finite volume approach has been applied to study the MHD mixed convective flow and entropy production in a wavy enclosure. The non-Newtonian Bingham model has been used to investigate the flow phenomena depending on different parameters. Local and average Nusselt numbers, velocity and temperature distributions, streamline and isotherms have been analyzed in this study, including entropy production with the effect of various pertinent parameters, including the Reynolds number (Re = 10, 100, 200, and 400), Bingham number (Bn = 0, 2, 5, 10, 50, 100, and 200), Hartmann number (Ha = 0, 10, and 20), volume fraction (φ = 0, 0.02, and 0.04), Prandtl number (Pr = 6.2) and Grashof number (Gr = 100). The Richardson number Ri has been computed using Re and Gr. Ri is the controlling parameter for the mixed convection phenomenon. In this study, two different flow characteristics have been considered: one is dominating forced convection with the values of Ri = 0.000625, 0.0025, 0.01, and one is dominating mixed convection by considering the value of Ri = 1. The buoyancy effect causes natural convection, and the lid's motion induces forced convection inside the chamber. The numerical results are discussed in the following sections.
The Bingham number (Bn) is the yield stress ratio to viscous stress. Yield stress usually slows down the flow. Figure 4 shows the dimensionless axial velocity for flow with different Bn and Ha at Re = 10, 100, 200, and 400. The presence of wavy walls and a temperature differential excite the flow domain and induce movement inside the cavity. The flow direction is clockwise and always complies with the produced buoyancy force. At Re = 400, the high-velocity component is owing to the delivered inertia from the surface to the flow, which causes large axial velocity. When an object experiences a force, the item accelerates, which in turn has an impact on its velocity. The inertia force and velocity have a proportionate connection. An increased Re results in increased fluid momentum depletion and lower resistance near the fluid path. The impact of forced fluid movement in layers away from hot and cold sources is mitigated by lowering the fluid penetrability. A high Re also contributes to greater flow velocities within the cavity. Circulation in a cavity with a higher Re and a lower Ri near the walls produces stronger velocity parameters; however, as you approach the cavity's center, the fluid velocity has become negligible caused by a change in the liquid direction, and the same thing happens in the upper half of the cavity. Adding Ha decreases the u-velocity in the cavity, because increases in Ha improve the Lorenz force, resulting in a reduction in velocity distributions. As the Ha rises, the amplitude of the highest horizontal velocity decreases.

Temperature Distributions for various Bn, Re, Ha and φ
The impact of different parameters such as Re, Bn, φ and Ha on temperature distribution (θ) is shown in Figures 5 and 6. Because yield stress may alter the fluid's flow characteristics and heat transmission capabilities, the Bingham number can have an impact on temperature in fluid systems. Particularly, the yield stress can alter how the fluid behaves in areas of high and low shear, which can affect how quickly heat is transferred there. From Figure 5, it is seen that temperature distribution is linear at lower Re. In this case, the convection process is lower. That is why θ decreases linearly. As Re increases, θ changes from linear to chaotic. When Re increases, the inertia force activates and accelerates the flow. Due to the inertia's higher force, the temperature distribution looks chaotic at higher Re. At Re = 400, while the domination of forced convection is maximum (lowest Ri). However, as Ha increases, the temperature increases before x = 0. However, it shifts its trend from x = 0, and θ decreased as Ha increased. Ha induces the magnetic effect on the flow. The magnetic effect opposes the fluid flow; hence, θ decreases as Ha increases.  Moreover, Figure 6 shows that as Bn grows, the temperature distribution becomes linear due to its high viscosity. That means increasing Bn increases the viscosity of the fluid, which opposes the fluid flow, and that is why the change in θ becomes linear as Bn grows. As a result, the overall temperature of the system decreases as Bn increases. From Figure 6a-f, it can be seen that θ is higher for the higher values of Re and φ. The increment happens due to the inertia force and nanoparticles effect. However, overall, θ also increases as φ increases. Adding 4% of nanoparticles accelerates the temperature distribution process.

Effects of Ha, Re, and φ on Streamlines at a Fixed Bn
The effects of Re, Bn, and φ on streamline are presented in Figure 7. At Re = 10, one central elliptical vortex appears on the upper portion of the cavity due to the effect of the mixed convection flow. As Ha increases, the flow trend changes due to the magnetic effect, and at Ha = 20, two secondary vortices appear in the lower portion of the cavity. Due to the weak buoyancy forces, the flow pattern for Re = 100 shows one large vortex on the upper side of the chamber and a secondary vortex at the bottom. Pressure increases in the top left corner while decreasing along the upper right corner. The movement from left to right toward the top is what causes the fluid to recirculate clockwise toward the left. The created top vortex does not completely fill the area due to the mentioned pressure disparity. The shear force is the main cause of the observed secondary lower vortex. At Re = 200 and Re = 400, the central vortex seemed to be shifted to the top right corner of the enclosure, and streamlines are more intensively distributed close to the top wall, which means the area where the fluid is moving faster. The coarse dispersion of streamlines away from the top wall shows that the fluid is moving in a slower area. It is due to weak fluid circulation at the center of the cavity and a high-velocity gradient close to the driven wall.   Figure 8 illustrates the streamlines for various Re, Ha, and φ to provide a more thorough picture of flow and heat transmission in a wavy square cavity. A magnetic field has been applied parallel to the x-axis. When the magnetic field is zero (Ha = 0), a clockwise revolving vortex with its center toward the right upper of the chamber has been discovered. The influence of mixed convection on the fluid flow and the associated heat transfer at various Re is discussed here. The streamlines of a fluid flow are affected by the Bingham number, especially in areas with high shear stresses. When a threshold stress is achieved in a fluid flow with a high Bn number, the fluid first behaves like a solid or exhibits strong resistance to deformation before starting to flow. The fluid's velocity field becomes complicated as a result of this behavior. According to the data, the flow phenomenon consists of one primary vortex on the central top of the cavity for Re = 10. When Re and Ha are increased, a secondary vortex appears at the bottom right corner. A further increment in Ha results in two independent cells that separate the cavity into two different chambers due to the magnetic field and buoyancy effect. The additional chamber found is a direct result of the shear force. At Ha = 20, secondary vortices with opposing, rotations may be seen toward the left bottom end of the chamber. The pattern is the same regardless of fluid type for various Re except Re = 10 at Ha = 20. The number of these cells is insignificant compared to the significant circulation cells. Forced convection is attempted to be balanced by a strong magnetic field. The stream function's value drops as the Ha increases. The weak convective flow in the cavity results from a greater magnetic field, which in turn causes a higher Lorentz force inside the fluid domain when the Ha values are increased.

Effects of Ha, Re, and φ on Isotherms at a Fixed Bn
The impacts of Bn, Re, and φ on isotherms are presented in Figure 9. In many flow field scenarios, Re aids in the prediction of flow patterns. When convection significantly contributes to heat transport, the isotherms become more curved as Re rises. In the scenario of Bn = 100 and Bn = 200 at Re = 10, it can be seen that the isotherms for Bingham fluids are parallel to the vertical walls, but the isotherms for Newtonian fluid flow (Bn = 0) are curved. Due to larger viscous effects in the Bingham fluid at the same nominal value of Re, this difference effectively indicates that the Re at which convection effects are felt is greater for Bingham fluids than for Newtonian fluids. Flows with low Re are likely to be laminar, whereas those with high Re are likely to be chaotic. Chaotic flow is caused by changes in the fluid's direction and speed, which can occasionally cross or even travel in the opposite direction of the flow's main direction. Chaotic flow takes place at high Re (Re = 400) and therefore is dominated by inertial forces, which tend to produce chaotic eddies, vortices, and other flow instabilities; laminar flow occurs at low Re (Re = 10) and is characterized by smooth, constant fluid motion. Temperature graphs show that adding solid nanoparticles to the base fluid increases the nanofluid's thermal conductivity parameter and results in a more even distribution of temperature inside the chamber. The left side wall is where isotherms are concentrated, indicating the location of the quickest fluid movement. Away from the top wall, isotherms are more widely distributed. 5.5. Effects of Ba, Re, and φ on Isotherms at a Fixed Ha Figure 10 shows a substantial response to thermal outlines with varying values of Re. When Re is small, as shown in Figure 10, the isotherms in the Bingham fluid scenario remain almost parallel to the vertical walls, indicating that the heat transport is either diffusion or conduction-driven. The thermal distribution shows that the generated heats accumulated on the left wall while the right wall stayed cold. In this instance, the dominant phenomenon is mixed convection. The fluid flow became disrupted as Re increased. The isotherms become more curved as Re rises. As a result, convection significantly contributes to heat transport. Because of its low density, it begins to compress at the heated left wavy wall. As seen in Figure 10a

Effects of Re and Bn on Unyielded Zone
The effect of Bingham number on unyielded regions is shown in Figure 11 for Bn = 0 to Bn = 50 at Re = 200 and Re = 400 with fixed Re, Ha, and φ. Unyielded zones arise near the lower part of the wavy chamber in the Bingham flow scenarios since stresses are less there because of the distance from the moving wall (top wall) [35]. When Bn increases, these regions enlarge and provide a little area for the flow, driving the vortex upward toward the moving wall. Because of the zero rates of strain situation inside the unyielded zone, the material motion is zero inside the unyielded zones. From the figure, it is seen that the unyielded regions increased as Bn increased.
Furthermore, when Bn rises, the flow area narrows, and the vortex moves toward the motion source (lid). The cavity contains two distinct unyielding zones for all values of the Bingham number, except for Bn = 0 and Bn = 5 (for Re = 200). The smaller portion is at the upper side of the chamber, toward the primary vortex's center or left, while the more significant part is at the lower side of the chamber. Due to the no-slip boundary condition, the velocity is quite low, close to the bottom wall. The upper unyielded portion did not touch any wall. The velocities are non-zero in this chamber portion because of the lid motion. Lid motion allows the fluids to move in the upper part of the chamber even in the presence of higher Bn. As a result, as the Bn increases, the unyielded zone grows but always retains room near the moving wall to form a yielded region.

Effect of Bn and Re on Local Nusselt Number (Nu)
Figure 12a-c represents the impact of different parameters on Nu such as Bn and Re. Nu is the ratio of heat transfer due to convection and heat transfer due to conduction at the boundary layer. Nu is a significant parameter that can help to improve heat transfer rates. It is a non-dimensional parameter equal to the non-dimensional temperature gradient at the surface. It is a measure of the convective heat transfer happening at the surface. Re and Pr mostly determine it. Here, an increased number of Nu was found due to the increased number of Re. That means Nu has a proportionate relationship with Re. As Re increased, the flow density and speed also increased. As a result, the heat transfer also increased. At Re = 10, the mixed convection phenomenon dominates the most; in this case, the heat transfer rate is the least. An increase in the Re value increases the forced convection domination; hence, the heat transfer rate gradually increases. The maximum heat transfer is found for maximum values of Re = 400. Moreover, potential growth is also found in Nu when φ increases from 0% to 4%. On the other hand, the local Nu values obtained for Bingham fluids are smaller than those found for Newtonian fluids (Bn = 0) with the same nominal Re, Ha, and φ values because Bingham fluids have weaker convection due to increased viscous effects. The lowest Nu has been found for the highest Bn. Nu also drops as Ha increases due to the magnetic effect on the flow. Magnetic effects oppose the fluid flow, and the heat transfer rate decreases. Hence, it can be concluded that Nu is proportional with φ and Re and inversely proportional with Bn and Ha.

Effect of Bn, Ha, Re, and φ on Average Nusselt Number (Nu)
The average Nusselt number (Nu) measures the average heat transfer rate of the fluid. In the mixed convective phenomenon based on Bingham fluid, Nu is dependent on different parameters such as Bn, Re, Ha, and φ. The effect on Nu due to Bn, Ha, φ, and Re is shown in Figures 13a-d and 14a-c and Table 2. While the isotherm and the local Nu distribution give precise information about heat exchange, the average Nu is frequently required in engineering applications to size the heat exchange equipment [76]. Furthermore, Nu is required in engineering process computations to estimate the heat transfer rate between the fluid and the heated surface or to assess one of the temperatures if the heat flux is known. Nu is predicted to be a function of Re, Bn, and Ha, explaining why Nu is calculated over a large range of these parameters. It is shown ( Figure 13) that when Bn grows, the obtained value of Nu drops gradually from the value that has been found for the Newtonian case (Bn = 0). This indicates that the convective heat transmission drops when Bn rises, because as Bn increases, the stress of the fluid increases, too. Higher stress retards the fluid flow. Similar results have been found in [76]. However, the Nu values obtained for Bingham fluids are smaller than those found for Newtonian fluids (Bn = 0) with the same nominal Re, Ha, and φ values because Bingham fluids have weaker convection due to increased viscous effects which indeed agreed with the results of the study of Turan et al. [77]. From Figure 14a-c, it is seen that Nu continuously grows as Re increases, but Nu drops as Ha increases. Since Re increases, the buoyancy effect's domination decreases, but the buoyancy effect increases as Ha increases. The buoyancy effect affects the fluid flow that causes a lower heat transfer rate due to the augmentation of Ha and increases the heat transfer rate as Re increases [78]. Potential growth is also found in Nu when φ increases from 0% to 4%. Hence, it can also be stated that Nu is proportional with φ and Re and inversely proportional with Ha. Table 2 shows the effect of those parameters on Nu. Nu is proportional to Re and φ and inversely proportional to Ha and Bn. This is because when Re increases, the fluid's velocity and density also increase. Due to this, the fluid's heat transfer rate increases; hence, Nu also increases for a given number of Ha and Bn. Similarly, if we add an extra 2% and 4% of nanoparticles to the fluid, Nu will increase, too. On the other hand, when Ha increases, the magnetic effect of the fluid also increases. An increasing amount of magnetic field retards the flow phenomenon. As a result, Nu falls for an increased number of Ha. It is also clear from Table 2 that for a larger value of Bn, the convection becomes too weak to affect the thermal transport. That means that the mean heat transfer decreases as Bn increases. It is worth noting that fluid flow can still occur, but this flow is not sufficient to impart any influence on thermal transport.

Entropy Production for Various Re, Ha, and φ
The effects of relevant factors on entropy production are shown in Table 3, where the effects of fluid friction (E F ) t , thermal effect (E T ) t , as well as magnetic field (E M ) t on overall entropy production (E s ) t and Bejan number (Be) are illustrated at Bn = 100 for various Ha, Re, and φ. Entropy production is a measure of energy lost and the low reliability of any system, which refers to the level of irreversibilities present throughout a process. It provides a numerical measurement of irreversibility. The total entropy depends on the summation of other quantities such as (E F ) t , (E T ) t , and (E M ) t . Furthermore, the ratio of (E T ) t to (E S ) t is known as the local Bejan number (Be) [79,80]. When Be > 0.5, the heat transfer irreversibility becomes dominant, and when Be < 0.5, other irreversibilities become dominant. For Be = 1, the irreversibility is caused by heat. Lower entropy production is required to improve the efficiency of energy conversion processes. Overall, reducing entropy generation can lead to more efficient and safe processes, making it an important goal in many fields. From the numerical data, it is observed that increasing Re increases the overall entropy formation (E s ) t . Higher Re induces higher fluid friction and heat transfer to flow. Higher fluid friction and heat transfer cause higher entropy production due to (E F ) t and (E T ) t . Hence, the overall entropy production is attenuated as well. Maximum (E S ) t is recorded at maximum value of Re (Re = 400). According to the obtained result, the overall entropy was attenuated by 110% when Re grows from Re = 10 to Re = 400 at Ha = 0 and φ = 0.04. However, increasing Ha decreases the flow circulation. As a result, the overall entropy reduces in maximum cases with the increment of Ha. At zero magnetic effect Ha = 0, (E M ) t is found 0 for all cases of φ. At Re = 10, the overall entropy reduces 1.01% when Ha increases from 0 to 20. (E s ) t increases when an extra 2% and 4% of nanofluid is introduced. The addition of nanoparticles accelerates the flow phenomenon. Hence, the entropy formation grows as well. The obtained data show that an increment of 1.33% happens when φ increases from 0 to 0.02 at Re = 10. However, the Bejan number (Be) avg shows similarities in the maximum case. It means that (Be) avg does not show any significant difference in its magnitude. The correlation equation for (E s ) t and Be avg can be formulated as below: Be avg = 0.9371 + 0.00015 Re + 0.00057 Ha + 0.00068Re φ − 0.0285Ha φ −1.782 × 10 −7 Re 2 + 0.000049Ha 2 .
These correlation equations are required, since they establish a mathematical link between every factor employed in this research. These formulas are very useful for determining the value of each feasible combination of variables. Figure 15 shows a response surface plot for (E s ) t , which is a visual depiction of the link between total entropy generation (E s ) t and the various parameters that influence it. The x-axis of the graphic indicates the Ha number, the y-axis is the Re number, and the color reflects the entropy generation. The graphic illustrates that increasing the Reynolds number improves entropy generation, demonstrating that this parameter has a beneficial influence on entropy creation. The graphic also demonstrates that the influence of volume fraction on entropy production is positive, which means a higher φ enhances the (E s ) t . Table 3. Effect of various Re, Ha and φ on entropy production for Bn = 100, Pr = 6.2 and Gr = 100.

Re
Ha

Response Surface Methodology of Nu
Response Surface Methodology (RSM) is a statistical approach for modeling and analyzing the interactions between single or multiple responses and input parameters. The primary purpose of RSM is to determine the best possible values of the variables being used to maximize the response variable(s). RSM is employed to examine how different factors affect heat transport in a situation of the average Nusselt number. It analyzes and quantifies the important factors that have an impact on the average Nusselt number. Designing operations to change the input parameters within a predefined range, evaluating the response factor or variables, and then combining a mathematical model with the results are common steps in the process. The response variable(s) may therefore be predicted using the model for any set of input parameters, including their ideal values that optimize the outcome. RSM is one of the helpful approaches for modeling multivariate situations, especially when the interest-generating behaviors are continuously influenced by the input parameters. It has a wide range of practical uses in technology, science, and biological processes, in which it is utilized to optimize methods and goods, decrease costs, and enhance quality. RSM has been used to build improved heat exchangers and improved systems for cooling, increasing the thermal efficiency of different equipment in the field of heat transfer.
The RSM's first goal is to develop a numerical and statistical estimation for the functional connection between the variables [81]. Although there are several RSM models available, a second-order RSM model is widely used because it considers all linear, square, and interaction components for estimating the response. The quadratic polynomial regression model in RSM is a sort of mathematical model that uses a second-order polynomial equation to represent the connection between the input parameters and the output responses, which is considered in the present study. When the influence of the input parameters on the reaction of the output variable(s) is predicted to be non-linear, and the linear model is insufficient to describe the complex nature of the connection, the quadratic polynomial regression model comes into action. The equation of the quadratic polynomial regression model depends on the involved variables. Here, in the RSM process, the number of variables involved as input is four, and there is one output response. The following equation is a representation of the quadratic polynomial regression model having four input variables and one output response.
where 10 ≤ Re ≤ 400, 0 ≤ Ha ≤ 20, 0 ≤ Bn ≤ 100, and 0 ≤ φ ≤ 0.04 are considered in this analysis. Finding the coefficients is the first step in creating a correct correlation between each of these variables (Re, Ha, Bn, and Φ) and the outcome response (Nu). In a series of tests, the outcome variable is calculated for numerous combinations of the input parameters to estimate the model's coefficients. The best-fit coefficient values are then determined by analyzing the data using regression analysis methods. Statistical measurements such as the R 2 value and the adjusted R 2 value are used to assess the model's goodness of fit. In a regression model, the R 2 value shows the amount of variance in a dependent variable that can be determined by the independent variables. A higher R 2 value indicates that the regression model is better fitted to the data. In this case, the measure R 2 value is 99.58% for Nu; this implies that this model is appropriate for determining the response function values. The following model summary table (Table 4) for Nu shows the predicted, adjusted, and actual R 2 values and standard deviation of the analysis: Other than that, an analysis of variance (ANOVA), 2D and 3D surface plot, standard error, and correlation plot are presented to assess the results. The analysis of variance (ANOVA) is a useful statistical method for comparing the means of multiple groups. It is used in studies and analyses of data to examine if variations within categories are statistically significant or just coincidental. ANOVA helps to investigate the impact of several variables on a single dependent variable. For instance, in this case, we examined if there are any significant changes in the Nu due to the effect of the considered parameters. It provided us with a more informed statistical basis for our conclusions. The findings of the statistical analysis of the present mixed convective study using the Bingham model are shown in Table 5. Here, degrees of freedom (DOF) refers to the number of independent elements of data available in the ANOVA test, which generates the mean square values that are used to examine the statistical significance of the differences between groups. The sum of squares (SS) in ANOVA refers to the portion of data differences that identify as a specific source of treatment or the error. It is calculated by the sum of squares for the residual error subtracted from the sum of squares for the model. The p-value and F-value are two essential indicators for determining the statistical significance of an analysis of variance (ANOVA) test's findings [82]. The p-value denotes the chances of detecting a test statistic that is as significant as the one computed from the data. In other words, a low p-value (p ≤ 0.05) indicates that the model is appropriate and significant. On the other hand, a higher F-value indicates that the differences between the groups are substantial. F-values larger than one are important in engineering applications [83]. In statistics, the importance of each term is determined through a comparison of the probability value at the F-value > 1 [84]. Here, the SS value is 112.32 and the F-value is 85.66 which is quite significant, and with the p-value, it concludes that the model for the Nu is statistically significant. The obtained value of the coefficient with the respective p-value is given in Table 6 below: From the above table, the coefficients with the p-value greater than 0.05 (highlighted in bold font in the table) will be omitted from the final regression model equation. The final regression model equation summarized the connection between the input parameters with the output response Nu; after computing the coefficients, Equation (63) can be written as follows:

Model Contour and Response Surface Plot
Based on a mathematical framework developed through response surface analysis, a response surface plot is an illustration of the connection between multiple independent variables and the response variable. The graphic depicts how the output response variable varies when the independent factors' values vary. The independent factors are commonly shown on the x and y axes, while the response variable is depicted on the z-axis in a response surface plot. A model contour plot is a 2D visual depiction of outcome variables as a function of two distinct variables. It behaves similarly to a response surface plot, and rather than a 3D surface, it displays the response variable as a sequence of contour lines in a 2D plot. Both independent variables are commonly shown on the x and y axes, while the response variable is represented by a sequence of contours in a model contour plot.   It is shown that increasing the Bn number causes a fall in the average Nu number because the fluid's yield stress limits its ability to transmit heat. Similar to this, an increase in the Ha number causes a decrease in the Nu number when the magnetic field restricts heat transfer rate and slows fluid motion. On the other hand, the geometry of the system and the boundary conditions determine the precise effects of these parameters on the Nu. In summary, the Bingham and Hartmann numbers have a considerable impact on the Nu in a variety of physical systems. Their individual impacts may be examined and shown using 2D and 3D surface plots as well as computer simulations and experimental tests. As Re increases, the average rate of heat transfer increases as the fluid flow becomes more chaotic, enhancing heat transmission. Furthermore, based on the system, a rise in the volume fraction can have a favorable or detrimental effect on Nu. In this case, increasing the φ improves the thermal conductivity of the fluid, which increases the rate of heat transmission.
The residual plot and the normal probability plot represent two popular diagnostic representations used to examine the hypotheses of a regression model. It represents a visualization that compares the predicted result of the dependent factor with the residuals (the variation among the actual and predicted values of the dependent parameter). Figure 19 depicts the residual plots for Nu, which give an indication of the RSM model's precision level. Figure 19a depicts the normal probability of Nu. The plotted points almost maintain a straight line, indicating that the residuals have a normal distribution. Figure 19b depicts the expected versus actual values for Nu. The residuals, as illustrated, fall close to a straight line, indicating that the regression models are well matched with the actual results.
As a result, we can infer in this situation that there is a substantial correlation between the average Nu number and the independent variable(s) represented on the correlation graph, and the regression line shows the accuracy level of the value of Nu. Figure 20 shows the standard error for Nu at different Re numbers. The fluctuation or uncertainty around the estimated response surface simulation are represented by the standard error in RSM. It shows the degree of variance that could be expected in the projected response at any specific spot inside the design area [85]. The response surface model's standard error is graphically represented by the standard error plot. The contour plot is used in this figure to show the standard error, where lines with equal standard errors are shown on the response surface. The contour lines demonstrate how, as the design space is explored, the standard error varies. The levels of error are very small in different cases of Re, which suggests that the model accuracy is significant and this error can be neglected. In summary, the standard error of the average Nusselt number with the effect of Ha, Bn, and Re shows the level of deviation of the whole process, and it can be determined using statistical techniques. The specific method depends on the assumptions and characteristics of the data.

Sensitivity Analysis
A sensitivity evaluation is a crucial tool in the mixed convective research of Bingham nanofluid because it explains how variations in various variables impact the study's result. Several factors can impact the heat transfer rate in the mixed convective heat transfer of Bingham nanofluid. We can determine the most pertinent variables and understand the influence of all of them on the rate of heat transfer by performing a sensitivity analysis. This aids in optimizing the system's design and enhancing its performance. Furthermore, sensitivity analysis can aid in finding the places where the heat transfer rate is most susceptible to parameter changes. In this study, we individually calculate the partial derivatives of the dependent parameter (Nu) with regard to the independent parameters Re, Ha, Bn, and φ. By setting the values of the independent factors at three levels (low, moderate, and high) (−1, 0, 1), we were able to obtain the desired outcomes.
The regression Equation (64) is used to determine the sensitivity. Sensitivity is measured by calculating the partial derivative of the output function (Nu) with respect to the independent factors [86]. The following is the response function's sensitivity functions: A positive sensitivity output shows that raising the value of the input variable enhances the value of the output function. In addition, a negative sensitivity value reflects a drop in the output function produced by raising the input variable [87]. Figure 21 and Table 7 show the sensitivity analysis for different Bn, φ, Ha, and Re numbers on the output response Nu. It depicts that the sensitivity of Re and φ is positive while the sensitivity of Ha and Bn is negative. Because Ha and Bn have a dampening impact on flow and heat transmission, their sensitivity is negative. When Ha is raised, the electromagnetic forces become stronger and dampen the fluid flow, causing the heat transfer rate to drop. As a consequence, the sensitivity of Ha is negative, implying that a rise in Ha results in a reduction in Nu. The yield stress increases as Bn increases, and the fluid becomes more resistant to flow. Due to the slowed fluid flow, the heat transfer rate decreases. As a result, the sensitivity of the Bingham number is negative, implying that increasing the Bingham number causes a drop in the heat transfer rate. On the other hand, Re and φ have a positive sensitivity on the Nu; because Re and φ have a direct influence on fluid flow and heat transfer, increasing these parameters can lead to an increase in the heat transfer rate. Among these parameters, φ has the highest sensitivity on Nu because the small change in φ (2%, 4%) enhances the heat transfer rate more in comparison with Re.

Conclusions
This study used the Bingham model to investigate the MHD mixed convection and entropy production of non-Newtonian nanofluid in a wavy enclosure numerically. Several non-dimensional factors, including Re, Bn, Ha, as well as φ, are shown to have a significant influence on the results. There are diagrams and tables showing the numerical data to provide visual and numerical comprehension of the findings. The investigation was also performed by the RSM method. The accuracy of the RSM method has been found to be 99.58%. The 2D contour and 3D surface plot presented have been obtained from the RSM. The sensitivity of the Nu to the input parameters also investigated using the RSM method. This analysis is conducted for various values of the aforementioned parameters, and the results are as follows: • For non-Newtonian Bingham fluids, the Hartmann number (Ha), Reynolds number (Re), Bingham number (Bn), and nanoparticle volume fraction (φ) all play an important role in the thermal velocity and temperature as well as the local and average rate of heat transfer in the wavy cavity. • Convective heat transmission in a wavy cavity increases as φ and Re increase and decreases as Bn and Ha increase. • The local rate of heat transfer increases as Re and φ increase. The highest magnitude of local Nu is found at Re = 400. This means that the local heat transfer rate is significant for the lower Richardson number Ri (High Re). • High Re (low Ri) causes greater axial velocity within the wavy cavity. However, the addition of Ha addition decreases u-velocity.

•
The addition of an extra 2% and 4% nanoparticle volume fraction increases the magnitude of ψ max . • At the highest Re and zero magnetic effect, the highest magnitude of Nu has been found when extra 4% nanoparticles are added to the simulation. • Re, Bn and φ have a significant impact on stream function, when Re and φ are progressively raised, ψ max tends to rise as well. However, ψ max is inversely proportional to Ha and Bn.

•
As Re grows, the overall entropy production (E s ) t rises dramatically. When Re grows from 10 to 100, 100 to 200, and from 200 to 400 at Ha = 10 and φ = 0.04, (E s ) t increases by 61.85%, 30.51%, and 31.93% respectively. • The total entropy production (E s ) t is reduced when Ha is increased. (E s ) t decreases by 7.9%, 8.2%, and 8.3%, respectively, as Ha grows from 0 to 20. • At φ = 0, the entropy production caused by the magnetic field is found to be 0. However, (E s ) t increases when an extra 2% and 4% of nanofluid are added to the system. • The local Bejan number(B e ) does not show any significant changes in its magnitude. It shows a little variation with the change of different parameters (Bn, Re, Ha, and φ).

•
The magnetic or Lorentz force reduces convective heat transport; as a result, Nu, θ, (E s ) t , ψ max and velocity are reduced as Ha increases. • The obtained correlation equation from the RSM method shows the relation of output responses to the input parameters.

Conflicts of Interest:
The authors declare no conflict of interest.