High Accuracy Simulation of Electro-Thermal Flow for Non-Newtonian Fluids in BioMEMS Applications

1  Introduction

Numerical methods play a vital role in solving differential equations, which are essential because some physical phenomena can be expressed as differential equations. Ordinary and partial differential equations are the two main types of differential equations. One independent variable makes up ordinary differential equations, while two or more independent variables are involved in partial differential equations. So, ordinary differential equations are one way of discretization, but partial differential equations require more than one way of discretization. Some differential equations exist that contain time-dependent terms. So, these time-dependent terms are approximated by numerical schemes. There exist two types of numerical schemes. Among these, some schemes are explicit, and others are implicit categories. The explicit scheme does not require linearization of the non-linear terms in differential equations, but the implicit schemes require linearization. Also, one additional iterative scheme can be used to solve difference equations discretized by the implicit schemes. However, some implicit schemes are unconditionally stable so that any step length can be used, but primarily explicit schemes are conditionally stable, so they have restrictions on step size. If appropriate step length is chosen for explicit schemes, the stable solution can be obtained; otherwise, the solution will diverge if the step size or involved parameter in the given differential equation does not meet the stability condition. If the finite difference schemes do space discretization, a stability analysis, called Fourier series analysis, exists to find the stability condition. Fourier series analysis can be employed for both explicit and implicit schemes. It gives an exact stability condition for linear differential equations but estimates the exact stability condition for non-linear differential equations. Also, explicit and implicit schemes are divided into multi-step or multi-stage schemes. Multi-stage Runge-Kutta-type schemes require information at a one-time level to find the solution at the next time level. These methods contain one or more predictor stages and one corrector stage.

The process known as electro-osmosis occurs when an electric field causes fluid to flow through a porous substance or across a surface. It happens when a fluid in contact with a charged surface is exposed to an electric field, which causes the fluid to move in reaction to the electric force. The boundary between a solid and a fluid causes the formation of an electric double layer. This layer comprises a Stern layer of adsorbed ions and a diffuse counter-ion layer in the fluid. When an electric field is generated perpendicular to the surface, the charged particles in the fluid move toward the electrodes, dragging the bulk fluid with them. Electro-osmotic flow is the term for the fluid flow produced by this movement. Electro-osmosis is a flexible and effective method for controlling fluid flow in various applications by utilizing the interaction between electric fields and charged surfaces or particles.

The presence of ions in viscoelastic liquids causes the generation of electric Coulomb forces, which exert an electric field on the ions [1]. Electro-osmosis efficiently modifies and regulates blood flow and structure [2]. Levine et al. [3] and the team conducted some of the earliest theoretical investigations in electrokinetic rheology. Misra et al. [4] conducted a theoretical investigation examining how electro-osmotic forces affected the flow of micropolar liquids through a vibrating tiny conduit. Regarding blood rheology, Jubery looked into the physical consequences of electrokinetic events [5]. The effects of an electric field on electro-osmotic flow across a T-junction were investigated by Dutta et al. [6].

Both Jeffrey [7] in 1915 and Hamel [8] in 2009 were the first to find an exact solution to the problem of the constant flow of an incompressible fluid inside two intersecting planes with a converging-diverging character. The two-dimensional flow of Newtonian and non-Newtonian fluids of the Jeffrey-Hamel type via an inclined wall channel has been the subject of numerous investigations since then. The constant two-dimensional Jeffrey-Hamel nanofluid flow in a non-Darcy permeable medium was investigated in [9], along with the thermal leap and changing fluid properties. With an oblique magnetic ground, variable thermal conductivity, heat sink/source influences, and two collateral sheets, Dinarvand et al. [10] computationally investigated the aqueous Fe3O4/CNTs binary nanofluid flow. Scholarly literature such as Garimella et al. [11], Harley et al. [12], etc., extensively details the Jeffery-Hamel flow in convergent-diverging channels, which is fundamental for Newtonian and non-Newtonian fluids. PJ developed the non-Newtonian fluid model. Carreau [13] provides a more comprehensive explanation of the properties of materials whose viscosity is affected by the shear rate. This model effectively depicts the thickening and thinning properties at different shear speeds.

