Electro-magneto-kinetic thermo-fluid–structure interactions of viscoelastic electrolytes through soft micro-confinements

Abstract

A coupled electro-magneto-hydrodynamic (EMHD) framework for a viscoelastic electrolyte (of Phan–Thien–Tanner (PTT) fluid rheology) flowing through a compliant micro-confinement with linearly elastic walls is developed. The flow is driven by a combination of an imposed pressure gradient and externally applied electric and magnetic fields. Closed-form perturbation solutions are obtained for the velocity, pressure, wall deformation, and temperature. Fluid–structure interactions (FSI) are examined against four parameters: Debye–Huckel parameter ([small kappa, Greek, macron]), Weissenberg number (Wi), Hartmann number (Ha), and electrical Reynolds number (S). We show that favourable pressure gradients drive wall contractions towards a converging channel, while adverse gradients cause wall expansion to a diverging geometry. Observations also show that [small kappa, Greek, macron] reduces the pressure requirement through electroosmotic pumping; Wi induces shear-thinning that flattens velocity profiles; Ha has a dual effect – assistive at low Ha and resistive at higher Ha due to Lorentz drag; and S is consistently assistive, lowering the required pressure drop and enhancing near-wall transport. Thermal behaviour is characterized using three parameters—the Biot number (Bi), the Peclet number (Pe), and the wall-to-fluid conductivity ratio (k_r): higher Bi improves cooling via wall–environment exchange, larger Pe increases axial thermal advection and raises fluid temperature, and higher k_r facilitates heat removal and limits thermal buildup. Collectively, the insights provide a systematic approach for regulating hydrodynamic resistance and thermal loading in deformable EMHD microsystems, with potential applications in bio-lab-on-chip technologies and microscale thermal management platforms.

---

1. Introduction

The advent of compliant microfluidic systems has introduced a new paradigm in the design and application of lab-scale flow technologies. Unlike traditional rigid substrates, channels fabricated from elastomers such as polydimethylsiloxane (PDMS) offer flexibility, transparency, and biocompatibility, making them indispensable for emerging micro-electro-mechanical systems (MEMS), micro total analysis systems (μTAS), and lab-on-chip (LOC) platforms.^1–3 Such deformable structures enable rapid and cost-effective prototyping while providing mechanical characteristics akin to soft biological tissues,^4,5 thereby advancing biomedical applications such as point-of-care (POC) diagnostics,^6,7 organ-on-chip (OOC) modelling,^8,9 and drug delivery.^10,11 Their adaptability also facilitates integration of active components—such as pumps, valves, and sensors—into a single device, broadening their impact across biotechnology, chemical processing, and environmental monitoring.^12–14

Studies have sought to unravel the fundamental mechanisms of fluid–structure interaction in deformable microchannels. The nonlinear pressure drop–flow rate relation (Δp–Q) is the hallmark of compliant microchannels, arising from the two-way coupling of hydrodynamic stresses and wall deformation. Early foundations were laid in studies of collapsible tubes^15,16 and extended to microscale PDMS channels, where Gervais et al.¹⁷ proposed one of the first analytical models by treating the walls as Hookean solids and linking channel height variation to Young's modulus and Poisson's ratio. While this framework captured essential trends, it required an empirical fitting parameter derived from coupled simulations, thereby limiting predictive accuracy. Subsequent theoretical developments sought more systematic closures. Dendukuri et al.¹⁸ and Panda et al.¹⁹ introduced generalized lubrication-type formulations for squeeze flows in deformable channels, whereas Christov et al.²⁰ derived asymptotic Δp–Q relations in bending-dominated regimes. Wang and Christov,²¹ through reduced-order one-dimensional models informed by three-dimensional fluid–structure computations, highlighted regimes with and without significant inertial effects. Parallel numerical studies, including full three-dimensional simulations, resolved the coupled elasticity–fluid response and confirmed the utility of dimensionally reduced descriptions in practical operating ranges.^22,23 Non-Newtonian fluid models^24,25 incorporating viscoelastic and shear-thinning rheology have revealed additional resistance modulation and wall–flow interplay, enriching the theoretical landscape.

Compliant microchannels fabricated from elastomers or hydrogels suffer from inherently low flow rates, since only modest pressure gradients can be applied before rupture occurs.²⁶ A common strategy to augment throughput is the use of electroosmotic flow (EOF), generated when an axial electric field drives the electrical double layer (EDL) at the channel walls. Early-phase theoretical models, such as those by Mukherjee et al.,²⁷ extended Dendukuri et al.¹⁸ to capture pressure-assisted EOF in deformable channels, while later studies²⁷ showed that variations in zeta potential and material parameters (slip length, ionic concentration, elastic modulus) could enhance load-bearing capacity. EOF has also been exploited for mixing^28,29 and stability analyses^30,31 in compliant channels. However, EOF is fundamentally constrained by the strength of the applied field—higher voltages lead to significant Joule heating and undesirable physicochemical alterations of transported species.³² To overcome these limitations, recent work has focused on electro-magnetohydrodynamic (EMHD) actuation, where transverse electric and magnetic fields produce an axial Lorentz force that enhances flow controllability. In rigid channels, EMHD transport has been successfully modelled for Newtonian and non-Newtonian fluids^33,34 and applied to imbibition dynamics and streaming potential-driven flows.³⁵ Extending such EMHD frameworks to compliant microchannels remains relatively unexplored, presenting opportunities for new analytical models that couple elasticity with multiphysics actuation.

Conventional microfluidic analyses typically rely on the Newtonian assumption, which proves insufficient when dealing with fluids of higher complexity—such as polymeric solutions, colloidal suspensions including gels and emulsions, biological fluids like blood, or even ionic and magnetically responsive liquids.³⁶ To overcome these shortcomings, a range of constitutive models has been proposed to capture both viscous and elastic contributions to stress, thereby offering more reliable predictions of flow behaviour in deformable microsystems. Early approaches included purely inelastic models such as the Ostwald–de Waele power-law^37,38 formulation and the Carreau–Yasuda model,³⁹ while viscoelastic responses were often approximated using simplified frameworks like the Maxwell⁴⁰ or Kelvin–Voigt⁴¹ models. These formulations, however, tend to remain valid only within narrow categories of nonlinear fluids, as they cannot adequately capture complex features such as creep, recovery, and stress relaxation. More recently, molecular network-based theories have gained prominence for their broader applicability to viscoelastic fluids. Of particular relevance is the Phan–Thien–Tanner (PTT) model,^42,43 an outgrowth of the Lodge–Yamamoto network theory,⁴⁴ which introduces tunable parameters into the constitutive relation and thereby provides a flexible means of representing a wide spectrum of non-Newtonian fluid behaviours. Although the PTT model has been widely employed to study viscoelastic transport in rigid geometries, such as microchannels^45,46 and capillaries,^47,48 its application to compliant microsystems has remained comparatively scarce. Even fewer studies have attempted to incorporate electro-magnetohydrodynamic (EMHD) effects into such settings, leaving a clear gap in the literature. This absence is particularly significant, given that coupling viscoelastic rheology with wall deformability and external electro-magnetic forcing introduces a unique class of nonlinear interactions that cannot be captured by models developed for rigid domains.

A critical gap in the current literature lies in the limited understanding of electro-magnetohydrodynamic (EMHD) transport of viscoelastic fluids within compliant microchannels. The interplay of wall deformability, nonlinear rheology, and coupled electro-magnetic forces remains largely unexplored, leaving key questions unanswered about their combined impact on flow resistance, stability, and control. In this work, we develop a semi-analytical framework to capture these coupled dynamics and identify governing parameters that dictate transport behaviour. By doing so, we aim not only to advance the fundamental physics of EMHD-driven viscoelastic flows but also to enable strategies for precision control in microfluidic technologies.

2. Mathematical formulation

2.1 The physical problem

The two-dimensional schematic of the physical problem under investigation is illustrated in Fig. 1. The system consists of a microchannel formed by two parallel plates of initial thickness D and length L, separated by an inter-plate distance of 2H, where H denotes the half-channel height. The reference coordinate system is placed at the channel inlet along the centerline of the flow domain. The plates are assumed to be infinitely wide in the transverse z-direction, thereby neglecting any gradients in that direction. Furthermore, since the channel height is much smaller than its length, the flow is predominantly one-dimensional, characterized by creeping flow conditions and a fully developed velocity profile.

Fig. 1 Schematic of viscoelastic fluid transport through a parallel-plate microchannel with deformable walls. The ionic sPTT fluid forms electrical double layers (EDL) at the polymer–fluid interfaces. The flow is driven by a combined action of an axial pressure gradient, electric fields applied along the x- and z-directions, and a transverse magnetic field along the y-direction. The resulting fluid flow induces wall deformation, which in turn alters the flow, leading to a two-way fluid–structure interaction (FSI).

The physical system under consideration is a two-dimensional microchannel, as depicted in the schematic in Fig. 1. The channel is formed by two parallel plates of length L and initial thickness D, separated by a distance of 2H, where H is the half-height of the channel. The entire system is sandwiched between two immovable, rigid plates separated by a distance 2M, thereby keeping the system size constant. A Cartesian coordinate system is established with its origin at the geometric center of the channel inlet. The model is based on several key assumptions. The plates are assumed to be infinitely wide in the transverse z-direction, thereby neglecting any gradients in that direction. Furthermore, since the channel height is much smaller than its length, the flow is predominantly one-dimensional, characterized by creeping flow conditions and a fully developed velocity profile.

The channel walls are composed of soft, deformable polymers commonly used in microfluidic device fabrication (e.g., polydimethylsiloxane or hydrogels). The plate material is modelled as a linearly elastic solid defined by distinct Lamé constants {λ_L,G}. Under isotropic fluid pressure, the channel walls undergo deformation, although the outer boundaries of the plates are assumed to be fixed and thus non-deformable. The working fluid is an electrically conducting, symmetric monovalent electrolyte. Its rheological behaviour is described as non-Newtonian, governed by the simplified Phan–Thien–Tanner (sPTT) constitutive model. The fluid's physical properties include density ρ, specific heat capacity c_p, dynamic viscosity μ, permittivity ε, relaxation time λ, extensibility parameter ε, non-affine slip parameter a_s, and mobility parameter a_m. When the electrolyte comes into contact with the polymeric walls, surface charges accumulate at the fluid–solid interface, giving rise to electrical double layers (EDLs) and a non-uniform ionic distribution within the channel.

The flow within the system is driven by a combination of external and internal forces: an axial electric field E_x (along the x-direction), a transverse magnetic field B_y (along the y-direction), a secondary electric field E_z (along the z-direction), and an imposed axial pressure gradient. Consequently, the overall flow results from the superposition of pressure-driven transport, electroosmotic motion, and Lorentz force effects. The pressure field induced within the channel walls contributes to their elastic deformation, which in turn modifies the fluid flow, leading to a two-way fluid–structure interaction (FSI).

Thermal effects are also considered: heat is generated primarily through Joule heating and viscous dissipation and is dissipated to the surroundings maintained at an ambient temperature T_a. Heat transfer occurs through both the fluid (thermal conductivity k_F) and the solid walls (thermal conductivity k_S), thereby coupling thermal, electrical, mechanical, and fluidic phenomena in the system.

2.2 Solution strategy and fundamental assumptions

The system under consideration involves strongly coupled electrostatics, hydrodynamics, structural mechanics, and thermal transport. Rather than employing a single, highly complex multi-parameter perturbation expansion for the entire system, we adopt a sequential, modular solution strategy. This approach relies on decoupling the physics where appropriate, followed by a targeted regular perturbation expansion to resolve the fluid's nonlinear rheology.

To ensure mathematical clarity, the following fundamental assumptions and linearizations governing the problem are formalized at the outset:

(i) Geometric scaling and lubrication limit: the microchannel is characterized by a small aspect ratio, α = H/L ≪ 1, having a low Reynolds number flow, Re ≪ 1. Under this strict creeping-flow regime, transverse gradients heavily dominate axial gradients (∂/∂y ≫ ∂/∂x) and the transverse velocity component is very small in comparison with axial velocity (u ≫ v). This asymptotic scaling allows us to treat the flow as locally fully developed at any cross-section.

(ii) Electrostatic weak-field approximation: the total electrostatic potential within the channel is assumed to decouple linearly as Ψ(x,y) = ψ(y) − E_xx. This explicitly assumes that the externally applied axial electric field (E_x) is sufficiently weak compared to the transverse field of the electrical double layer (EDL). Consequently, the applied field acts as a driving force for electroosmotic flow but does not distort the equilibrium Boltzmann distribution of the ions near the charged channel walls. The transverse potential distribution ψ(y), driven by the boundary wall potential ψ_w, represents the equilibrium electrical double layer (EDL). This state exists purely due to the physicochemical charge accumulation at the solid–fluid interfaces. It is fundamentally independent of the externally applied axial electric field (E_x) and the bulk fluid motion.

(iii) Thermal and electrical decoupling in the solid: the deformable polymeric channel walls (e.g., PDMS) are modelled as perfect electrical insulators. Because no electric current penetrates the solid domain, Joule heating within the solid walls is strictly zero. Furthermore, internal heat generation within the solid matrix due to mechanical strain energy is assumed to be negligible compared to the viscous and electromagnetic dissipation occurring in the fluid domain.

(iv) Targeted rheological perturbation expansion: after applying the aforementioned spatial and electrostatic linearizations, the core mathematical complexity remains localized in the simplified Phan–Thien–Tanner (sPTT) constitutive equation, which introduces a cubic nonlinearity. To resolve this, we employ a regular perturbation expansion strictly for the rheological variables. We define the small perturbation parameter as γ = 2εWi² ≪ 1 (where ε is the fluid extensibility parameter and Wi is the Weissenberg number) and employ the regular perturbation technique to expand axial velocity and shear stress only, which partitions the fully coupled nonlinear problem into a hierarchy of linear boundary value problems. The O(γ⁰) term represents the exact hydrodynamic response of an equivalent Newtonian fluid. This mathematical order physically captures the linear superposition of a pressure-driven flow, electroosmotic pumping (via E_x acting on the EDL), and magnetohydrodynamic (MHD) Lorentz forces acting on the bulk fluid. The O(γ¹) term isolates the purely rheological deviations caused by the fluid's complex microstructure. Physically, this order captures the elastic stress distribution and the shear-thinning characteristics governed by the polymer relaxation time (Wi) and network extensibility (quantified by ε).

Emergence of mixed-order interactions:

While the momentum equations are solved via the targeted perturbation in γ = 2εWi², the true complexity of the system—the “mixed” multiphysics—emerges in the fully coupled fields of total velocity (u), pressure (p), wall deformation (h), and temperature (Θ).

Because the microchannel is compliant, the local transverse wall deformation h is intrinsically tied to the local isotropic fluid pressure p(x). In turn, to conserve the prescribed volumetric flow rate (Q) throughout the channel, the total axial pressure gradient must be dynamically adjusted to balance the superimposed total velocity field u = u₀ + 2εWi²u₁ + O(2εWi²)². Therefore, the calculated deformation h is not merely a static mechanical output; it encapsulates the highly coupled interaction between the applied EMHD fields and the viscoelastic stress relaxation.

Similarly, the energy transport equations are driven by viscous dissipation (τ∂u/∂y) and Joule/MHD heating terms (σ_eB_y²u²) that scale nonlinearly with the velocity and stress fields. Consequently, the resulting temperature distributions Θ_F(x,y) and Θ_S(x,y) inherently reflect the fully coupled, mixed-order interactions of the electrical, magnetic, structural, and viscoelastic phenomena acting simultaneously across the domain.

By formalizing these bounds, the sequential derivations for electrostatics, fluid dynamics, wall deformation, and thermal transport can be clearly and rigorously established.

To visually summarize this sequential solution strategy and clarify the multi-domain interactions, a comprehensive flowchart of the coupled electro-magneto-kinetic, fluid–structure interaction (FSI) and thermal frameworks is presented in Fig. 2. This schematic outlines the progression from the foundational electrostatic and hydrodynamic formulations to the emergent wall deformation. It explicitly highlights the two-way FSI feedback loops that ultimately govern the thermal energy distribution within the system.

Fig. 2 Flowchart illustrating the sequential mathematical solution strategy and the coupled interactions between the electromagnetic, hydrodynamic, structural, and thermal domains. The FSI feedback loops demonstrate how emergent wall deformation dynamically recalculates the local flow and potential fields.

2.3 Electrostatics