Regarding Carreau nanofluid flow, Khan and Ali Shehzad investigate the effects of microstructure on the oscillating and periodically shifting configuration [14]. Akbar and Nadeem has utilized the Carreau fluid model to represent blood flow via a narrowing, stenotic artery [15]. Recent efforts have explored electro-osmotic transport in Carreau fluids using advanced numerical techniques to capture the non-linear rheological and electrokinetic behaviour under microchannel confinement [16]. The influence of stochastic fluctuations and porous medium resistance on electro-osmotic flow has also been investigated, revealing the critical role of random perturbations and energy dissipation in shaping transport behaviour in microfluidic systems [17]. According to Hayat et al. [18], non-Newtonian fluids have extraordinarily complex constitutive equations, increasing the number of terms and governing equation order. A numerical investigation on the steady incompressible laminar two-dimensional hybrid nanofluid flow upon a convectively warmed moving wedge with radiative transition has been carried out by Berrehal et al. [19] based on the Carreau model of blood viscosity, which is a non-linear model concerning shear rate. For the flow in circular pipes and thin slits, Sochi [20] used two separate fluids, the Carreau, and the Cross fluids, to study the analytical solutions. About Taylor’s famous paint scraping problem, which provides a model for studying wall-driven corner flow caused by an oblique plane moving at a constant speed, Chaffin and Rees [21] studied the behaviour of the inertia-less limit of a Carreau fluid in this kind of system. Recent research has investigated electro-osmotic flow in complicated shapes and viscoelastic fluids. For example, the work by [22] looks at electro-osmotic peristaltic streaming of a fractional second-grade viscoelastic nanofluid with carbon nanotubes in a ciliated tube. This shows how electrokinetics and non-Newtonian behaviour interact in fractional-order modelling frameworks. Using fractional viscoelastic models, researchers have looked at electro-osmotic transport in non-Newtonian fluids [23]. Most of these works are around peristaltic movements and certain shapes.

In contrast, our study looks at the electro-thermal flow of Carreau fluids with full Multiphysics coupling, which is important for BioMEMS applications. Researchers have also investigated magnetized non-Newtonian nanofluids in biomedical settings. For example, the study in [24] did a thermal analysis of blood-based nanofluid flow with a couple of stresses in a vertical microchannel under magnetic effects. This showed how the interaction of the magnetic field and the microstructural fluid behaviour affected temperature regulation, an important factor for thermal control in bio-microdevices. In other studies, optimization methods have been used on complicated non-Newtonian models. For instance, the study in [25] used response surface methodology to optimize the flow of Eyring–Powell fluids with Cattaneo–Christov heat flux and cross-diffusion effects. This shows how advanced thermal models can be combined with transport phenomena to give precise control in non-linear fluid systems.

In this paper, we proposed a modified exponential integrator-based numerical method for solving time-dependent partial differential equations generated in the electro-osmotic flow of Carreau fluid over a stationary plate. Comprising two explicit phases, the suggested method is: the first uses an exponential time integrator to manage stiff linear components; the second uses a Runge-Kutta-type technique to capture the non-linear dynamics. We use a compact finite difference method that can provide fourth- or sixth-order accuracy to improve spatial accuracy. This hybrid computing system performs strongly in stiff, non-linear domains and second-order time precision. Electric potential, Helmholtz-Smoluchowski velocity, and thermodiffusive forces’ influences are included to create the mathematical formulation of the EOF problem. Using the suggested method, the governing equations are non-dimensionalized and solved numerically. The simulation findings show the notable impact of electrokinetic factors on flow behaviour, proving the scheme’s correctness and efficiency in capturing the intricate interactions in non-Newtonian electro-osmotic transport. In this paper, we make the following contributions.