As discussed in Section 2.1, when the viscoelastic electrolyte comes into contact with the charged polymeric channel walls, electrical double layers (EDLs) are established at the fluid–solid interfaces, resulting in charge accumulation near the walls and non-uniform charge gradient within the channel. The ionic distribution within the electrical double layer is described using the Boltzmann relation, subject to the following assumptions:⁴⁹ (a) the ionic species are maintained in thermal equilibrium with the surrounding solvent, (b) the dominant interactions between the ions themselves and between ions and the charged channel walls arise from Coulombic forces combined with thermal motion, (c) the local ionic concentration adjusts according to the prevailing electrostatic potential, (d) the potential at any location is treated as an averaged mean-field potential resulting from all other charges in the system, (e) ions are modelled as point charges in a dilute ideal solution, (f) excluded-volume effects and finite ion size are neglected, (g) the dielectric permittivity of the electrolyte is assumed spatially uniform and unaffected by local electric fields, and (h) the electrolyte is treated as a continuous medium at the macroscopic scale.

Under these assumptions, the equilibrium ionic distribution follows the Boltzmann law. Denoting n₊ and n₋ as the local cationic and anionic number densities, n₀ as the bulk ionic density in an electrically neutral solution, z as the ionic valency, e as the elementary charge, k_B as the Boltzmann constant, and [Doublestruck T] as the absolute temperature, we have n_± = n₀ exp(∓zeψ/k_B[Doublestruck T]), where ψ is the local electrostatic potential. The corresponding net volumetric charge density in the EDL is given by ρ_q = (n₊ − n₋)ze = −2zen₀ sinh(zeψ/k_B[Doublestruck T]). The electrostatic potential distribution within the channel can be explained using the Poisson equation, ∇²ψ = −ρ_q/ε. Substituting the expression of ρ_q into the Poisson equation and employing the lubrication approximation (valid for long, shallow channels) yields the classical Poisson–Boltzmann equation as d²ψ/dy² = 2zen₀ sinh(zeψ/k_B[Doublestruck T])/ε. For small wall potentials (|ψ_w| ≤ 25 mV), the hyperbolic sine function can be linearized in accordance with the Debye–Huckel linearization as sinh(zeψ/k_B[Doublestruck T]) ≈ zeψ/k_B[Doublestruck T]. Defining the Debye–Huckel parameter (inverse of the Debye length), , the electrostatic potential field can be described using the equation

Eqn (1) can be solved using the following boundary conditions:

At the channel walls (y = h): ψ = ψ_w (2.a)

At the channel centreline (y = 0): dψ/dy = 0 (2.b)

Using reference parameters, y_ref = H, h_ref = H κ_ref = 1/H and ψ =ψ_w, eqn (1) and (2) can be non-dimensionalized to obtain

At the channel walls (ȳ = h̄): [small psi, Greek, macron] = 1 (4.a)

At the channel centreline (ȳ = 0): d[small psi, Greek, macron]/dȳ = 0 (4.b)

Eqn (3) and (4) can now be solved to obtain the closed-form potential distribution across the channel height:

It must be noted that since the walls are deformable, h̄ is not a constant but a function of x̄, which makes the electrostatic potential a function of both x̄ and ȳ.

2.4 PTT constitutive relationship

We now introduce the constitutive formulation for the viscoelastic electrolyte, modelled using the Phan–Thien–Tanner (PTT) framework. The general form of the constitutive relation for a non-affine PTT fluid is expressed as⁵⁰where τ denotes the shear stress tensor, λ is the relaxation time of the viscoelastic fluid, u is the velocity vector, μ is the coefficient of dynamic viscosity and D is the rate of strain tensor. For low shear rates, the stress coefficient function can be approximated by a linear relation,with ε being the extensibility parameter. Under the assumption of a fully developed flow, the constitutive equation reduces to

To render the formulation dimensionless, we introduce reference quantities: the characteristic stress, τ_ref = μU_HS/H, and velocity, u_ref = U_HS where U_HS = −εE_xψ_w/μ is the Helmholtz–Smoluchowski velocity representing the electrokinetic slip velocity induced by the electrostatic force on the charged species. In addition, the Weissenberg number is defined as Wi = λU_HS/H, which quantifies the ratio of elastic to viscous timescales of the fluid in the microsystem.

Using these scales, the constitutive relation in non-dimensional form becomes

2.5 Flow dynamics

To model transport within the microchannel, the governing momentum equations for the viscoelastic electrolyte are employed. Before proceeding, several simplifying considerations are introduced. (i) Under the lubrication approximation, the incompressible continuity equation indicates that the transverse velocity component is negligible compared with the axial velocity. (ii) Because the flow occurs at the microscale, viscous forces dominate over inertia, and with the additional assumption of moderate electric, magnetic, and pressure fields, inertial contributions to the momentum balance can be omitted. (iii) The axial velocity is primarily a function of the y-coordinate, though a weak dependence on the x-direction persists due to the deformable channel boundaries. (iv) A scaling analysis of the transverse momentum equation further shows that ∂p/∂y = 0, implying that pressure within any given cross-section remains uniform and varies only along the axial direction.

The next step is to evaluate the body forces acting on the viscoelastic electrolyte in the presence of externally applied fields. The imposed fields are given by E = E_xê_x − E_zê_z and B = B_yê_y. The electric current density can be expressed as J = σ_e(E + u × B), and the Lorentz force becomes F = ρ_QE + J × B. Incorporating these contributions into the momentum balance along the axial direction yields (refer to Appendix A for details)

To analyse the system more conveniently, we introduce characteristic parameters, x_ref = L, p_ref = λ_L + 2G, and the standard non-dimensional numbers: Hartmann number, , electric Reynolds number, and the elastic-to-viscous force (EVF) number, ξ = H²(λ_L + 2G)/μU_HSL. Using these scales, eqn (10) can be written in non-dimensional form as

Eqn (9) and (11) form a coupled nonlinear system with respect to ū and [small tau, Greek, macron]_xy. Closed-form solutions cannot be obtained directly due to the cubic nonlinearity introduced by the sPTT model. To address this, we employ a regular perturbation approach with the small parameter 2εWi² where the extensibility parameter is fixed at ε = 0.1 and the Weissenberg number is restricted to Wi ≤ 1. This choice ensures that the perturbation parameter remains sufficiently small for a first-order expansion to be accurate. The velocity and shear stress are expanded as power series,

ū = ū₀ + (2εWi²)ū₁ + O(2εWi²)² (12)

[small tau, Greek, macron]_xy = [small tau, Greek, macron]_xy,0 + (2εWi²)[small tau, Greek, macron]_xy,1 + O(2εWi²)² (13)

Substituting these into eqn (9) and (11) and collecting terms of equal order in 2εWi² yield the following governing equations,

The boundary conditions are specified as follows: (i) no slip at the walls (ȳ = h̄): ū = 0 (i.e., ū_0,1 = 0) and (ii) centre-line symmetry (ȳ = 0): ∂ū/∂ȳ = 0 (i.e., ∂ū_0,1/∂ȳ = 0) and [small tau, Greek, macron]_xy = 0 (i.e., [small tau, Greek, macron]_xy,0,1 = 0). Eqn (14)–(17) form linear ordinary differential systems under these constraints, and the resulting expression for velocity profile evaluates as

The dimensionless volumetric flow rate is then obtained as

Since the flow rate is linked to the imposed pressure gradient, eqn (19) can be inverted to evaluate the pressure variation along the channel. The pressure field will subsequently be used to determine the deformation of the compliant channel walls under isotropic loading.

2.6 Channel wall deformation

We now turn to the analysis of structural deformation in the compliant microchannel walls induced by the internal fluid flow. Although the fluid motion is primarily axial, the isotropic nature of pressure results in a net lateral displacement of the channel walls. At any location, the pressure acts equally in all directions, thereby driving deformation predominantly along the transverse y-direction.

To model the solid response, the channel walls are treated as linear, homogeneous, and isotropic elastic materials, characterized by constant Lamé parameters {λ_L,G}. The governing equilibrium equation, in the absence of body forces and under steady conditions, is given by ∇·τ = 0, where the stress–strain relationship for a linear elastic solid takes the classical form τ = λ_L(∇·δ)I + G[(∇δ) + (∇δ)^T]. Here, δ represents the displacement vector, while the two terms correspond to (i) volumetric deformation due to pressure changes and (ii) shear distortion response of the material, respectively. Using slenderness wall scaling arguments (D ≪ L) and thick wall considerations (H ∼ D), it follows that significant deformation occurs primarily in the y-direction, leading to the simplified Navier–Lamé equation for the wall (refer to Appendix A for details):

Eqn (20) is solved subject to rigid outer boundary conditions and normal stress continuity at the solid–fluid interface. The latter condition represents the two-way fluid–structure interaction (FSI) coupling at the fluid–solid interface.

Outer solid boundary (y = M): δ = 0 (21.a)

Solid–fluid interface (y = h): −p = (λ_L + 2G)∂δ/∂y (21.b)

Introducing reference scales, δ_ref = H, M_ref = H, D_ref = H and d_ref = H, the governing equation and boundary conditions reduce to

and

Solving eqn (22) with conditions (23) yields

[small delta, Greek, macron](x̄,ȳ) = p̄(M̄ − ȳ) (24)

Since we are concerned only with the deformation of the inner fluid–solid interface, i.e., at ȳ = h̄, the expression simplifies. Also, since [small delta, Greek, macron] = h̄ − H̄ = D̄ − d̄, the wall inner boundary deformation becomes

Eqn (25) is now substituted into eqn (19) to establish a differential relation between pressure, p̄, and pressure gradient, dp̄/dx̄. The resulting formulation is solved numerically by imposing the inlet pressure as atmospheric [p̄(x̄ = 0) = 0] and then adjusting the outlet pressure to maintain the desired flow. Once the pressure distribution is obtained, the velocity field at each axial location can be determined using eqn (18).

2.7 Thermal analysis

We now turn to the analysis of heat generation and transport in the microfluidic channel. Since the working fluid is electrically conducting and simultaneously subjected to external electric and magnetic fields, three distinct mechanisms of heat generation arise: (i) viscous dissipation – represented by τ_xy(∂u/∂y), accounting for energy loss due to fluid shear, (ii) magnetohydrodynamic (MHD) heating – given by σ_eB_y²u², corresponding to the interaction between moving charged fluid particles and the applied magnetic field, and (iii) Joule heating – expressed as σ_e(E_x² + E_z²) arising from electrical resistance encountered by current conduction within the fluid. The heat produced in the fluid is conducted into the surrounding solid walls, which subsequently dissipate it to the environment maintained at an ambient temperature T_a. Two simplifying assumptions are introduced to facilitate tractable modelling: (a) since the wall deformation has previously been assumed to be small, internal heat generation within the solid walls can be neglected and (b) the channel walls are considered perfect electrical insulators, thereby eliminating Ohmic heating within the solid medium.

Let T_F and T_S denote the temperature fields in the fluid and solid domains, and k_F and k_S their respective thermal conductivities. The governing steady-state energy equations are then formulated as follows:

For the fluid domain,

For the solid domain,

Eqn (26) and (27) are solved subject to boundary conditions in which heat is dissipated from the outer boundary of the solid domain to the surroundings. Continuity of temperature and heat flux is maintained at the solid–fluid interface, and symmetry is assumed along the centre-line.

Outer solid boundary (y = M): −k_S∂T_S/∂y = h_e(T_S − T_a) (28.a)

Solid–fluid interface (y = h): k_S∂T_S/∂y = k_F∂T_F/∂y (28.b)

Solid–fluid interface (y = h): T_S = T_F (28.c)

At the centre-line (y = 0): ∂T_F/∂y = 0 (28.d)

We now attempt to non-dimensionalize eqn (26)–(28). We introduce the dimensionless temperatures as and , where ΔT_r is the reference temperature evaluated using the Joule heating scale J = σE_x²H²/k_FΔT_r. The following non-dimensional parameters are also defined: Peclet number, Pe = ρc_pU_HSH²/k_FL, Brinkman number, Br = μU_HS²/k_FΔT_r, Biot number, Bi = h_eH/k_S, relative thermal conductivity ratio, k_r = k_S/k_F, the electric field strength ratio, α = E_z/E_x, and axial thermal gradient in the fluid, .

The dimensionless energy transport equations can then be written as:

For the fluid domain,

For the solid domain,

with the associated dimensionless boundary conditions:

Before progressing further, it is important to note that both and are functions of both x̄ and ȳ. To facilitate analytical progress, particularly for the fluid temperature field, additional simplifying assumptions are required. We rely on our geometric scaling: the deformed wall is defined as h(x) = H + δ(x). Given the earlier assumption of small wall deformation (δ/H ≈ O(0.1)), combined with the lubrication approximation (H ≪ L), the axial variation of the channel geometry is very small (δ/L ≪ 1). Because the local velocity and internal heat generation fields evolve so gradually along the x-axis, we posit that the variation of fluid temperature in the axial direction can be approximated as linear. Accordingly, the axial thermal gradient is taken to be constant. Based on this reasoning, we adopt the following ansatz for the fluid temperature field,

Substituting eqn (32) into the governing equations and boundary conditions, closed-form solutions for both the solid and fluid temperatures may be derived. Since the complete analytical expressions are lengthy and cumbersome, they are not explicitly reported here. Instead, their functional dependencies are summarized as follows.

For the fluid domain,

where i = 0, 1, 2, 3, and 4 and j = 1, 2, 3, 4, 5, and 6.

For the solid domain,

We have now completed the mathematical formulation which provides a complete framework for describing fluid flow, wall deformation, and thermal transport, which we now utilize to analyse results and discuss their physical implications.

3. Results and discussion

3.1 Justification for using first-order solutions in the perturbation method

In eqn (12) and (13), first-order accurate solutions for the velocity and shear stress distributions have been assumed. To rigorously assess the validity of truncating the regular perturbation expansion at this order, the velocity profile at the channel exit is evaluated up to the second-order, O(2εWi²)². The zeroth-order solution reflects Newtonian behavior, while the inclusion of higher-order terms captures the viscoelastic characteristics of the sPTT fluid. Physical parameters have been adopted from existing literature,^26,33 with ε = 0.1 and Wi = 1 in order to represent extreme conditions and ensure the robustness of the analysis.

To verify that this truncation holds across different physical regimes, the validation was expanded across varying multiphysics parameters: specifically, varying the electrokinetic parameter ([small kappa, Greek, macron]), Hartmann number (Ha), and electrical Reynolds number (S). As shown in Fig. 3, the incorporation of the first-order term introduces essential viscoelastic effects, leading to a noticeable deviation from the zeroth-order Newtonian profile. Crucially, however, the addition of the second-order term results in only a marginal adjustment. The first-order and second-order velocity profiles nearly overlap across all tested parameter sets. Therefore, to maintain a balance between mathematical simplicity and analytical accuracy, first-order expansions for velocity and shear stress are considered highly accurate and strictly sufficient for the present analysis.

Fig. 3 Changes in velocity distribution at the channel exit as higher-order terms are included in the perturbation method. The zeroth-order solution represents Newtonian fluid behaviour, while the inclusion of higher-order terms highlights the rheological effects of the PTT fluid.

3.2 Model validation

Our first step is to compare our results against the existing literature in order to validate our mathematical model, and for that, we adopt a two-step benchmarking approach. First, we verify the fluid rheology by comparing our results with those of Alves et al.,⁴⁵ who studied pressure-driven flow of a PTT fluid in a rigid parallel-plate microchannel without applied fields. To match their setup, we assign a very high EVF number (ξ), effectively making the channel rigid, turn off all external fields, and consider a non-electrolytic fluid to eliminate electrokinetic effects. This yields an excellent match for both Newtonian (Wi = 0) and viscoelastic (Wi = 1) cases, as can be observed from Fig. 4(a). Second, we validate the inclusion of electroosmotic and electromagnetic effects by comparing our results with those of Chakraborty et al.,³³ who analyzed Newtonian fluid flow in a rigid microchannel under combined electroosmotic and Lorentz forces. We again set a high ξ and choose Wi = 0 to replicate Newtonian behavior. Our results show close agreement across all Hartmann numbers tested (Ha = 0, 0.5, and 1), as evident from Fig. 4(b). Together, these comparisons establish strong confidence in our model's ability to capture the complex flow behaviour of viscoelastic PTT fluids in compliant microchannels under external fields.

Fig. 4 Validation of our model against previous literature. (a) Results are compared against those of Alves et al.⁴⁵ for PTT fluid flow in a rigid microchannel. To match conditions, we use a very high elastic modulus, disable applied fields, and consider purely pressure-driven flow to validate fluid rheology. (b) Results are compared against those of Chakraborty et al.³³ for Newtonian fluid flow in a rigid microchannel. To match conditions, we use a very high elastic modulus and set Wi to zero.

3.3 Fluid–structure interaction (FSI) analysis: channel wall deformation and velocity profile distribution