1.    We develop a modified two-stage explicit exponential integrator combining exponential and Runge–Kutta techniques for solving non-linear time-dependent PDEs.

2.    We incorporate a high-order compact finite difference scheme to enhance spatial accuracy, achieving fourth/sixth-order precision.

3.    We model the electro-osmotic flow of Carreau fluid with integrated effects of magnetic field, porous media, heat, and mass transfer.

4.    We captured the influence of variable thermal conductivity and reaction kinetics on flow, temperature, and concentration profiles.

5.    We analyze oscillatory boundary conditions, demonstrating their impact on electrokinetic transport and thermal diffusion. The study provides a detailed parametric investigation (e.g., We,M,Fsme,γ) and their influence on velocity, temperature, and concentration, which has direct implications for the design and tuning of BioMEMS and microfluidic devices where precision control is crucial.

6.    We verify that the system may address complicated fluid behavior without linearization or iterative solvers by maintaining nonlinearity.

The rest of the paper is organized as follows. Section 2 constructs the numerical scheme for solving time-dependent partial differential equations generated in the electro-osmotic flow of Carreau fluid over a stationary plate. Section 3 presents a stability analysis for the proposed scheme. Section 4 presents a problem formulation of the electro-osmotic flow of the Carreau fluid model across a stationary plate. Empirical results are provided in Section 5. Section 6 concludes the paper.

2  Proposed Exponential Time Integrator Scheme

An explicit predictor-corrector scheme is proposed for solving time-dependent partial differential equations. The first stage of the scheme is the predictor stage, while the second stage of the scheme is called the corrector stage. The whole domain is divided into small parts to apply the scheme. First, the solution is found at an arbitrary time level, and then the actual solution will be found at the next time level. For constructing the scheme, consider the following time-dependent partial differential equation.

∂p∂t=F(p,∂p∂x,∂p∂y,∂2p∂y2)(1)

and initial and boundary conditions are given as:

p(0,x,y)=f1(x,y),p(t,0,y)=f2(t,y),p(t,x,0)=f3(t,x),p(t,x,L)=f4(t,x)

where L is a finite number and fi′s are functions of spatial and temporal coordinates.

2.1 Reformulation for Exponential Integrator

For constructing the scheme, Eq. (1) can be written as:

∂p∂t=−2p+P(2)

where P=F+2p. This allows the linear part −2p to be handled by exponential integration.

2.2 Predictor Stage

The predictor or first stage of the scheme can be written as:

p¯l,mn+1=pl,mne−2Δt+(−e−2Δt+1)2Pl,mn(3)

where Δt is temporal step size and Pl,mn=F(pl,mn,∂p∂x|l,mn,∂p∂y|l,mn,∂2p∂y2|l,mn)+2pl,mn. We use an exponential integrator to estimate an intermediate solution p¯l,mn+1 at the next time level. This gives a prediction of the solution using known quantities at time level n.

2.3 Corrector Stage

The second stage, or corrector stage of the scheme, contains three parameters whose values will be computed by matching terms of the Taylor series expansion of the equation given by:

pl,mn+1=apl,mn+bp¯l,mn+1+c(eΔt−1)F(p¯l,mn+1)(4)

The constants a,b,c are determined using Taylor series expansion to ensure second-order accuracy in time.

Rewrite Eqs. (3) and (4) as:

p¯l,mn+1=pl,mne−2Δt+(−e−2Δt+1)2{∂p∂t|l,mn+2pl,mn}=pl,mn+(−e−2Δt+1)2∂p∂t|l,mn(5)

pl,mn+1=apl,mn+bp¯l,mn+1+c(eΔt−1)∂p¯∂t|l,mn+1(6)

Expanding pi,jn+1 using Taylor series expansion as:

pl,mn+1=pl,mn+Δt∂p∂t|l,mn+(Δt)22∂2p∂t2|l,mn+O((Δt)3)(7)

By putting Eqs. (5) and (7) into Eq. (6).

pl,mn+Δt∂p∂t|l,mn+(Δt)22∂2p∂t2|l,mn=apl,mn+b(pl,mn+(−e−2Δt+1)2(∂p∂t|l,mn))+c(eΔt−1){z∂p∂t|l,mn+(−e−2Δt+1)2(∂2p∂t2|l,mn)}(8)

By equating coefficients of pl,mn,∂p∂t|l,mn and ∂2p∂t2|l,mn on both sides of Eq. (8) to yield.

1=a+bΔt=b(−e−2Δt+1)2+c(eΔt−1)(Δt)22=c(eΔt−1)(−e−2Δt+1)2}(9)

By solving a system of linear Eq. (9) the values for a,b and c are:

a=(1−e−2Δt)(1−e−2Δt−2Δt)+2(Δt)2(1−e−2Δt)2b=2Δt(1−e−2Δt−Δt)(1−e−2Δt)2c=(Δt)2(1−e−2Δt)(eΔt−1)}(10)

2.4 Spatial Derivative Approximation

Let F=a1p+a2∂p∂x+a3∂p∂y+a4∂2p∂y2 in Eq. (1), then the first and second stages of the proposed scheme can be expressed as:

p¯l,mn+1=pl,mn+(−e−2Δt+12){a1pl,mn+a2∂p∂x|l,mn+a3∂p∂y|l,mn+a4∂2p∂y2|l,mn}(11)

pl,mn+1=apl,mn+bp¯l,mn+1+c(eΔt−1){a1p¯l,mn+1+a2∂p¯∂x|l,mn+1+a3∂p¯∂y|l,mn+1+a4∂2p¯∂y2|l,mn+1}(12)

So far in this work, a time discretizing scheme is constructed, and now a space discretizing scheme is applied to Eq. (1). We use high-order compact finite difference schemes for spatial derivatives: First derivative in x: matrix form via M1−1N1, First derivative in y: matrix form via M2−1N2 and Second derivative in y: matrix form via M3−1N3.

2.5 Final Predictor-Corrector Form (Matrix-Based)

To implement the space discretizing scheme, matrices are provided as:

p¯l,mn+1=pl,mn+(−e−2Δt+12)[a1pl,mn+a2M1−1N1pl,mn+a3M2−1N2pl,mn+a4M3−1N3pl,mn](13)

pl,mn+1=apl,mn+bp¯l,mn+1+c(eΔt−1)[a1p¯l,mn+1+a2M1−1N1p¯l,mn+1+a3M2−1N2p¯l,mn+1+a4M3−1N3p¯l,mn+1](14)

where Mi′s and Ni′s are matrices generated by the coefficients of the left- and right-hand sides of the following equations.

β1p′|l−1,mn+p′|l,mn+β1p′|l+1,mn=c∘(pl+1,mn−pl−1,mn)2Δx+c1(pl+2,mn−pl−2,mn)4Δx(15)

β1p′|l,m−1n+p′|l,mn+β1p′|l,m+1n=c∘(pl,m+1n−pl,m−1n)2Δx+c1(pl,m+2n−pl,m−2n)4Δx(16)

β2p″|l,m−1n+p″|l,mn+β2p″|l,m+1n=c2(pl,m+1n−2pl,mn+pl,m−1n)(Δx)2+c3(pl,m+2n−2pl,mn+pl,m−2n)4(Δy)2(17)

where c∘=2/3(β1+2),c1=1/3(4β1−1),c2=4/3(1−β2),c3=1/3(10β2−1).