We begin our analysis by examining the coupled fluid–structure interactions in the compliant microchannel, with emphasis on how four governing parameters influence both velocity distribution and wall deformation: (a) Debye–Huckel parameter [small kappa, Greek, macron], (b) Weissenberg number Wi, (c) Hartmann number Ha, and (d) electrical Reynolds number (S). Throughout this study, the flow rate is prescribed as constant. Consequently, variations in electrokinetic and electromagnetic body driving forces are balanced by corresponding adjustments in the pressure gradient, which may even reverse in direction to conserve volumetric flux. Unless otherwise specified, the baseline parameter set is adopted from the existing literature^26,33 as [small kappa, Greek, macron] = 100, Ha = 1, S = 1, ξ = 10 and Wi = 1, with normalized channel height H̄ = 1, wall thickness D̄ = 0.5 and flow rate Q̄ = 1.

3.3.1 Effect of the Debye–Huckel parameter. Fig. 5(a) illustrates the variation of pressure along the channel length for different values of the Debye–Huckel parameter, [small kappa, Greek, macron]. Unlike rigid channels, where the pressure variation is typically linear, the present results display a distinctly non-linear profile, a direct outcome of the compliant nature of the channel walls. In such confinements, part of the pressure energy is consumed in deforming the elastic walls, thereby introducing curvature into the axial pressure distribution. For [small kappa, Greek, macron] = 50, electroosmotic pumping is relatively weak, and a favorable pressure gradient is required to maintain the prescribed volumetric flow rate, Q̄ = 1. Consequently, the pressure decreases monotonically along the channel length. In contrast, when [small kappa, Greek, macron] is increased to 75 and 100, electroosmotic forcing becomes dominant, reducing the need for pressure-driven assistance and eventually driving the fluid against an adverse pressure gradient. This is manifested in the positive slope of the pressure distribution, where the outlet pressure exceeds the inlet pressure, a hallmark of EO-dominated transport in compliant microchannels.

Fig. 5 Variation of (a) distribution of axial pressure and (b) deformation of the top wall at the solid–fluid interface along the channel length, corresponding to three different values of the Debye–Huckel parameter, [small kappa, Greek, macron] (= 50, 75, and 100).

The corresponding wall deformation profiles are presented in Fig. 5(b). Since wall deflection is governed by the local fluid pressure according to eqn (25) as h̄ = (1 + M̄p̄)/(1 + p̄), the deformation trends closely mirror the pressure distribution. For [small kappa, Greek, macron] = 50, the negative pressure along the channel induces a contraction of the channel walls, with the deformation amplitude growing towards the outlet. Conversely, at higher values of [small kappa, Greek, macron] (75 and 100), the positive pressure along the channel drives wall expansion, with the degree of deformation being more pronounced for [small kappa, Greek, macron] = 100. The nearly monotonic nature of the deformation profiles stems from the nonlinear but smoothly varying pressure field, confirming that the structural response of the channel is directly coupled to the electrohydrodynamic forcing.

Next, we observe the velocity distribution within the channel at three different sections: (a) at the channel inlet, (b) at the channel mid-length, and (c) at the channel exit, for three distinct values of the Debye–Huckel parameter ([small kappa, Greek, macron] = 50, 75, and 100), as shown in Fig. 6(a)–(c). Fig. 6(a) shows no noticeable distinction in the flow profile with a change in [small kappa, Greek, macron]. This is well expected because the channel exit is marked by atmospheric pressure with fixed channel dimensions, and therefore there has not yet been a re-adjustment between pressure-gradient, electrokinetic and electromagnetic forces to cause major changes to the flow behaviour. Now, proceeding to Fig. 6(b), we notice a higher centre-line peak velocity for [small kappa, Greek, macron] = 50, in line with expectations based on the converging nature of the channel with increasing x̄.

Fig. 6 Velocity profile evolution across the upper half of the channel height for three values of the Debye–Huckel parameter, [small kappa, Greek, macron] (= 50, 75, and 100), shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

The velocity distributions across the channel height at three streamwise locations, (a) at the channel inlet, (b) at the channel mid-length, and (c) at the channel exit, for three distinct values of the Debye–Huckel parameter ([small kappa, Greek, macron] = 50, 75, and 100) are shown in Fig. 6. At the inlet (x = 0), all three cases collapse onto a nearly identical profile, since the pressure is uniform at p̄(x = 0) = 0 and the wall position is close to the undeformed state. We further notice that the velocity profiles at this section already exhibit a considerably full profile, which is characteristic of electroosmotically driven flows. At mid-length, differences arise for different values of [small kappa, Greek, macron]. At [small kappa, Greek, macron] = 50, the velocity profiles are conventional, with pressure-driven and electroosmotic forces acting in the same direction, producing a forward flow peaking near the centre-line. However, for higher values of [small kappa, Greek, macron] (75 and 100), the growing electroosmotic dominance introduces a striking departure, where we see that near the channel walls, electroosmotic force sustains forward motion, while in the core, the adverse pressure gradient drives a counterflow, signifying velocity reversal at the centre-line region. This plug-like structure, enhanced by electroosmotic effects, becomes even more prominent at the channel exit, where deformation of the compliant walls further modifies the flow domain. The marked crossover zones in Fig. 6(b) and (c) capture these transitions, indicating the points where pressure-driven and electroosmotic contributions exactly balance, causing the velocity to change sign. These velocity features are therefore a direct manifestation of the nonlinear pressure distribution and compliant-wall deformation trends discussed earlier.

3.3.2 Effect of the Weissenberg number. The impact of viscoelasticity, quantified through the Weissenberg number (Wi), on the pressure field is presented in Fig. 7(a). In all cases, the pressure decreases monotonically along the channel, but the variation is nonlinear owing to the compliant nature of the walls discussed earlier. A systematic trend with Wi is observed: higher viscoelasticity leads to a moderate pressure drop. The Newtonian case (Wi = 0) experiences the steepest decline, whereas for Wi = 1, the reduction is considerably smaller. This behaviour is consistent with the constitutive relation [eqn (9)] [small tau, Greek, macron]_xy + (2εWi²)[small tau, Greek, macron]_xy³ = ∂ū/∂ȳ, which indicates that increasing Wi lowers the effective shear stress, reflecting the shear-thinning tendency of the fluid. Consequently, a smaller pressure gradient is sufficient to sustain the flow. Moreover, the elastic stresses inherent to viscoelastic fluids contribute to momentum transport, thereby easing the reliance on pressure forces.

Fig. 7 Variation of (a) distribution of axial pressure and (b) deformation of the top wall at the solid–fluid interface along the channel length, corresponding to three different values of the Weissenberg number Wi (= 0, 0.5, and 1).

The wall deformation trends in Fig. 7(b) closely follow the pressure behaviour, since pressure directly drives wall deflection. For all cases, the deformation increases monotonically along the streamwise direction, with the maximum occurring at the channel exit. However, the extent of deformation diminishes with increasing Wi. The Newtonian fluid produces the largest deflection, while viscoelastic fluids, especially at Wi = 1, exhibit noticeably smaller deformation. This reduction can be attributed to the lower pressure forces required to sustain the flow in viscoelastic cases, as part of the stress is borne by the fluid's elastic component.

The velocity distributions across the channel height at three axial locations are depicted in Fig. 8. At the inlet (x = 0), all three cases nearly collapse onto a single profile, as the pressure field is uniform and the wall remains close to its undeformed state, as already explained in the previous sub-section. However, as we move downstream, deviations become apparent. At mid-length (x = L/2), the Newtonian case exhibits a higher core velocity compared to the viscoelastic cases, while both Wi = 0.5 and Wi = 1 show visibly flattened profiles. This flattening is a hallmark of viscoelasticity, wherein the elastic stresses resist velocity gradients. By the channel exit (x = L), this distinction sharpens further: the Newtonian fluid sustains the most sheared profile, whereas the viscoelastic cases present more uniform velocity distributions across the channel height. The reason is also due to channel wall deformation. For the Newtonian case, the wall deforms more, and there is more converging section, so a higher centre-line velocity is required to maintain the desired flow rate, while for Wi = 1, the channel is less converging, and so the centre-line velocity is low.

Fig. 8 Velocity profile evolution across the upper half of the channel height for three values of the Weissenberg number Wi (= 0, 0.5, and 1) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

The velocity distributions across the channel height at three axial locations are depicted in Fig. 8 for Wi = 0, 0.5, and 1. At the inlet (x = 0), all three cases nearly collapse onto a single profile, as the pressure field is uniform and the wall remains close to its undeformed state, as already explained in the previous sub-section. Moving downstream to x = L/2, the Newtonian case (Wi = 0) displays a higher core velocity, whereas Wi = 0.5 and Wi = 1 produce visibly flatter profiles—an expected signature of viscoelasticity, where elastic stresses oppose shear and redistribute momentum across the gap. By the exit (x = L), these differences intensify: the Newtonian fluid retains the most sheared profile and the largest centre-line speed, while the viscoelastic cases exhibit markedly fuller, more uniform profiles with reduced peaks. This behaviour is tied directly to wall compliance. For Wi = 0, the larger pressure drop drives greater channel contraction, effectively creating a more converging passage; with the volumetric flow rate prescribed, the reduced area forces a higher centre-line velocity. In contrast, higher Wi alleviates the pressure demand, suppresses wall deformation, and preserves a less-converging geometry, so the peak velocity diminishes even when the profile becomes flatter. The cross-over zones in Fig. 8(b) and (c) mark the locations where Newtonian and viscoelastic profiles intersect. Up to the crossover (toward the core), the Newtonian case maintains larger velocities; beyond it (closer to the wall), the flatter viscoelastic profiles slightly exceed the Newtonian values. These intersections quantify the momentum redistribution caused by elastic stresses: they reduce core speeds while elevating near-wall velocities just enough to conserve the prescribed discharge.

3.3.3 Effect of the Hartmann number. The influence of the Hartmann number (Ha), which quantifies the effect of an externally applied magnetic field, on the pressure distribution, wall deformation, and velocity profiles is illustrated in Fig. 9 and 10. In the governing momentum balance [eqn (11)], the magnetic field introduces two competing contributions: an assistive term +HaS and a resistive term proportional to −Ha²ū. The relative dominance of these terms dictates the observed flow behaviour. Starting with the pressure distribution (Fig. 9a), the baseline case (Ha = 0) corresponds to a flow unaffected by magnetic forcing, where the pressure decreases monotonically along the channel length as both pressure-gradient and Lorentz force drive the flow at the required discharge. When the Hartmann number is increased to Ha = 0.5, the assistive magnetic force term dominates, leading to an exhibition of a weakly adverse gradient relative to the baseline case. This reduction in the required pressure drop maintains the volumetric flow rate, Q̄ = 1. However, at a higher field strength (Ha = 1), the resistive Lorentz term outweighs the assistive contribution, demanding a significantly larger pressure drop to sustain the same flow rate.

Fig. 9 Variation of (a) distribution of axial pressure and (b) deformation of the top wall at the solid–fluid interface along the channel length, corresponding to three different values of the Hartmann number Ha (= 0, 0.5, and 1).

Fig. 10 Velocity profile evolution across the upper half of the channel height for three values of the Hartmann number Ha (= 0, 0.5, and 1) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

The corresponding wall deformation (Fig. 9b) reflects these pressure variations. In the absence of a magnetic field (Ha = 0), the wall undergoes a moderate monotonic collapse along the channel, consistent with the pressure-driven loading. When the Hartmann number is increased to Ha = 0.5, the assistive magnetic force introduces a weak adverse pressure gradient, which locally increases the isotropic pressure in the channel as we progress axially, and consequently at the channel exit, the wall even exhibits a slight expansion relative to its undeformed state. By contrast, at Ha = 1, the resistive Lorentz force dominates, leading to a strong favourable pressure gradient in the channel and therefore the channel pressure reduces as we move towards the channel exit, leading to pronounced wall collapse, producing a more convergent channel geometry toward the outlet.

The effect of these coupled pressure and wall dynamics on the velocity distributions is depicted in Fig. 10. At the inlet (x = 0), the velocity profiles for all cases nearly collapse, as the pressure field is uniform and wall deformation is negligible. At the mid-length section (x = L/2), differences begin to emerge. We see that the marked crossover zones identified in Fig. 10(b) capture the transition from magnetic forcing dominated flow (near the walls) to pressure-dominated effects (near the core region). Near the walls, the case Ha = 0.5 exhibits the largest local velocities because the assistive term +HaS outweighs the resistive term −Ha²ū, delivering momentum to the fluid even while the axial pressure trend is weakly adverse; conversely, at Ha = 1, the resistive Lorentz damping reduces near-wall speeds. In the channel core, the ordering is reversed: Ha = 1 shows the highest centerline velocity (followed by the baseline Ha = 0 and then Ha = 0.5), since the strong resistive magnetic forcing at Ha = 1 necessitates a much larger favourable pressure gradient to maintain the prescribed throughput. This steeper favourable gradient, together with the pronounced wall collapse at Ha = 1, accelerates the core flow. At the outlet (x = L), these tendencies intensify (Fig. 10c): for Ha = 0.5, the locally adverse pressure gradient and slight wall expansion amplify near-wall velocities but depress the core, whereas for Ha = 1, the large favorable pressure gradient and severe wall constriction concentrate flow in the core, producing a sharply peaked profile with suppressed near-wall velocities due to magnetic damping. The crossover loci therefore shift and sharpen downstream, marking the spatial transition between magnetically assisted near-wall transport and pressure-driven core acceleration.

3.3.4 Effect of the electrical Reynolds number. The electrical Reynolds number S (proportional to the transverse electric field E_z) modulates the assistive electromagnetic contribution +HaS in the axial momentum balance. As S increases, the net body forcing shifts toward the assistive side, altering the pressure field (Fig. 11a). For low S (S = 1), the pressure decays more strongly along the channel as at low strengths of transverse electric field, a highly favourable pressure gradient is required to support the volumetric discharge. Whereas, for larger values of S, the pressure-gradient progressively weakens that decline and eventually produces a pressure rise toward the outlet. This change in the sign and magnitude of the axial pressure gradient follows directly from the increased electromagnetic assistance that must be balanced by the pressure term under a fixed volumetric flux.

Fig. 11 Variation of (a) distribution of axial pressure and (b) deformation of the top wall at the solid–fluid interface along the channel length, corresponding to three different values of the electrical Reynolds number S (= 1, 2, and 3).

Wall deformation (Fig. 11b) responds directly to these pressure adjustments. Increased S that drives the pressure upward downstream produces slight wall expansion near the exit, while lower S yields the expected contraction. In short, stronger transverse-field forcing tends to inflate the compliant channel downstream by raising the local isotropic pressure.

The velocity field (Fig. 12) reflects the competing effects of electromagnetic assistance, pressure redistribution and changing geometry. At the inlet, all profiles coincide. By mid-length, a crossover appears: larger S elevates near-wall forward velocity (electromagnetic assistance concentrated near boundaries), whereas the core begins to feel the pressure adjustment. At the outlet, the differences are most pronounced – cases with higher S show enhanced near-wall velocities but reduced centerline speeds relative to low S cases, consistent with a downstream pressure rise that suppresses bulk motion while the boundary regions remain driven by the electromagnetic forcing. The marked crossover loci therefore indicate the transverse location where transverse-field-induced forcing and axial pressure effects balance and the local velocity changes ordering.

Fig. 12 Velocity profile evolution across the upper half of the channel height for three values of the electrical Reynolds number S (= 1, 2, and 3) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

3.4 Thermal analysis: temperature profile distribution

In this section, we turn our attention to the thermal response of the microchannel system. As outlined previously, heat generation within the fluid arises from three distinct mechanisms: (i) viscous dissipation, represented by τ_xy(∂u/∂y), (ii) magnetohydrodynamic (MHD) heating, given by σ_eB_y²u², and (iii) Joule heating – expressed as σ_e(E_x² + E_z²). The generated heat is conducted into the compliant solid walls and ultimately dissipated to the surroundings maintained at ambient temperature T_a. To characterize the temperature distribution, three non-dimensional groups are considered: (a) the Biot number of the wall, Bi, which governs wall–environment heat exchange, (b) the Peclet number of the fluid, Pe, which sets the relative strength of advection to conduction in the fluid and (c) the ratio of thermal conductivity of solid and fluid, k_r, which controls conduction across the fluid–solid interface. While velocity and wall deformation naturally influence the temperature field, our focus here is restricted to isolating the role of these thermal parameters, with hydrodynamic effects implicitly accounted for in the computed profiles. Additionally in the following discussions, the flow occurs due to a combination of pressure-gradient and Lorentz force, resulting in a converging profile of the channel axially. Unless stated otherwise, the default parameter set employed in Section 3.4 is: Bi = 1, Br = 0.1, k_r = 2, A = 1, Pe = 1, J = 1, and α = 1. These default parameters remain unaltered unless variations are introduced only when examining the specific influence of a parameter.