The proposed methodology is an explicit two-stage predictor-corrector method designed to solve time-dependent partial differential equations with high accuracy and computational efficiency. While the second stage improves the answer using a Runge-Kutta-type corrector, the first stage uses an exponential integrator to manage the stiff linear component of the governing equations efficiently. The method obtains second-order time precision using analytical determination of the weighting coefficients via Taylor series expansion. Compact finite difference approximations are used to provide high-order spatial accuracy; they can provide fourth- or sixth-order accuracy depending on the formulation. The method is also designed in matrix form for quick execution and scaling to more dimensional issues. Its clear character makes it especially appropriate for simulating non-linear and electrokinetically driven non-Newtonian flows, such as those modelled by the Carreau fluid under electro-osmotic circumstances with heat and mass transfer influences, since it removes the need for iterative solvers.

3  Stability Analysis

The Fourier series analysis serves as a criterion for determining the stability conditions of finite difference schemes. The study provides precise conditions for linear partial differential equations and assesses the stability conditions for non-linear differential equations. To employ this stability analysis, consider the following transformations.

M1elIψ1+mIψ2=β1e(l+1)Iψ1+mIψ2+elIψ1+mIψ2+β1e(l−1)Iψ1+mIψ2(18)

N1elIψ1+mIψ2=c∘(e(l+1)Iψ1+mIψ2−e(l−1)Iψ1+mIψ2)2Δx+c1(e(l+2)Iψ1+mIψ2−e(l−2)Iψ1+mIψ2)4Δx(19)

M2elIψ1+mIψ2=β1elIψ1+(m+1)Iψ2+elIψ1+mIψ2+β1elIψ1+(m−1)Iψ2(20)

N2elIψ1+mIψ2=c∘(elIψ1+(m+1)Iψ2−elIψ1+(m−1)Iψ2)2Δy+c1(elIψ1+(m+2)Iψ2−elIψ1+(m−2)Iψ2)4Δy(21)

M3elIψ1+mIψ2=β2elIψ1+(m+1)Iψ2+elIψ1+mIψ2+β2elIψ1+(m−1)Iψ2(22)

N3elIψ1+mIψ2=c2(elIψ1+(m+1)Iψ2−2elIψ1+mIψ2+elIψ1+(m−1)Iψ2)(Δy)2+c3(elIψ1+(m+2)Iψ2−2elIψ1+mIψ2+elIψ1+(m−2)Iψ2)4(Δy)2(23)

By employing transformations (18)–(23) into the scheme’s first stage and simplifying it.

p¯l,mn+1=pl,mne−2Δt+(−e−2Δt+12){a1+a2(c∘Isinψ1+2c1Isinψ12Δx(2β1cosψ1+1))+a3(c∘Isinψ2+2c1Isinψ22Δy(2β1cosψ2+1))+a4(c2(cosφ2−1)+2c3(cosφ2−1)(Δy)2(2β2cosψ2+1))+2}pl,mn(24)

Eq. (24) can be written as:

p¯l,mn+1=(γ1+Iγ2)pl,mn(25)

where γ1=e−2Δt+(−e−2Δt+12){a4(c2(cosφ2−1)+2c3(cosφ2−1)(Δy)2(2γ2cosφ2+1))+a1+2} and β2=(−e−2Δt+12){a2(2c∘sinφ1+c1sinφ12Δx(2γ1cosφ1+1))+a3(2c∘sinφ2+c1sinφ22Δy(2γ1cosφ2+1))}.

Now employing the transformations (18)–(23) into second stage of scheme that results in:

pl,mn+1=apl,mn+bp¯l,mn+1+c(eΔt−1)[a1+a2(c∘Isinψ1+2c1Isinψ12Δx(2β1cosφ1+1))+a3(c∘Isinψ2+2c1Isinψ22Δy(2β1cosψ2+1))+a4(c2(cosψ2−1)+2c3(cosψ2−1)(Δy)2(2β2cosψ2+1))]p¯l,mn+1(26)

By rewriting Eq. (26) as:

pl,mn+1=γ3pl,mn+(γ4+Iγ5)p¯l,mn+1(27)

where γ3=a

γ5=c(eΔt−1){a2(c∘sinψ1+2c1sinψ12Δx(2β1cosψ1+1))+a3(c∘sinψ2+2c1sinψ22Δy(2β1cosψ2+1))}

γ4=b+c(eΔt−1){a4(c2(cosψ2−1)+2c3(cosψ2−1)(Δy)2(2β2cosψ2+1))+a1}

By inserting Eq. (25) into Eq. (27) results in:

pl,mn+1=γ3pl,mn+(γ4+Iγ5)(γ1+Iγ2)pl,mn(28)

Eq. (28) can be rewritten as:

pl,mn+1=(γ6+Iγ7)pl,mn

where γ6=γ3+γ4γ1−γ5γ2,γ7=γ1γ5+γ2γ4.

The amplification factor for this case can be written as:

|pl,mn+1pl,mn|2≤γ6+γ7≤1(29)

If the scheme meets condition (29), it will remain stable. The condition (29) can be satisfied by choosing temporal and spatial step size values and involved parameters in the given differential equations.

In this work, the proposed scheme solves the convection-diffusion system. To do so, consider the following matrix-vector equation.

∂f∂t=A1∂f∂x+A2∂f∂y+A3∂f∂y2(30)

where f is a vector and Ai′s are matrices.

The spatial components in Eq. (30) are discretized utilizing a compact methodology, while the proposed method discretizes the time-dependent term.

f¯l,mn+1=fl,mne−2Δt+(e−2Δt−12)[A1M1−1N1fl,mn+A2M2−1N2fl,mn+A3M3−1N3fl,mn−2fl,mn](31)

fl,mn+1=afl,mn+bf¯l,mn+1+c(eΔt−1)[A1M1−1N1f¯l,mn+1+A2M2−1N2f¯l,mn+1+A3M3−1N3f¯l,mn+1](32)

Theorem 1: The proposed computational scheme and compact spatial discretization converge for the vector-matrix Eq. (30).

Proof 1: The proof of this theorem begins with the following exact scheme.

F¯l,mn+1=Fl,mne−2Δt+(−e−2Δt+12)[A1M1−1N1Fl,mn+A2M2−1N2Fl,mn+A3M3−1N3Fl,mn+2Fl,mn](33)

Fl,mn+1=aFl,mn+bF¯l,mn+1+(eΔt−1)[A1M1−1N1F¯l,mn+1+A2M2−1N2F¯l,mn+1+A3M3−1N3F¯l,mn+1](34)

After subtracting Eq. (31) from (33) and let fi,jn−Fi,jn=Ei,jn , etc.

F¯l,mn+1=Fl,mne−2Δt+(−e−2Δt+12)[A1M1−1N1El,mn+A2M2−1N2El,mn+A3M3−1N3El,mn+2El,mn](35)

Upon taking ‖⋅‖∞ on both sides of Eq. (35), it is obtained:

E¯n+1=Ene−2Δt+|−e−2Δt+12|[‖A1M1−1N1‖∞En+‖A2M2−1N2‖∞En+‖A3M3−1N3‖∞En+En](36)

Rewrite Eq. (36) as:

E¯n+1=λ1En(37)

where λ1=e−2Δt+|−e−2Δt+12|[‖A1M1−1N1‖∞+‖A2M2−1N2‖∞+‖A3M3−1N3‖∞+2].

Now, subtracting (32) from (34) yields.

El,mn+1=aEl,mn+bE¯l,mn+(eΔt−1)[A1M1−1N1E¯l,mn+1+A2M2−1N2E¯l,mn+1+A3M3−1N3E¯l,mn+1](38)

By applying norm ‖⋅‖∞ on both sides of Eq. (38) is:

En+1=aEn+bE¯n+1+c|eΔt−1|[‖A1M1−1N1‖∞E¯n+1+‖A2M2−1N2‖∞E¯n+1+‖A3M3−1N3‖∞E¯n+1](39)

By using inequality (37) in inequality (39), the resulting inequality can be expressed as:

En+1≤λ2En+R(O((Δt)3,(Δx)6,(Δy)6))(40)

where λ2=a+bλ1+c|eΔt−1|{A1M1−1N1+A2M2−1N2+A3M3−1N3+1}δ1.

By using n = 0 in inequality (40).

E1≤λ2E0+R(O((Δt)3,(Δx)6,(Δy)6))(41)

Since E0=0, the resulting inequality (41) is:

E1≤R(O((Δt)3,(Δx)6,(Δy)6))(42)

In inequality (40), suppose n = 1.

E2≤λ2E1+R(O((Δt)3,(Δx)6,(Δy)6))≤(λ2+1)R(O((Δt)3,(Δx)6,(Δy)6))(43)

If this is continued, then for a finite number of n.

En≤(λ2n−1+…+λ2+1)R(O((Δt)3,(Δx)6,(Δy)6))=(1−λ2n)1−λ2R(O((Δt)3,(Δx)6,(Δy)6))(44)

When the infinite limits apply, then the infinite geometric series ..+λ2n+…+λ2+1 will converge if |λ2|<1.□

4  Problem Formulations

Think about the non-Newtonian, incompressible, laminar, erratic flow across the fixed plate. The y∗-axis is perpendicular to the x∗-axis, which is taken vertically along the plate. The flow is driven by temperature and concentration gradients, causing thermal and solutal buoyancy forces. The plate surface is warmer and more concentrated than the surrounding fluid. Assume that the concentration and temperature at the plate are higher than the ambient concentration and temperature outside the plate. A transverse magnetic field of strength B∘ is applied perpendicular to the plate. It introduces a Lorentz force, which influences fluid motion through magnetohydrodynamic (MHD) effects. This is relevant in MHD flow control, cooling of microelectronics, and biofluid manipulation. The fluid obeys the Carreau model, capturing shear-thinning behaviour relevant to fluids like blood or polymeric solutions. An electric field is applied, inducing electro-osmotic motion through the electric double layer (EDL) formed near the charged surface. This mechanism is essential in lab-on-a-chip systems and micro-pumps. Porous medium effects (Darcy and Forchheimer resistance) are included, simulating real-world channels or membranes. Fig. 1 illustrates the geometry of the problem. This geometrical setup is physically relevant for BioMEMS, environmental microdevices, micro-reactors, and electrokinetic separation systems, offering a comprehensive simulation framework for complex multiphysics microflows.

Figure 1: Geometry of the Problem

Think about the assumptions about the boundary layer; the equations that regulate the flow phenomena are expressed as [26]:

∂u∗∂t∗=ν(1+(Γ∂u∗∂y∗)2)m−12∂2u∗∂y∗2⏟Non−Newtonianviscousdiffusion+ν(m−1)Γ2∂2u∗∂y∗2(∂u∗∂y∗)2(1+(Γ∂u∗∂y∗)2)n−32⏟Carreaufluidcorrectionterm−(σB∘2ρ+νk)u∗⏟Magnetic+Darcy resistance−b¯u∗2⏟Forchheimerresistance+ρeEx⏟Electro−osmoticbodyforce+gβT(T−T∞)⏟Thermalbuoyancy+gβc(C−C∞)⏟Solutalbuoyancy(45)

∂T∂t∗=1ρcp∂∂y∗(k(T)∂T∂y∗)⏟Heatconductionwithvariablethermalconductivity+q′′′⏟Volumetricheatgeneration(46)

∂C∂t∗=DB∂2C∂y∗2⏟Massdiffusion−k1(C−C∞)⏟Firstorderchemicalreaction(47)