3.4.1 Effect of the Biot number (Bi). The effect of the Biot number on the temperature field is presented in Fig. 13 at three axial locations: inlet (x = 0), mid-length (x = L/2) and channel exit (x = L). Two consistent observations can be drawn. First, for a fixed axial location, an increase in the Biot number leads to a reduction in the overall temperature within the channel. Second, for a given Biot number, the temperature rises progressively as the flow advances downstream.

Fig. 13 Temperature profile evolution across the upper half of the channel height in both solid and fluid domains for three values of the Biot number Bi (= 1, 2, and 3) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

The first observation may be understood by recalling that the Biot number represents the ratio of internal thermal resistance within the wall to the external thermal resistance offered by convective exchange with the surrounding environment. A higher Biot number thus signifies stronger heat transfer from the microsystem into the ambient through wall convection. This enhanced dissipation pathway reduces the accumulation of thermal energy in the fluid domain, leading to an overall reduction in channel temperature across all axial locations.

The second observation, an axial rise in temperature, can be ascribed to two complementary mechanisms. The most direct factor is continued volumetric heat generation along the channel length, arising from viscous dissipation, magnetohydrodynamic (MHD) heating due to the Lorentz force acting on charged fluid elements, and Joule heating from current conduction under the applied electric fields. These effects are cumulative, and their persistence ensures a steady increase in local fluid temperature downstream. A secondary contribution arises from axial convection of thermal energy, wherein the velocity field continuously transports generated heat downstream. This convective effect amplifies with distance from the inlet, resulting in higher temperatures near the channel exit compared to the entrance, even for the same Biot number.

The second observation, an axial rise in temperature, can be ascribed to two complementary mechanisms. The most direct factor is continued volumetric heat generation along the channel length, arising from viscous dissipation, magnetohydrodynamic (MHD) heating due to the Lorentz force acting on charged fluid elements, and Joule heating from current conduction under the applied electric fields. These effects are cumulative, and their persistence ensures a steady increase in local fluid temperature downstream. A secondary contribution arises from axial convection of thermal energy, wherein the velocity field continuously transports generated heat downstream. This convective effect amplifies with distance from the inlet, resulting in higher temperatures near the channel exit compared to the entrance, even for the same Biot number.

3.4.2 Effect of the Peclet number (Pe). Fig. 14 presents the effect of the Peclet number (Pe) on the temperature distribution at three axial positions. Two clear trends emerge. First, increasing Pe from 1 to 3 elevates the overall system temperature across the channel. Since the Peclet number represents the ratio of advective to diffusive heat transport, a higher Pe implies stronger axial convection relative to transverse diffusion. As a result, the heat generated within the fluid – arising from viscous dissipation, MHD heating, and Joule heating – is predominantly advected downstream rather than diffused through the walls. This limits heat removal by wall conduction and leads to greater retention of thermal energy within the fluid domain.

Fig. 14 Temperature profile evolution across the upper half of the channel height in both solid and fluid domains for three values of the Peclet number Pe (= 1, 2, and 3) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

Second, for any fixed Peclet number, the system temperature rises steadily from the channel inlet to the exit. This increase reflects the continuous nature of internal heat generation and its downstream accumulation under convective transport. The effect is more pronounced at higher Peclet numbers, where the stronger advective flux carries larger amounts of thermal energy towards the exit. Additionally, the downstream region of the compliant channel experiences enhanced viscous stresses, which further intensify dissipative heating. Together, these mechanisms explain the observed monotonic rise in fluid temperature along the channel length.

3.4.3 Effect of the solid–fluid thermal conductivity ratio (k_r). Fig. 15 illustrates the influence of the solid-to-fluid thermal conductivity ratio (k_r) on the temperature distribution along the channel. Two principal observations can be made. First, an increase in k_r leads to a noticeable reduction in fluid temperature across all axial positions. Since k_r is defined as the ratio of wall to fluid thermal conductivity, higher values correspond to walls with stronger heat conduction capacity. Such walls are more effective in channelling the internally generated heat from the fluid into the solid domain, where it is subsequently dissipated to the surroundings. Conversely, when k_r is lower, the weaker wall conduction restricts heat removal, thereby elevating the overall system temperature. Across all k_r, temperatures still rise from the inlet to the exit due to sustained heat generation and convection, but the rise is noticeably tempered at higher k_r owing to enhanced wall conduction.

Fig. 15 Temperature profile evolution across the upper half of the channel height in both solid and fluid domains for three values of the solid–fluid thermal conductivity ratio k_r (= 1, 2, and 3) shown at three axial positions, (a) x = 0 (inlet), (b) x = L/2 (mid-length), and (d) x = L (exit). Owing to symmetry about the centre-line (y = 0), only the upper half of the profile is displayed.

4. Conclusions

This work developed a coupled electro-magnetohydrodynamic (EMHD) model for a single viscoelastic electrolyte confined within a deformable parallel-plate microchannel. The formulation integrates electroosmotic forcing, magnetic effects, transverse electric field contributions, wall elasticity, and viscoelasticity through the simplified PTT model. Perturbation-based closed-form solutions allow us to systematically characterize the flow behaviour, wall deformation, and conjugate heat transfer.

The fluid–structure interaction nonlinearly alters the pressure distribution, as part of the pressure load deforms the compliant walls—favourable pressure gradients drive wall collapse into a converging channel shape, whereas adverse gradients induce wall expansion and a diverging geometry. Electrokinetic forcing ([small kappa, Greek, macron]) and viscoelastic shear-thinning (Wi) reduce the pressure burden and flatten the velocity profile due to reduced pressure gradient requirement for driving the flow. Magnetic forcing (Ha) reveals a dual role: at low Ha, assistive terms reduce the pressure burden and induce slight wall expansion, while at high Ha, resistive Lorentz forces dominate, generating stronger favorable gradients and enhanced wall collapse. In contrast, the electrical Reynolds number (S) always acts assistively, thereby consistently lowering the required pressure drop and enhancing near-wall transport.

Thermal performance is governed by three internal heat sources—viscous dissipation, Joule heating, and magnetohydrodynamic heating—counteracted by heat conduction through the walls and convective losses to the surroundings. A higher Biot number enhances cooling by strengthening wall–environment heat exchange, whereas a larger Peclet number increases fluid temperature by intensifying axial advection. The wall-to-fluid thermal conductivity ratio further modulates heat removal, with more conductive walls promoting dissipation and reducing thermal build-up within the channel.

In summary, this work establishes how electrical, magnetic, viscoelastic, and elastic mechanisms collectively dictate flow and thermal performance in compliant EMHD microsystems. The results provide clear guidelines for tuning hydrodynamic resistance and thermal management by adjusting field strengths, wall properties, and material conductivity, offering useful strategies for lab-on-chip, biofluidic, and microscale energy devices.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5sm00876j.

References