Momentum Eq. (45): This equation is derived from the modified Navier-Stokes equation for a Carreau fluid: ∂u∗∂t∗ represents the time evolution of horizontal velocity, the unsteady inertial term, the first term on the right-hand side is non-Newtonian viscous diffusion, the second term is Carreau-based non-linear viscosity, σB∘2ρ: represents the Lorentz force due to the applied magnetic field, νk is flow resistance through porous medium (Darcy term), b¯u∗2 represents Inertial resistance Forchheimer’s term for non-linear drag, ρeEx: showing the electric body force inducing electro-osmotic motion, gβT(T−T∞): represents buoyancy from temperature gradient and gβc(C−C∞): represents the buoyancy from the concentration gradient.

Energy (Temperature) Eq. (46): ∂T∂t∗: represents the time-dependent temperature change, k(T)=k∞(1+ε1T−T∞Tw−T∞) is variable thermal conductivity and q′′′=kuRρx∗νcp(A∗u(Tw−T∞)+B∗(T−T∞)) represents the internal heat generation due to joule heating (electrokinetic effect) and viscous dissipation or other sources.

Concentration Eq. (47): ∂C∂t∗: represents the time-dependent change in species concentration, DB is mass diffusivity (Brownian diffusion coefficient), ∂2C∂y∗2: represents the diffusive transport in the vertical direction and k1(C−C∞) is the rate of chemical reaction depleting species, assumed first-order (e.g., absorption or degradation). Subject to initial and boundary conditions.

u∗=0,T=T∞,C=C∞ when t∗=0u∗=0,T=T∞+ϵ1(Tw−T∞)cosw∗t∗,C=C∞+ϵ1(C−C∞)cosw∗t∗ when y∗=0u→0,T→T∞,C→C∞ when y∗→∞}(48)

at initial state t=0 fluid is at rest and ambient conditions. At the plate y∗=0, the wall velocity is zero, whereas temperature and concentration oscillate harmonically. When far away from the plate y∗→∞ the fluid properties approach ambient values. This setup simulates time-varying wall conditions in biological tissues or surface-treated membranes exposed to periodic thermal or concentration cycles. The velocity boundary condition at the wall is set to reflect no-slip or oscillatory slip induced by electro-osmotic effects, which is consistent with charged microchannel surfaces exposed to external electric fields. The imposed time-dependent slip velocity mimics scenarios such as alternating-current electro-osmosis (ACEO) or pulsatile pumping mechanisms in lab-on-a-chip devices. The thermal boundary condition assumes a temperature at the wall higher than the ambient temperature, often encountered in thermally regulated microfluidic devices, thermal cycling chambers (e.g., for PCR applications), or localized Joule heating in electrokinetic setups. Where u∗ and v∗ are a horizontal and vertical component of velocity, respectively, T is the temperature of the fluid, C is concentration, respectively, βT is the coefficient of thermal convection, βc is the coefficient of solutal convection, g is gravity, Γ is a material, fluid parameter, ν denotes kinematic viscosity, σ is electrical conductivity, b¯ is non-Darcian parameter, k1is reaction rate parameter, k denotes permeability constant, Ex is electric field component, ρ is the density of the fluid, ρe is the total ionic density. Let B∘ is the strength of the magnetic field applied transversely to the plate, cp is the specific heat capacity, k∞ represents liquid thermal conductivity and ε1 represents smaller parameters elaborating temperature characteristics for thermal-dependent conductivity.

Consider the following transformations [26] to convert (45)–(49) into dimensionless partial differential equations.

y=y∗LR,u=u∗uR,v∘=v1uR,w=tRw∗,t=t∗tR,ϕ=ϕ∗ζ,θ=T−T∞Tw−T∞,ϕ=C−C∞Cw−C∞}(49)

where uR=(νgβTΔT)1/3,LR=(gβTΔTν2)−1/3 and tR=(gβTΔT)−2/3ν1/3. These non-dimensionalizations simplify the equations and reveal governing parameters like Weissenberg number, magnetic parameter, Darcy number, etc.

Derive Dimensionless Governing Eqs. (50)–(52): Eqs. (45)–(47) can be rewritten as dimensionless governing equations by using transformations (49).