1. F. Meshkinfam and G. Rizvi, A MEMS-based drug delivery device with integrated microneedle array—design and simulation, J. Biomech. Eng., 2021, 143(8), 081010.
2. R. Keçili and C. M. Hussain, Green micro total analysis systems (GμTAS) for environmental samples, Trends Environ. Anal. Chem., 2021, 31, e00128.
3. Y. Zhao, X. Lv, X. Li, N. Rcheulishvili, Y. Chen, Z. Li and Y. Deng, Microfluidic actuated and controlled systems and application for lab-on-chip in space life science, Space:Sci. Technol., 2023, 3, 0008.
4. E. Song, Y. Huang, N. Huang, Y. Mei, X. Yu and J. A. Rogers, Recent advances in microsystem approaches for mechanical characterization of soft biological tissues, Microsyst. Nanoeng., 2022, 8(1), 77.
5. S. C. Inbody, B. E. Sinquefield, J. P. Lewis and R. E. Horton, Biomimetic microsystems for cardiovascular studies, Am. J. Physiol.: Cell Physiol., 2021, 320(5), C850–C872.
6. D. Liu, Y. Wang, X. Li, M. Li, Q. Wu, Y. Song and C. Yang, Integrated microfluidic devices for in vitro diagnostics at point of care, Aggregate, 2022, 3(5), e184.
7. S. M. Yang, S. Lv, W. Zhang and Y. Cui, Microfluidic point-of-care (POC) devices in early diagnosis: A review of opportunities and challenges, Sensors, 2022, 22(4), 1620.
8. M. M. Keshtiban, M. M. Zand, A. Ebadi and Z. Azizi, PDMS-based porous membrane for medical applications: design, development, and fabrication, Biomed. Mater., 2023, 18(4), 045012.
9. F. Alghannam, M. Alayed, S. Alfihed, M. A. Sakr, D. Almutairi, N. Alshamrani and N. Al Fayez, Recent progress in PDMS-based microfluidics toward integrated organ-on-a-chip biosensors and personalized medicine, Biosensors, 2025, 15(2), 76.
10. A. Shakeri, S. Khan and T. F. Didar, Conventional and emerging strategies for the fabrication and functionalization of PDMS-based microfluidic devices, Lab Chip, 2021, 21(16), 3053–3075.
11. Y. Liu, G. Yang, Y. Hui, S. Ranaweera and C. X. Zhao, Microfluidic nanoparticles for drug delivery, Small, 2022, 18(36), 2106580.
12. K. Ma, S. Wu, Y. Zheng, M. Shao, J. Zhang, J. Wu and J. Zhang, Development and integration of carbon–polydimethylsiloxane sensors for motion sensing in soft pneumatic actuators, Actuators, 2024, 13(8), 285.
13. N. Mustakim, M. Pandey and S. W. Seo, Surface photografted thermoresponsive hydrogel microvalves on PDMS/silicon hybrid membrane for light-actuated localized chemical release, Sens. Actuators, B, 2025, 433, 137559.
14. J. Kim, A. Shanmugasundaram, C. Park and D. W. Lee, A PDMS diaphragm-based bioreactor integrating mechanical and electrical stimulation for enhanced cardiomyocyte functionality and maturation, Sens. Actuators, A, 2025, 387, 116468.
15. T. J. Pedley, The fluid mechanics of large blood vessels, Cambridge University Press, Cambridge, 1980.
16. S. I. Rubinow and J. B. Keller, Flow of a viscous fluid through an elastic tube with applications to blood flow, J. Theor. Biol., 1972, 35(2), 299–313.
17. T. Gervais, J. El-Ali, A. Günther and K. F. Jensen, Flow-induced deformation of shallow microfluidic channels, Lab Chip, 2006, 6(4), 500–507.
18. D. Dendukuri, S. S. Gu, D. C. Pregibon, T. A. Hatton and P. S. Doyle, Stop-flow lithography in a microfluidic device, Lab Chip, 2007, 7(7), 818–828.
19. P. Panda, K. P. Yuet, D. Dendukuri, T. A. Hatton and P. S. Doyle, Temporal response of an initially deflected PDMS channel, New J. Phys., 2009, 11(11), 115001.
20. I. C. Christov, V. Cognet, T. C. Shidhore and H. A. Stone, Flow rate–pressure drop relation for deformable shallow microfluidic channels, J. Fluid Mech., 2018, 841, 267–286.
21. X. Wang and I. C. Christov, Theory of the flow-induced deformation of shallow compliant microchannels with thick walls, Proc. R. Soc. A, 2019, 475(2231), 20190513.
22. O. Ozsun, V. Yakhot and K. L. Ekinci, Non-invasive measurement of the pressure distribution in a deformable micro-channel, J. Fluid Mech., 2013, 734, R1.
23. D. Chakraborty, J. R. Prakash, J. Friend and L. Yeo, Fluid-structure interaction in deformable microchannels, Phys. Fluids, 2012, 24(10), 102002.
24. V. Anand, J. David Jr and I. C. Christov, Non-Newtonian fluid-structure interactions: Static response of a microchannel due to internal flow of a power-law fluid, J. Non-Newtonian Fluid Mech., 2019, 264, 62–72.
25. E. Boyko and I. C. Christov, Non-Newtonian fluid-structure interaction: Flow of a viscoelastic Oldroyd-B fluid in a deformable channel, J. Non-Newtonian Fluid Mech., 2023, 313, 104990.
26. A. Roy and P. Dhar, Fluid–structure-interactive elasto-and thermo-hydrodynamics of electrokinetic binary fluid flows in compliant micro-confinements, Phys. Fluids, 2024, 36(3), 032008.
27. U. Mukherjee, J. Chakraborty and S. Chakraborty, Relaxation characteristics of a compliant microfluidic channel under electroosmotic flow, Soft Matter, 2013, 9(5), 1562–1569.
28. S. Mitra, S. Maity, S. Sutradhar and D. Bandyopadhyay, Electroosmosis with augmented mixing in rigid to flexible microchannels with surface patterns, Ind. Eng. Chem. Res., 2019, 59(9), 3717–3729.
29. Y. Wei, Y. Chen, J. Xu and J. Li, Induced charge electro-osmotic mixing performance of viscoelastic fluids in microchannels with an electrically conductive plate, Phys. Fluids, 2023, 35(8), 083102.
30. X. Wang and I. C. Christov, Reduced modelling and global instability of finite-Reynolds-number flow in compliant rectangular channels, J. Fluid Mech., 2022, 950, A26.
31. R. Patne and V. Shankar, Instability induced by wall deformability in sliding Couette flow, Phys. Fluids, 2020, 32(11), 114102.
32. A. Roy and P. Dhar, Micro-imbibition electro-magneto-hydrodynamics of viscoelastic fluids with interactive streaming potential, J. Non-Newtonian Fluid Mech., 2022, 310, 104936.
33. S. Chakraborty and D. Paul, Microchannel flow control through a combined electromagnetohydrodynamic transport, J. Phys. D: Appl. Phys., 2006, 39(24), 5364.
34. Z. Y. Xie and Y. J. Jian, Rotating electromagnetohydrodynamic flow of power-law fluids through a microparallel channel, Colloids Surf., A, 2017, 529, 334–345.
35. A. Roy and P. Dhar, Capillary orientation and morphology dictated oscillatory electro-magneto-imbibition of viscoelastic electrolytes, Langmuir, 2024, 40(45), 23788–23805.
36. M. T. Shaw and W. J. MacKnight, Introduction to polymer viscoelasticity, John Wiley & Sons, 2018.
37. W. Ostwald, The scientific foundations of analytical chemistry: treated in an elementary manner, Macmillan and Company, 1900, p. 105.
38. R. S. Morrell and A. Waele, Rubber, resins, paints and varnishes, D. Van Nostrand Company, 1920.
39. P. J. Carreau, I. F. MacDonald and R. B. Bird, A nonlinear viscoelastic model for polymer solutions and melts—II, Chem. Eng. Sci., 1968, 23(8), 901–911.
40. J. Clerk Maxwell, On the dynamical theory of gases, Philos. Trans. R. Soc. London, 1867, 157, 49–88.
41. W. Voigt, Ueber innere Reibung fester Körper, insbesondere der Metalle, Ann. Phys., 1892, 283(12), 671–693.
42. N. P. Thien and R. I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech., 1977, 2(4), 353–365.
43. N. Phan-Thien, A nonlinear network viscoelastic model, J. Rheol., 1978, 22(3), 259–283.
44. A. S. Lodge Mathematics Research Centre, U.S. Army Technical Report No. 1092.
45. M. A. Alves, F. T. Pinho and P. J. Oliveira, Study of steady pipe and channel flows of a single-mode Phan-Thien–Tanner fluid, J. Non-Newtonian Fluid Mech., 2001, 101(1–3), 55–76.
46. M. Trivedi, S. Maurya and N. Nirmalkar, Numerical simulations for electro-osmotic flow of PTT fluids in diverging microchannel, Mater. Today: Proc., 2022, 57, 1765–1769.
47. H. S. Gaikwad, A. Roy and P. K. Mondal, Autonomous filling of a viscoelastic fluid in a microfluidic channel: Effect of streaming potential, J. Non-Newtonian Fluid Mech., 2020, 282, 104317.
48. V. N. Phan, C. Yang and N. T. Nguyen, Analysis of capillary filling in nanochannels with electroviscous effects, Microfluid. Nanofluid., 2009, 7(4), 519.
49. J. H. Masliyah and S. Bhattacharjee, Electrokinetic and colloid transport phenomena, John Wiley & Sons, 2006.
50. L. L. Ferrás, A. M. Afonso, M. A. Alves, J. M. Nóbrega and F. T. Pinho, Electro-osmotic and pressure-driven flow of viscoelastic fluids in microchannels: Analytical and semi-analytical solutions, Phys. Fluids, 2016, 28(9), 093102.

---