Optimization and Sensitivity Analysis of Non-Isothermal Carreau Fluid Flow in Roll Coating Systems with Fixed Boundary Constraints: A Comparative Investigation

1  Introduction

The Classical Navier-Stokes equations are insufficient to describe the rheological complex fluids, motivating the development of constitutive models for non-Newtonian behavior. Interest in non-Newtonian fluid mechanics has grown because of extensive applications in medical devices, manufacturing, energy, and process engineering [1]. Examples include sludge, polymers, crude oil, cosmetics and soap solutions, petroleum products, drilling fluids, and coatings. The role of non-Newtonian fluids in enhancing coating performance offers significant benefits for achieving a uniform, thin-coated layer on a surface. The constitutive relationship between stress and deformation rate tensor is shown using several rheological models. There are many fluid models available, including those by Jeffrey, Sisko, Casson, differential types, grade n fluids (e.g., Phan Thien and Tanner), Maxwell, elastic-viscous, Oldroyd’s family, and Carreau [2,3]. The non-linear properties of non-Newtonian fluids present difficulties in deriving solutions by analytical and numerical techniques. The complexity of these fluids requires alterations to the governing formulas. The major governing equations are formulated by adjusting the momentum equation to incorporate the stress tensor obtained from the constitutive equation. To transform the set of equations, the lubrication approximation method (LAT) is applied, and specific boundary conditions have been established in alignment with the specified problem. This article covers numerical solutions through the finite difference method (FDM) [4], the Midrich boundary value problem (BVP) method [5], and analytical solutions using the Perturbation method (PM) [6].

The process of applying a consistent layer of liquid to a moving substrate is called roll coating (RC) [7]. This process is extensively employed to apply accurate and uniform coatings on various substrates, including magnetic electronic media, adhesive tapes, foils, plastic films, printed materials, paper, and paperboard. This procedure enhances the functional and aesthetic properties of products, including gloss, printability, and smoothness, while also improving their durability, performance, and resistance to environmental influences. Fig. 1 illustrates a schematic of the roll-coating setup analyzed in this study. Where a rotating roll moves upward, entraining a non-Newtonian Carreau fluid into the narrow nip region created between the roll and a stationary constrained wall. The confined gap induces significant shear, pressure, and heat-transfer effects that shape the film development. As the roll extends beyond the wall, a uniform coating layer adheres to its surface, influenced by geometric parameters, like roll radius, wall position, and nip gap, governing the velocity, pressure, temperature, and coating uniformity.

Figure 1: Diagram exhibiting the RC configuration specific to the current problem [8]

Significant studies integrating analytical, experimental, and computational methodologies regarding the roll-coating (RC) process have been reported in Greener and Middleman [9] investigated the RC process of a viscoelastic fluid utilizing the conventional LAT and derived the solution via the perturbation approach. They observed a decline in the pressure with a slight increase in coating thickness. Coyle [7] previously conducted a thorough examination of the fluid dynamics of RC, focusing on steady flows, rheology, and stability, and clarifying the critical flow characteristics, stability factors, and rheological influences involved in these processes. Hinter and White [10] examined the water flow occurring between two rolls in their research findings. The lubricating principle was employed to confirm that the results aligned with the experimental data.

Previous investigations have significantly contributed to the understanding of roll-coating flows. Greener and Middleman [11] established analytical and experimental frameworks linking pressure distribution with coating thickness. However, their work omitted to take surface tension, free surfaces, and contact line behaviour into consideration. Owens et al. [12] investigated misting in RC and determined that it’s non-linearly dependent on time for relaxation as a result of filaments and bead breakdown. Additionally, Lécuyer et al. [13] incorporated roll deformation and free-surface dynamics, demonstrating how roll compliance modifies nip hydrodynamics. Ascanio et al. [14] developed validated computational strategies for high-speed and discrete gravure coating of complex fluids. Echendu [15] implemented analytic and computational methodologies to examine the flow attributes of coatings for a variety of fluid types, including Newtonian and non-Newtonian fluids [16]. These summaries of earlier achievements clarify prior progress and highlight the need for the present study, which integrates non-isothermal Carreau fluid behavior with fixed boundary constraints and optimization analysis.

The study of non-Newtonian fluids has become increasingly significant due to their diverse applications across engineering, industrial processes, and commercial manufacturing [17–20]. Early numerical analyses by, to explore the importance of fluid flow characteristics in the reverse roll coating process (RRCP), Jang and Chen [21] conducted this research by considering the volume of fluid-free surface (VOFFS) and finite-volume (FVM) methodologies. Additionally, Shahzad et al. [22] examined the behavior of a non-isothermal couple stress fluid within an RRC configuration, taking into account the slip conditions present at the surfaces of the rolls. Devisetti and Bousfield [23] investigated the impacts of fluid absorption in forward roll coating (FRC), highlighting the significance of rheological qualities in enhancing industrial efficiency, especially for non-Newtonian lubricants. Bhatti et al. [24] utilized advanced numerical techniques such as the generalized differential quadrature and Newton-Raphson methods to precisely solve non-linear flow and heat transfer equations for magnetized nanofluids. Their approach showed significant stability and precision in predicting velocity, temperature, and concentration distributions over rotating surfaces. Recent research by Farooq et al. [25] has explored nanofluid flows over stretching surfaces under magnetic fields and slip conditions, emphasizing the influence of melting heat, stratification, and nanoparticle additives on heat and mass transfer, particularly through machine learning and neural network-based modeling techniques for enhanced predictive accuracy.

Response surface methodology (RSM) is a statistical technique employed to assess the sensitivity of models to the influencing variables. In recent years, it has been utilized by numerous researchers to address diverse physical problems, with extensive procedures and discussions available in [26,27]. The primary aim of RSM is to optimize responses, making it especially valuable for addressing fluid flow problems where the goal is to enhance desired flow characteristics. More broadly, RSM plays a critical role in enhancing product quality, minimizing costs, and promoting innovation across various industries, including manufacturing, chemical engineering, and coating technologies. The effective implementation of this approach emphasizes the enhancement of achievable outcomes. Consequently, statistical analysis is vital for identifying the relative importance of critical parameters and thereby improving model precision. Furthermore, the strategic application of optimization algorithms to enhance goal functions constitutes a key component of modern decision-making processes.

1.1 Research Gap and Objectives

An extensive examination of the literature reveals a significant gap: Although RC processes with various non-Newtonian fluids have been thoroughly studied [8,28–30], the particular scenario of a non-isothermal Carreau fluid under fixed constraining boundary conditions remains unexplored. The primary objective of our study is to examine how a non-isothermal CRF impacts the heat transfer and flow in the RC process. This configuration is highly relevant for industrial applications, including a stationary bar or wall to regulate control coating thickness. The lack of studies on this particular combination limits the optimization and understanding of such processes. Consequently, the primary objective of this work is to address this gap by rigorously examining the influence of a non-isothermal Carreau fluid on heat transfer and flow dynamics within an RC process featuring a fixed constrained wall. Accordingly, this study fills the identified gap by inquiring:

•   How effectively does the non-isothermal Carreau fluid model represent the flow characteristics and coating results with fixed boundary conditions?

•   How does the model predict the velocity and temperature distributions under non-isothermal conditions?

•   How are the pressure distribution and pressure gradient developed within the coating film, and what do they offer about the flow structure and film formation?

•   To what degree can perturbation and numerical methods function as a dependable prediction framework for optimizing roll-coating?

•   In what way does the model quantify engineering factors, including coating thickness, power input, and Nusselt number, and evaluate their dependence on the perturbative and gravity parameter for industrial coating applications?

•   Why and how do fixed walls alter flow resistance, and what consequences does this have for machine wear and product quality?

The primary aims of our study are:

•   Develop analytical (via perturbation approach) and numerical (using BVP Midrich and Finite Difference methods) solutions for the flow and temperature profiles for non-isothermal Carreau fluid.

•   Assess the impacts of key physical parameters on critical engineering quantities, including pressure gradient, coating thickness, separation force, and power input.

•   Conduct optimization and sensitivity analysis with RSM to determine the optimal process parameters.

•   The findings provide new physical insights into how Weissenberg and gravity parameters affect flow behavior, energy consumption, and coating performance, establishing a prediction framework for industrial applications.

Finally, the validity of our results will be confirmed through comprehensive cross-verification between the analytical and numerical methods, as well as comparison with current literature for limited cases. We anticipate excellent agreement, which will validate the robustness of our integrated methodology.

1.2 Motivation for the Work

RC is a crucial process in industries such as packaging, printing, and surface protection, where uniform application of liquid films is essential for product quality and material efficiency. The established model and results are directly applicable to roll-to-roll printing, polymer and packaging film production, textile finishing, pharmaceutical and biomedical coatings, as well as decorative or anti-corrosion surface treatments, where controlling film thickness and minimizing mechanical wear are essential. In practical situations, the procedure becomes more complex when addressing non-Newtonian fluids, whose shear-thinning and viscoelastic characteristics significantly influence coating thickness, heat transfer, and energy requirements. Despite its significance, precise modeling of RC with constrained walls remains challenging due to the non-linear characteristics of governing equations and the existence of constrained boundary conditions. By combining perturbation analysis, numerical verification, and statistical optimization, this work provides a comprehensive framework for investigating RC with non-Newtonian fluids. The findings not only advance theoretical understanding of coating flows under constrained boundary conditions but also offer practical insights for improving coating uniformity, optimizing process efficiency, and minimizing material waste in industrial applications. The novelty of this study lies in formulating a complete non-isothermal Carreau fluid model for roll coating under fixed boundary constraints, an unaddressed configuration in the literature.

2  Geometry and Mathematical Formulation

A schematic of the RC system under investigation is presented in Fig. 1. The setup consists of a rotating roller of radius R positioned against a stationary wall. The roller rotates clock-wise with a constant angular velocity (ω). A steady, smooth stream of an incompressible non-Newtonian fluid coats the roller, and the flow is analyzed under the corresponding non-isothermal conditions and prescribed boundary constraints. The nip separation is the distance between the roll and the wall. As shown in Fig. 2, it is assumed that the flow is moving towards the x-axis and perpendicular along y-axis.

Figure 2: Schematic illustration of the RC process and its variables related to our problem

The coating thickness (λ) is measured at the uppermost point of the roller. Here, g denotes the gravitational acceleration, h(x) is the distance between the roller and the solid surface, and U represents the linear speed of the roller. The shape of the roller surface is described by the geometric function h(x)=1+x22. Because the physical setup exhibits symmetry, the coordinate system shown in Fig. 2 is chosen as the most suitable for the analysis. In this configuration, the y-axis is oriented perpendicular to the fluid layer (pointing outward), while the x-axis is aligned with the direction of the fluid motion. Presented below are the governing mass, momentum, and energy equations expressed in their component form. These formulations correspond to a steady-state flow regime in which body forces are assumed to be absent.

Continuity equation [3]:

∂u¯∂x¯+∂v¯∂y¯=0(1)

Momentum equations [31]:

ρ(v¯∂v¯∂y¯+u¯∂u¯∂x¯)=μ0∂∂x¯{2∂u¯∂x¯[1+Γ2{2(∂v¯∂y¯)2+(∂u¯∂y¯+∂v¯∂x¯)2+2(∂u¯∂x¯)2}]n−12}−∂p∂x¯+μ0∂∂y¯{(∂u¯∂y¯+∂v¯∂x¯)[1+Γ2{2(∂v¯∂y¯)2+(∂u¯∂y¯+∂v¯∂x¯)2+2(∂u¯∂x¯)2}]n−12−F}(2)

ρ(v¯∂v¯∂y¯+u¯∂u¯∂x¯)=μ0∂∂x¯{2∂u¯∂x¯[1+Γ2{2(∂v¯∂y¯)2+(∂u¯∂y¯+∂v¯∂x¯)2+2(∂u¯∂x¯)2}]n−12}−∂p∂y¯+μ0∂∂y¯{(∂u¯∂y¯+∂v¯∂x¯)[1+Γ2{2(∂v¯∂y¯)2+(∂u¯∂y¯+∂v¯∂x¯)2+2(∂u¯∂x¯)2}]n−12−F}(3)

Energy equation [32]:

ρCp(u¯∂T¯∂x¯+v¯∂T¯∂y¯)=(∂2T¯∂x¯2+∂2T¯∂y¯2)+τ¯xx∂u¯∂x¯+τ¯xy(∂u¯∂y¯+∂v¯∂x¯)+τ¯yy∂v¯∂y¯(4)

Extra stress tenor of Carreau Fluid model:

τ¯xy={(∂u¯∂y¯+∂v¯∂x¯)[1+Γ2{2(∂v¯∂y¯)2+(∂u¯∂y¯+∂v¯∂x¯)2+2(∂u¯∂x¯)2}]n−12}(5)

Boundary conditions (BCs) [8]:

•   u(x,h(x))=U imposes the no-slip condition at the moving roll surface y=h(x), where the fluid velocity matches the roll’s linear speed U.

•   u(x,0)=0 enforces the no-slip condition at the stationary wall y=0, where the fluid velocity is zero.

•   T(x,h(x))=T0 sets a constant temperature T0 at the roll surface.

•   T(x,0)=T1 sets a different constant temperature T1 at the wall, establishing the thermal gradient that drives heat transfer.

•   p(xs)=0 and dpdxx=xs are the Swift-Stieber or separation conditions applied at the unknown detachment point x=xs. They state that the pressure must equal the ambient pressure (set to zero) and the pressure gradient must vanish at the flow separation point where the fluid film splits from the wall. These conditions are crucial for determining the location of the separation point and the final coating thickness.

2.1 Rheological Modeling

The current study explores the rheological aspects characterized by the viscoelastic CRF model. The equations defining the CRF model are expressed as follows:

τ¯=−pI+μ(γ)A¯1,(6)

μ=μ∞+(μ0−μ∞)[1+(Γγ⋅)2]n−12,(7)

we consider the most realistic situations μ0≫μ∞. Therefore, μ∞ is regarded as zero; Eq. (5) in its modified form may be articulated as:

τ¯=−pI+μ0[1+(Γγ⋅)2]n−12A¯1(8)

It is essential to emphasize that within the context of the CRF model, a flow behavior index value (n) exceeding 1 is indicative of shear-thinning (pseudoplastic) behavior. Conversely, a value below 1 characterizes shear-thickening (dilatant) behavior. In two-dimensional non-isothermal, incompressible, steady laminar flow, the velocity profile can be articulated in the following form:

V¯=[u¯(x,y),v¯(x,y)](9)

2.2 Non-Dimensional Form via LAT

On the basis of geometrical data from LAT, our simulations suggest that a significant dynamical event happens in the RC towards the nip. Moreover, LAT is used in those flow geometries in which one dimension is markedly smaller in magnitude relative to others. The surface of the roll is almost parallel to the H0<<R for a short distance in any direction from the nip (i.e., in any ± direction). Hence, we can deduce that the flow is nearly parallel in the nip region. The basic assumption of the LAT can be applied. As, u¯>v¯ and ∂∂y¯>∂∂x¯. We quickly conduct an order of magnitude analysis to determine the characteristic scale for the velocity, pressure, and temperature.

The following scales for variation in x¯,y¯andu¯ can be identified as: x¯∼L, y¯∼H0, u¯∼U. After considering the mass conservation Eq. (1) and the relations defined previously, we obtain v¯U∼H0L<<1. This relation shows that the length of the longitudinal characteristic, denoted by L=2RH0 is larger than the transversal velocity v¯ by order of magnitude. Considering the plausible claim that the rolls nip region experienced the most significant dynamical phenomena. The material flows along the x-direction, and no velocity can be detected in the y-direction. The next step is to use the magnitude to generate the following set of non-dimensional attributes [8]:

{u=u¯U,μ¯0=μ0μ,x=x¯(RH0),T¯=T−T0ΔTc,τ¯xy=H0μ0Uτxy,y=y¯H0,v=v¯δU,F=ρgH02μU,p=p¯H0μU(H0R)12.(10)

Applying Eq. (10) to Eqs. (1)–(5) and yields the following equations in dimensionless form.

∂2u∂y2+3(We)(n−12)(∂u∂y)2∂2u∂y2=dpdx+F,(11)

(We)2Pe(u∂T∂x+v∂T∂y)=(We)2∂2T∂x2+∂2T∂y2+PrEc[u∂p∂x+v∂p∂y]+ψ,(12)

ψ=Brμ[2(We)2((∂u∂x)2+(∂v∂y)2)+((We)2∂v∂x+∂v∂y)2],(13)

∂2T∂y2+[1+We.(∂u∂y)2]n−12.Br(∂u∂y)2=0,(14)

τxy=∂u∂y[1+2We{1+(∂u∂y)2}]n−12,(15)

{u=0,y=0,T=0,andu=1,y=h,T=1.(16)

2.3 Validity Limits of the Derived Model

“The present formulation is based on the lubrication approximation, valid when the film thickness is much smaller than the roll radius (H0R) and the flow remains steady, laminar, and predominantly unidirectional. The perturbation expansion assumes a low Weissenberg number (We≪1), where elastic effects are weak, allowing truncation at first order. The temperature field is treated under the assumption of small to moderate thermal gradients and constant fluid properties, so that variations of viscosity and density with temperature are negligible. Additionally, inertial, surface tension, and free-surface deformation effects are ignored, which is justified for thin films at moderate roll speeds. Therefore, the derived model is valid for thin-film, weakly elastic, and moderately non-isothermal Carreau fluid flows. For higher We or large temperature gradients, higher-order and non-linear coupling effects would need to be included”.

3  Methodology

To obtain practical results for the present analysis, firstly, the PDEs are turned into ODEs by applying certain similarity transformations [8]. The governing equations were first solved analytically by applying the perturbation method, which provided approximate expressions for the flow variables. To validate these results, the problem was reformulated numerically using the finite difference method (FDM) along with Newton’s iteration. The associated non-linear algebraic problem was later solved using Midrich BVP-DSolve in Maple to verify the consistency and stability of the derived solutions. Additionally, RSM was utilized to investigate the impact of key factors and to develop prediction models. A sensitivity study was carried out employing regression-based methodology to identify the most influential variables affecting the system performance. The combined application of perturbation theory, numerical methods, and statistical modeling not only enhances the accuracy of the results but also offers a comprehensive framework for examining the system behavior under different parametric conditions. This hybrid approach ensures both mathematical rigor and practical reliability, hence enhancing the overall credibility of the conclusions. The following flow diagram depicts the full study technique, detailing each stage from model formulation to performance evaluation.

3.1 Asymptotic Solution via Perturbation Method

The perturbation technique is an analytical strategy for solving non-linear ordinary differential equations (ODEs). This approach is frequently employed by engineers to address various real-world difficulties. This is done to guarantee that at least one unknown is shown within a series of smaller factors. This research examines the influence of a very small Weissenberg number (We) on shear thickening or thinning, which is reliant upon its numerical values. when We>0, the fluid substance shows shear thickening behavior. Shear-thinning behavior is displayed by the fluid material when We<0. Examining the rheological behaviour at low We is crucial in coating applications, as the flow properties of shear-thickening or shear-thinning fluids significantly affect film uniformity and the efficiency of the coating process. Thus, We serves as the small dimensionless parameter controlling the perturbation expansion in this study. Eqs. (11) and (14) are challenging to solve in closed form as it is a non-linear differential equations. Consequently, the PM can be utilized to determine its analytical solution. Treating We<<1 as the perturbation parameter, we apply the conventional perturbation procedure and express the variables through the following series expansions:

{u(y)=We0u0(y)+We1u1(y)+…,λ=We0λ0+We1λ1+…,T(y)=We0T0(y)+We1T1(y)+…,x=We0x0+We11x1+…,P(x)=We0p0(x)+We1p1(x)+…(17)

Eqs. (11) and (14) are modified by inserting the relationships in Eq. (17) as given below

d2(u0We0+u1We1+…)dy2+3(We)(n−12)(d(u0We0+u1We1+…)dy)(d2(u0We0+u1We1+…)dy2)=d(P0We0+P1We1+…)dx+(F0We0+F1We1+…),(18)

d2(T0We0+T1We1+…)dy2+[1+We.(d(u0We0+u1We1+…)dy)2]n−12+Br.(d(u0We0+u1We1+…)dy)2=0.(19)

By collecting terms of identical orders in We, the following boundary value problems (BVPs) are obtained.

A. Zero-order velocity problem

dp0dx=d2u0dy2−Fo.(20)

Eq. (20) is subject to the following BCs:

u0(0)=0,u0(h)=1.(21)

In accordance with the BCs specified in Eq. (21), the solution to Eq. (20) is as follows:

u0(y)=(dp0dx+F0)y22−(h2F0+h2dp0dx−2)y2h.(22)

In the above formulation, u0 is treated as an implicit function of x while remaining an explicitly defined function of y.

Q0=∫xbh(x)u0(x,y)dy.(23)

The symbol Q0 is used to represent the baseline volumetric flow rate as established previously. By applying the principle of mass conservation, the pressure gradient in Eq. (22) can be determined, leading to the following expression.

dp0dx=−(1+x022)3F0−6−3x02+12λ0(1+x022)3.(24)

Assuming Eq. (23) xb→−∞, the derivation of the zeroth-order pressure distribution proceeds as follows:

p0(x)=14(x02+2)2(−18(x02+2)2(λ0−23)2arctan⁡(x022)−9(x02+2)2(λ0−23)π2−4(F0x04+(4F0+9λ0−6)x02+4F0+30λ0−12)x0).(25)

By symmetry considerations, the velocity must reduce to zero at the point of separation (xs), (xs,1/2h(xs)). Thus, it is simple to find that

λ0=−196(Fx4+4Fx2+4F−24)(x2+2),(26)

by substituting F=0, in Eq. (26), one can figure out that,

λ0=14x2+12.(27)

Substituting xs from the preceding relation into Eq. (25), and imposing the condition p0=0 at xs=x, yields the following result. The thickness of coating λ0=0.6056, is determined by using the Newton-Raphson technique at the point of separation xs=0.6544, and F=0.01. The solutions presented earlier correspond closely to those established by Middleman and Greener [8] in the absence of the gravity force that is F=0.

B. First-order velocity problem

d2u1dy2+3(n−1)[du0dy]2(d2u0dy2)2=dp1dx−F1.(28)

For Eq. (28) is subject to the following BCs:

u1(0)=0,u1(h)=0.(29)

The solution to Eq. (28) under the boundary conditions (29) is as follows:

u1(y)=116h3|h6np03y−3h5ny2p03+4h4ny3p03−2h3ny4p03−h6p03y+3h5y2p03−4h4y3p03+2h3y4p03−6h4np02y+12h3ny2p02−8h2ny3p02+6h4p02y−12h3y2p02+8h2y3p02+8p1y2h3+12h2np0y−12hny2p0−12h2p0y+12p0y2h−8ny+8y+(h4np02−h4p02−4h4p1+4n−4)y8h3|,(30)

dp1dx=340h7|F3h9n−F3h9−18F2h7n+36F2h6nλ0+18F2h7−36F2h6λ0+128Fh5n−432Fh4nλ0+432Fh3nλ02−128Fh5+432Fh4λ0−432Fh3λ02+160h4λ1−336h3n+1536h2nλ0−2592hnλ02+1728nλ03+336h3−1536h2λ0+2592hλ02−1728λ03|,(31)

in Eq. (30), P0=dp0dxandP1=dp1dx are the zeroth-order and first-order pressure gradient, respectively. The u1 is an implicit function of x and an explicit function of y. From the application of mass conservation, the first-order pressure profile in Eq. (30) is obtained. Integrating Eq. (31) with p(x)=0 as a boundary condition as xb→−∞ provides the pressure distribution p1(x) of the first order. The Swift-Stieber BCs for pressure at the transition point xs=x are utilized for calculating the value of λ1.

3.1.1 Temperature Distribution

By simplifying and substituting Eq. (14) into Eq. (17), we get a zeroth- and first-order temperature distribution by equating the corresponding powers of We:

d2T0dy2+Br.(du0dy)2=0,(32)

d2T1dy2=−Br[2(du0dy)(du1dy)+n−1(du0dy)42].(33)

A. Zero-order temperature problem

The zeroth-order temperature Eq. (32), has the dimensionless BCs:

T0(0)=0,T0(h)=1.(34)

Taking into consideration the BCs given in Eq. (34), the solution to Eq. (32) is given below:

T0(y)=|−Br.(dp0dxh2−2dp0dxyh−2)4192h4p02+(Brh4dp02dx2+12dp1dxh2+4Br)y12h3dp0dx+Br.(h8dp04dx4−8h6dp03dx3+24h4dp02dx2−32dp0dxh2+16)192h4dp02dx2|(35)

B. First-order temperature problem

The Eq. (33) is the first-order temperature equation followed with the subsequent dimensionless BCs.

T1(0)=0,T1(h)=0.(36)

The first-order solution employs a similar methodology to derive the equations for pressure profile, pressure gradient, and detachment locations, quite comparable to the zero-order case. Through the superposition of zero- and first-order solutions, the complete analytical expressions for velocity and temperature are established. The availability of pressure distributions, velocity profiles, and pressure gradients permits expedited determination of critical operating parameters, namely Nusselt number, power input, and the force of roll splitting.

4  Operating Factors

After establishing the pressure gradient, velocity profile, and pressure distributions, engineering parameters such as power input and separation force can be efficiently computed.

4.1 The Force of Roll Separation (Sf)

Sf=ρgH0μU=∫−∞xsp(x)dx(37)

4.2 The Power Input (Pw)

Pw=∫−∞xsτxy(x,1)dx(38)

4.3 The Nusselt Number (Nu)

Nu=|∂T∂y|y=h(39)

5  Optimization Section

The present investigation requires thirteen runs and twelve degrees of freedom, together with three separate parameter levels, for the use of RSM. The independent variables employed in this study include λ, Pw, and Sf. The variables and their corresponding levels are presented in Table 1. The table presents the low (−1), center (0), and high (1) values of several factors.

Response Surface Method (RSM)

RSM is an optimization method utilized to identify the combination of input parameters that minimizes or maximizes the target function (see Flowchart 1) [29]. The RSM findings provide a practical framework for enhancing roll-coating processes. RSM formulates regression-based models for coating thickness, power input, and separation force to determine the ideal combination of process parameters, including the Weissenberg number and body-force term, that ensures uniform coating while reducing energy consumption. The response surfaces and contour plots illustrate the trade-offs between coating quality and power requirements, allowing process engineers to choose operating conditions that preserve the desired film thickness and stability while minimizing mechanical wear and energy use.

Flowchart 1: Flow chart describing RSM

A central composite design (CCD) that is face-centered is the subject of this study. The CCD approach illustrates the correlation between many predictor factors and one or more responder variables. The inquiry includes two input factors, We, F, and three response variables, Pw, λ and sf. For determining the model’s suitable response, the use of a central composite face (CCF) incorporates quadratic, interaction, and linear components. Table 2 presents the experimental design and its corresponding responses, whereas the comprehensive quadratic model has linear, squared, intercept, and interaction.

ℜ=γ0+γ1A+γ2B+γ3A2+γ4B2+γ5AB+E,(40)

where the variables used for regression analysis are denoted as γ0…γ5, while physical factors are indicated by A and B, E signifies the residual error, and the ℜ-response variables are referred to as objective functions.

A. Analysis of variance (ANOVA)

To find out a significant difference in the means of three or more groups, statisticians utilize a process called analysis of variance (ANOVA). We computed and compared F-ratios and p-values. The following part provides a concise overview of the calculating methods mentioned by Montgomery [33]. By controlling for variability, this strategy assists researchers in understanding the effects of one or more factors. Total sum of squares (SST), treatment sum of squares (SSTreatment), and error sum of squares (SSE) can be calculated respectively as follows:

{SST=∑i=1k∑j=1ni(Yij−Y¯)2,SSTreatment=∑i=1Kni(Y¯i−Y¯)2,SSE=∑i=1K∑j=1ni(Yij−Y¯i)2,(41)

where K is the number of groups, ni is the number of observations in group i, Y¯ is the overall mean, and Yij is the jth observation obtained under the i in this context. The following equations can be used to compute the mean square (MS), F0, and p-value:

MS=SSDf,F0=MSTreatmentsMSE,p−value=P(F(K−1,n−1)>calculatedF),(42)

where, Df is the degree of freedom, if Fo>Fα,K−1,k(n−1) at α=0.05 is the level of significance then the corresponding effect is considered statistically significant. An alternate method is p-value. We retained the terms with p-value less than 0.05 because we believed they were significant. To ensure the model was adequate, residual analysis was performed. The experimental design consists of thirteen runs and incorporates twelve degrees of freedom (Df). The outcomes for the variables are documented according to the coded values (refer to Table 2). The model’s accuracy is evaluated via an analysis of variance.

The results presented employ adjusted MS, Df, SST, p-value, and F-value as statistical estimators for the analysis. The F-value serves as a quantifiable metric for assessing variation in data concerning the mean. The ANOVA analysis indicates that the input data exhibit higher precision when F-values exceed one.

ANOVA evaluated the model’s adequacy and variable significance via the F-statistic and p-values. The exceptional model F-value (5487.26) and p-values below 0.1000 (Table 3) confirm the model’s high significance and the statistical relevance of its terms. The ANOVA results, confirming the statistical significance of the selected independent variables (Sf, Pw and λ), are summarized in Tables 4 and 5. A higher F-value corresponds to a more influential effect. A p-value less than 0.05 indicates statistical significance, whereas a p-value greater than this value suggests that the model terms are not significant.

A larger coefficient of determination (R2-value) signifies a better fit of the regression model to the observed data. A control variable in this model was the R2. Referring to Table 6 the model reports higher R2 values of 99.97% for λ, 99.62% for Pw and 99.99% for Sf which suggests that the regression model provide an excellent relationship between the independent variables and the dependent parameters.

In this study, the value of R2 adjusted (adj) for λ is 99.96% Pw is 99.35% and for Sf was 99.98% found. There is a strong connection between the anticipated and observed values, as denoted by the values of R2 (adj.), R2 being near to 100%. This means that the regression model works well.

B. Regression analysis and model accuracy

Multiple linear regression models that are applicable to the selected range of values are used to estimate λ, Sf and Pw in this subsection. A total of nine sets of values were selected from the available list of [0.01,0.09], [0.1,0.9] for F and We, respectively. The ultimate accuracy of the regression model is ascertained by the ANOVA tables, with the formulae for parameters λ,Pw and Sf specified as follows:

λ=0.613759+0.03592A−0.7764B−0.00008A*A+0.836B*B−0.2828A*B,(43)

Pw=0.99464+0.1372A+3.654B+0.0199A*A−15.57B*B−1.289A*B,(44)

Sf=1.05054−0.3201A+10.758B−0.0014A*A−27.86B*B+6.030A*B.(45)

Residual analysis is essential for validating regression model accuracy by diagnosing violations of key statistical assumptions. The Goodness of Fit plot evaluates predictive alignment, while the Normal Probability Plot tests residual normality. Collectively, these graphs identify model inadequacies, refine structure, and ensure robust statistical inferences. Their interpretation is critical for optimizing performance and justifying the real-world applicability of the regression framework.

The residual plots may provide additional evidence of the model’s accuracy. Fig. 3a–c illustrates a major connection between the observed and fitted values as shown in the residuals vs. fitted values graphs. The highest variances are seen to be as little as 0.0005, 0.2, and 0.004, respectively, for λ, Pw and, Sf. Therefore, the RSM model is both accurate and significant. A two-dimensional visual aid in Fig. 3a–c the normal probability plots indicate that all points align closely with the straight line, exhibiting little variation, the results demonstrate that the data is normally distributed, which illustrates the link between two independent variables and two outcome variables is called a model contour plot. Fig. 4a–f depicts the response surface and model contour plot, which show how several independent factors, like We and F in our example, affect the λ,Pw and Sf Contour plots, together with their corresponding three-dimensional surface representations, illustrate the variation of a function or dataset across the parameter space.

Figure 3: (a). Residual plots for λ. (b). Residual plots for λ. (c). Residual plots for Sf

Figure 4: (a–f). Plots in two and three-dimensions showing the impact of surface response on Pw (a,b), λ (c,d) and Sf (e,f)

Contour plots and their corresponding 3-D surface representations illustrate how a function or dataset varies across the domain. In a response surface plot, the x-axis and the y-axis, represent the input variables, while the z-axis depicts the resulting response. Each graph clarifies the evolving relationships between inputs and outputs, demonstrating tendencies such as saturation, diminishing returns, or direct growth. In Fig. 4b, the λ increases with increasing values of We and F, while the rate of growth diminishes at higher λ values. In Fig. 4d, Pw increases consistently when We and F, grows, demonstrating a clear positive relationship. Similarly, Fig. 4f expresses a positive relationship for increasing the parameters like We and F.

C. Sensitivity analysis

Sensitivity analysis is the initial and important stage in solving optimization problems, as it provides insights into the increasing or decreasing tendencies of the objective function relative to the design parameters. This approach determines essential criteria and prioritizes them according to their significance [29]. The sensitivity of the design objective function to a design variable is defined as the partial derivative of that function with respect to its variables. The estimation of sensitivity is conducted by evaluating the partial derivatives of the response functions (Sf and λ,Pw) in relation to the independent variables (B and A). The pertinent effective parameters are employed to assess these derivatives, thereby establishing the sensitivity functions.

We compute the partial derivative, or sensitivity function, utilizing the pertinent effective parameters.

{∂(λ)∂A=0.03592−0.2828B+0.00016A,∂(λ)∂B=−0.7764+1.672B−0.2828A,∂(Pw)∂A=0.1372−1.289B+0.0398A,∂(Pw)∂B=3.654−1.289A−31.14B,∂(Sf)∂A=−0.3201+0.0028B+6.030A,∂(Sf)∂B=10.758+6.030A−55.72B,(46)

where ∂(λ)∂A,∂(λ)∂B,∂(Pw)∂A,∂(Pw)∂B,∂(Sf)∂A,∂(Sf)∂B are the functions of sensitivity. Table 7 illustrates the sensitivity of the output to variations in input parameters (F, We) across several levels. Through the analysis of derivatives, engineers can acquire essential insights into the impact of parameter variations influence key output characteristics. Sensitivity is evaluated using the regression Eq. (46).

The sensitivity of the key parameters (λ, Pw, Sf) was evaluated at distinct levels of A (−1, 0, 1), as shown in Figs. 5–7. This analysis reveals that a positive sensitivity corresponds to an enhancement of the output function.

Figure 5: λ sensitivity analysis at (a) A = −1, (b) A = 0, and (c) A = 1, respectively

Figure 6: Pw sensitivity analysis at (a) A = −1, (b) A = 0, and (c) A = 1, respectively

Figure 7: Sf sensitivity analysis at (a) A = −1, (b) A = 0, and (c) A = 1, respectively

6  Validation through Numerical Methods

In this section, the results obtained from numerical and analytical approaches are systematically compared. The numerical solutions are computed using the Finite Difference Method (FDM) alongside the BVP-Dsolve solver (Midrich method) implemented in Maple software, enabling an accurate and efficient resolution of the governing equations. The advantages of the BVP Midrich solver are summarized in Flowchart 2. The approximate solution for the temperature and velocity profile is achieved through the PM.

Flowchart 2: Midrich BVP advantages

For a variety of non-linear boundary value problems, there is a numerical approach called the BVP Midrich. A variant of Euler’s midpoint technique is employed in this strategy. Up to a value of 1 × 10−6, this method demonstrates absolute error convergence. FDM is one of the most effective numerical methods for solving nonlinear coupled ODEs (12) and (15) with suitable boundary conditions. If relevant, it adjusts the differential equations by replacing derivatives with physical domain or time-based interval finite difference approximations. There is a relationship between the number of nodes n and the value of the step size h. The relationship between them is Δy=b−an, A greater number of nodes n leads to a shorter step interval h, resulting in higher accuracy outcomes. The flow chart for the FDM is expressed below in Flowchart 3. The next finite difference mathematical equations are employed to replace Eqs. (12) and (15), respectively:

dudy=ui+1−ui−12.Δy,d2udy2=ui+1−2.ui+ui−1Δy2,(47)

dTdy=Ti+1−Ti−12.Δy,d2Tdy2=Ti+1−2.Ti+Ti−1Δy2.(48)

Flowchart 3: Flow chart explaining FDM

The resulting non-linear system was solved iteratively using the Newton iteration approach. The specifics of the linearized system’s construction are left out to keep things brief. We investigate to ensure the analytical and approximate solutions are accurate before giving the results and discussing them. For this purpose, Figs. 8–16 illustrate the velocity and temperature results obtained from both numerical and analytical solutions, considering parameters such as F,and We (non-Newtonian parameter), and Br (Brickman number). The numeric calculations were conducted at fixed values of dpdx=−1.6,We=0.1. The discrepancy observed when contrasting the analytical solution with the approximate solution is illustrated in Figs. 17 and 18.

Figure 8: Comparison of u(y) at n = 0.5, We=0.1,0.5,0.9

Figure 9: Comparison of u(y) at n = 2, We=0.1,0.5,0.9

Figure 10: Comparison of u(y) at n = 0.9 F=0.01,0.05,0.09

Figure 11: Comparison of u(y) at n = 2 F=0.01,0.05,0.09

Figure 12: Comparison of T(y) at n = 0.5, We=0.1,0.5,0.9

Figure 13: Comparison of T(y) at n = 2, We=0.1,0.5,0.9

Figure 14: Comparison of T(y) at n = 0.9 F=0.01

Figure 15: Comparison of T(y) at n = 2 F=0.01

Figure 16: Comparison of T(y) at Br=0.1,0.5,0.9

Figure 17: Comparison of errors between analytic and numerical solutions

Figure 18: Analysis of the error behavior in the analytical vs. numerical solutions

The suggested numerical method yields both the exact and approximate solutions, and the error is visible in both cases. Solid lines indicate FDM findings, solid boxes show results from the integrated BVP approach, and dashed lines show analytical results in the numerical results. The use of numerical and analytical solutions is justified because of their clearly superior performance in terms of accuracy. Three curves, each in black, red, and blue, make up each design. The solutions for one parameter value are shown by the black curve, and the solutions for two more parameter values are shown by the red and blue curves. The visuals demonstrate a robust interrelationship within all proposed solutions. This reinforces confidence in the reliability of all results. Tables 8 and 9 present an analytical and numerical comparison of velocity and temperature measurements at We = 0.1, F = 0.01, and the corresponding absolute errors are also reported. The tables demonstrate a decrease in the absolute error when comparing the analytical and numerical results. It proves the validity of our technique and the related computations. This proves that the methods we used to achieve the results are reliable. The influence of physical characteristics on several important variables will be discussed in the next section using tables and graphs. The analytic and numerical solutions exhibit convergence.

Validation of the Present Model

The accuracy of the present formulation was verified by comparing results with the classical Newtonian model of Middleman [8]. As shown in Figs. 19–22, the velocity, temperature, and pressure profiles show excellent agreement with the established Newtonian data at We=0. When We,F→0, the present non-Newtonian model exactly reduces to Middleman’s [8] Newtonian solution, confirming its mathematical consistency. This close agreement validates both the analytical and numerical approaches employed in the current study.

Figure 19: The We effect on u(y) at n = 0.5, x = 0

Figure 20: The We impact on dpdx

Figure 21: The We impact on P(x)

Figure 22: The We effect on T(y) at n = 0.5, x = 0

7  Results and Discussion

The mathematical modeling of CRF flow is analyzed under fixed boundary constraints to establish a rigorous framework for a non-Newtonian, incompressible fluid traversing the narrow gap between the roll and substrate. This study systematically examines the influence of essential governing parameters on the spatial distributions of velocity, pressure, temperature, and other vital engineering variables. The research clarifies the complex relationship between flow dynamics and thermal transport, yielding essential insights into the fundamental physical mechanisms, offering a predictive basis for enhancing CRF processes. The governing equations are simplified with the application of LAT. Closed-form solutions are derived analytically utilizing PM for the velocity profile, the temperature distribution for the pressure distribution, and λ. The values of the parameters utilized in the simulations were selected based on their relevance to typical industrial roll-coating processes and to maintain numerical stability within physically realistic limits. The 0.1≤We≤0.9 was chosen to represent weak to moderate elastic effects commonly observed in bio-based non-Newtonian fluids such as sodium alginate. The body-force parameter 0.01≤F≤0.09 corresponds to practical ranges of gravitational or pressure-driven effects in coating gaps. The pressure gradient was fixed at dpdx=−1.6, while. Br=1. These parameter choices ensure that the predicted flow conditions align with experimentally realizable roll-coating systems involving viscous, shear-thinning fluids. The We characterizes the ratio of elastic to viscous forces in a non-Newtonian fluid. As We increases, elastic effects become more dominant, and the polymer chains in the fluid tend to stretch and align along the flow direction. This molecular alignment reduces the apparent viscosity, demonstrating the shear-thinning nature of the Carreau fluid. Due to the reduced viscosity, the fluid near the moving roll experiences less resistance, which allows more fluid to be entrained into the coating gap. Consequently, the coating or film thickness increases with higher We, while the velocity and pressure gradients decrease because of the enhanced elastic resistance within the flow. Thus, a rise in We amplifies shear-thinning behavior and results in a thicker coating layer on the roll surface.

Numerous flow parameters are presented through tables and graphical representations. Figs. 19 and 22–38 show how different parameters affect the velocity and temperature of the RC flow. Likewise, Figs. 20, 21, 39 and 40 show the pressure distribution and gradient. Moreover, the results for the xs, λ, Sf, Nu, and Pw are also carefully tabulated and examined in depth. The ranges of all physical variables employed in the RC model are summarized in Table 10. Lastly, when We,Br,F→0, the present approach is analogous to the Middleman [8] solution (Newtonian model) under non-isothermal conditions. The comparison between the numeric and analytic solutions is illustrated in the Figs. 8–16.

Figure 23: The We effect on u(y) at n = 0.5, x = 0.5

Figure 24: The We effect on u(y) at n = 2, x=0

Figure 25: The We effect on u(y) at n = 2, x=0.5

Figure 26: The F effect on u(y) at n = 0.5, x=0

Figure 27: The F effect on u(y) at n = 0.5, x=0.5

Figure 28: The F effect on u(y) at n = 2, x=0

Figure 29: The F effect on u(y) at n = 2, x=0.5

Figure 30: The We effect on T(y) at n = 0.5, x = 0.5

Figure 31: The We effect on T(y) at n = 2, x = 0

Figure 32: The We effect on T(y) at n = 2, x = 0.5

Figure 33: The F effect on T(y) at n = 0.5, x = 0

Figure 34: The F effect on T(y) at n = 0.5, x = 0.5

Figure 35: The F effect on T(y) at n = 2, x = 0

Figure 36: The F effect on T(y) at n = 2, x = 0.5

Figure 37: The Br effect on T(y) at x = 0

Figure 38: The Br effect on T(y) at x = 0.5

Figure 39: The F impact on dpdx

Figure 40: The F impact on P(x)

A. Velocity profile

The velocity profile shows a significant enhancement with increasing We and F at n=0.5, attributable to shear-thinning behavior and reduced viscous resistance. This effect is particularly valuable in high-speed coating applications such as paper manufacturing and flexible electronics, where uniform velocity profiles are critical for achieving consistent coating quality.

The results for the dimensionless velocity profiles as a function of y for the relevant parameters during the RCP are presented in Figs. 19 and 23–29. An improvement in the peak velocity inside the fluid flow is noted when the value of We is increased in the interval [0.1,0.9]. Higher We values amplify the fluid’s resistance to deformation, leading to a thicker coating layer and increased flow near the roll surface. An opposite trend can be seen in the graphs for We=[0.1,0.9] at n>1. Figs. 19 and 23 demonstrate how the non-Newtonian parameter effects on velocity at n=0.5, x=0,0.5. These illustrations show that for n<1, the speed of increasing argument values is negligible, but for Figs. 24 and 25 n>1, the opposite is true. For the CRF model, n > 1 corresponds to shear-thickening behavior, whereas n < 1 corresponds to shear-thinning behavior.

For a shear-thinning fluid (n<1), an increasing We enhances elastic effects, causing a significant reduction in apparent viscosity near the walls under high shear. This enables the fluid to slip more easily along the boundaries, resulting in a blunted, plug-like velocity profile where the core flows with greater uniformity and higher speed. Consequently, the overall fluid through ought increases while requiring less power input for the coating process. While a shear-thickening fluid (n>1), an increasing We dramatically increases the apparent viscosity near the walls under high shear. This viscous resistance effectively constricts the flow, forcing more fluid to move through the center of the gap. This behavior leads to higher power requirements and potential flow instabilities in the coating process. Physically, a higher We value indicates a thicker liquid since it increases the ratio of viscous forces to elastic forces, which in turn increases the concentration of the liquid. The velocity profile increases steeper for n<1 and generally flattens out for n>1 because the flow becomes more resistive. Figs. 26 and 27 display the fluid velocity by changing F from 0.01 to 0.09 at n=0.5 and x=0,0.5. The fluid velocity increases at large values of the corresponding parameter F. Moreover, velocity is zero near the roll tip, while by increasing, F maximum change occurs in velocity at its center both in the plane and the coater. Contrary behaviour can be witnessed in Figs. 28 and 29. The following is an analysis of how the given parameters affect the pressure gradient and pressure profile.

B. Pressure profile and gradient

The pressure distribution exhibits a strong dependence on rheological parameters, diminishing with increasing We but augmenting with rises F. These findings have direct implications in automotive paint, shops, and industrial coating lines, where pressure management is essential for minimizing separation forces and reducing wear on mechanical components. Newtonian [8] outcomes are attained as We→0,F→0,Br→0. The upward rotation of the roll creates a pressure gradient dpdx inside the fluid confined between the roll and the wall, which fundamentally governs the flow dynamics. The distribution of pressure and its gradient are essential in actual manufacturing, as they directly influence the necessary rolling force and power consumption of the roll.

Understanding how parameters like We and F influence pressure enables operators to optimize processing conditions, prevent defects like blistering or uneven thickness, and significantly reduce energy costs. Consequently, these figures provide essential practical insight for effective and high-quality manufacturing. point. For specific values of F in the interval [0.01, 0.09], one can observe that dpdx is declined by increasing F and We in the interval [0.1, 0.9]. Figs. 20 and 39 depict the pressure gradient (dpdx) as a function of the axial coordinate x for the relevant parameters. The symmetric configuration can be discovered at the nip zone x=0. A region where the pressure gradient is negative (dpdx<0), called the nip region, dpdx rises uniformly to a higher value, and declines gradually until it reaches zero at the separating.

Figs. 21 and 40 depict the pressure distribution for various values of We and F, respectively. It has been seen that the magnitude of the pressure distribution declines when We increases, while for enhancing the values of F the pressure distribution increased. The pressure escalates along the negative x-axis, attaining its peak at the critical point x=−1.5, after which it begins to decrease. The F directly influences the pressure distribution and load-carrying capacity within the coating film. As F increases, the body force acting along the flow direction intensifies, resulting in a higher-pressure accumulation in the nip region between the roll and the wall. This augmented pressure gradient accommodates a larger portion of the applied load, thereby enhancing the load-carrying capability of the fluid film. Physically, a more gravitational or body-force effect pushes more fluid toward the converging zone, thickening the film and elevating local pressure. Consequently, higher Fvalues yield enhanced hydrodynamic lift and improved bearing action, while lower F values result in decreased pressure and thinner coating layers.

C. Temperature profile

The theoretical investigation of the RC procedures for fluids that are non-Newtonian is greatly affected by temperature. Temperature significantly influences the rheological properties of fluids and the coating process as a whole. The temperature distribution is significantly influenced by both rheological and operational parameters. The roll has zero temperature at the edge, but reaches its maximum temperature in the middle. The rise in temperature with increasing We, F, and Br, highlights the importance of thermal management in processes such as pharmaceutical tablet coating and food packaging, where precise temperature control is crucial for product quality and safety. The influence of relevant parameters such as, We and F on a dimensionless temperature distribution is illustrated in Figs. 22 and 30–38 at different positions in the nip area (x=0,0.5) in the RC process at Br=1.

The temperature distribution is observed in Figs. 22 and 30 to decline at various coating process locations with increasing the values of We by using value of n=0.5 to conclude the results for n<1.

While an opposite trend is witnessed in Figs. 31 and 32 for n>1, a specific value for n=2 is used for these results. In Figs. 33–38, it is clearly observed that the temperature distribution exhibits a strictly increasing trend with respect to the specified parameter F, within the range [0.1, 0.9]. This indicates that as the parameter increases, the temperature at all points within the domain consistently rises, reflecting a direct and monotonic relationship.

Figs. 37 and 38 illustrate how Br affects the temperature distribution. It is evident that the temperature profile increases monotonically with Br. Figs. 12–15 show that the numerical temperature solution confirms the same trend.

Additional quantities of engineering significance, including xs, λ, Sf, Nu, and Pw, are presented in Tables 11 and 12 of this analysis. In Table 11, it is observed that as the value of F increases, the values of xs and λ decrease, while Sf, Nu, and Pw, increase. The effects of this We have been examined in Table 12. It has been noted that all quantities of interest, with the exception of Nu, have shown an increase as the value of it is raised. The results of this study align with the conclusions drawn by Sullivan and Middleman [8] as all additional parameters tend toward zero.

D. Streamlines

Streamlines are a fundamental concept in fluid mechanics, describing the path that fluid particles adhere to in a steady flow. They are defined as curves that are consistently tangent to the velocity vector of the flow, indicating that a fluid element traverses along a streamline without intersecting it.

A fundamental attribute of streamlines is their non-intersection, as the velocity at any location in the flow field has a unique direction. The flow behavior of the fluid can be illustrated by creating streamlined patterns. This mathematical formulation ensures the automatic satisfaction of the continuity equation. They illuminate the overall flow structure, encompassing the presence of areas with elevated or diminished velocity. The streamline plots in Figs. 41–44 illustrate that the flow displays symmetry about the vertical y-axis, so confirming the chosen geometric domain and boundary conditions. This symmetry signifies that the physical forces and boundary effects are equilibrated on both sides of the axis, a typical feature in well-posed fluid dynamics problems.

Figure 41: Streamline for We=0.1 at F = 0.01

Figure 42: Streamline for We=0.1 at F = 0.05

Figure 43: Streamline for We=0.5 at F = 0.01

Figure 44: Streamline for We=0.9 at F = 0.01

8  Conclusions

This study effectively addresses a significant gap in the literature by formulating a comprehensive model for non-isothermal Carreau fluid flow in a roll coating process with a fixed constrained wall, a configuration pertinent to industrial thickness control. The primary contributions and engineering results are summarized as follows:

(a)   Innovative Modeling and Validation: We developed a novel hybrid analytical-numerical framework for this previously unexplored configuration. The excellent agreement between our perturbation solutions and numerical results (FDM, BVP), as vividly demonstrated in Figs. 8–18, confirms the robustness and accuracy of our model. Furthermore, Figs. 19–22 validate our approach by showing perfect agreement with the classical Newtonian model for the limiting case (We→0).

(b)   Key Engineering Outcomes and Design Guidance: The parametric analysis, illustrated through extensive figures, provides direct design insights:

Coating Thickness Control: Figs. 19 and 23 show that for shear-thinning fluids (n<1), increasing the We promotes a uniform, plug-like velocity profile and significantly enhances coating thickness. Recommendation: Use shear-thinning fluids and operate at higher We to achieve thicker coatings with less power.

Energy and Load Management: Figs. 20, 21, 39 and 40 reveal that pressure and separation force can be managed by tuning We and F. Higher F increases load capacity but also the separation force. Recommendation: To minimize wear and energy consumption, optimize We and F parameters using the provided RSM models.

Thermal Management: Figs. 22 and 30–38 indicate that temperature rises with the Br and F. Recommendation: In temperature-sensitive processes, carefully control Br and F to prevent product degradation.

(c)   Practical Optimization Framework: We employed Response Surface Methodology (RSM) to transform our findings into a practical optimization tool. The high coefficients of determination Table 6 and the clear trends in the 3-D response surfaces (Fig. 4) allow engineers to identify the optimal combination of We and F to maximize coating thickness while simultaneously minimizing power input and separation force.

(d)   Future work and its applicability: In the future, we will be able to validate this study for a variety of engineering studies, both experimental and theoretical, that include complex structures. (Newtonian and non-Newtonian fluids) by applying slip/no-slip, MHD, and many other fluid properties. We will extend the present model to multi-dimensional and transient analyses using advanced numerical methods and integrate machine learning for faster and more accurate predictions. The extended model will help optimize coating thickness, flow uniformity, and heat transfer in practical non-Newtonian roll-coating systems.

(e)   Limitations of the Present Study: Although the present study establishes a comprehensive analytical-numerical framework for non-isothermal Carreau fluid flow in roll coating with a fixed constrained wall, several inherent limitations remain. The formulation relies on the lubrication approximation and assumes a thin, steady, laminar film with predominantly unidirectional flow, which may not fully capture high-speed or inertia-dominated coating conditions. The perturbation methodology is applicable only for small We, restricting its accuracy for strongly elastic or highly non-linear Carreau fluids. Thermal effects are modeled with constant properties and moderate temperature gradients, whereas real coating materials often exhibit temperature-dependent viscosity, density, and conductivity. In addition, the model excludes free-surface deformation, surface tension, slip effects, roll elasticity, and fluid-substrate interactions, all of which can influence industrial roll-coating performance. Moreover, the RSM-based optimization is performed inside a constrained parametric domain, and the resultant predictions should not be generalized without additional validation. These considerations outline the applicability of the existing framework and propose natural directions for further work, including transient simulations, temperature-dependent material characteristics, free-surface dynamics, roll compliance, and enhanced multi-physics integration.

In summary, this paper provides a validated theoretical mode, along with explicit directions and a practical framework for optimizing industrial roll coating processes. The results indicate that by strategically choosing the fluid’s rheological characteristics (shear-thinning behavior) and operating conditions (moderate to high We), significant enhancements in coating uniformity, process efficiency, and equipment durability can be achieved.

Acknowledgement: Not applicable.

Funding Statement: This study is supported by the Talent Project of Tianchi Young-Doctoral Program in Xinjiang Uygur Autonomous Region of China (No. 51052401510), and Natural Science Foundation General Project (Grant Number 2025D01C36) of the Xinjiang Uyghur Autonomous Region of China. Also, This study received financial support from the National Natural Science Foundation of Xinjiang Province (Grant Nos. 2022TSYCTD0019 and 2022D01D32), the China Scholarship Council (CSC) (Grant No. 2021SLJ009915).

Author Contributions: The concept was developed by Fateh Ali, M. Zahid and Mujahid Islam, and subsequently incorporated into a mathematical model by Mujahid Islam and Fateh Ali. The problem was solved by Mujahid Islam. The manuscript was drafted by Fateh Ali, Xinlong Feng and Sana Naz Maqbool, with supervision provided by Xinlong Feng. Mujahid Islam and Fateh Ali analyzed the underlying physics and contributed to the discussion. Finally, Xinlong Feng and M. Zahid reviewed, edited, and assisted with English language corrections. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Relevant data for this study are accessible from the corresponding authors upon reasonable inquiry.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

Nomenclature

References

1. Hermia J. Blocking filtration. application to non-Newtonian fluids. In: Mathematical models and design methods in solid-liquid separation. Berlin/Heidelberg, Germany: Springer; 1985. p. 83–9. doi:10.1007/978-94-009-5091-7_5. [Google Scholar] [CrossRef]

2. Awais M, Salahuddin T, Muhammad S. Evaluating the thermo-physical characteristics of non-Newtonian Casson fluid with enthalpy change. Therm Sci Eng Prog. 2023;42:101948. doi:10.1016/j.tsep.2023.101948. [Google Scholar] [CrossRef]

3. Toghraie D, Esfahani NN, Zarringhalam M, Shirani N, Rostami S. Blood flow analysis inside different arteries using non-Newtonian Sisko model for application in biomedical engineering. Comput Meth Programs Biomed. 2020;190:105338. doi:10.1016/j.cmpb.2020.105338. [Google Scholar] [PubMed] [CrossRef]

4. Liu TJ, Lin HM, Hong CN. Comparison of two numerical methods for the solution of non-Newtonian flow in ducts. Num Meth Fluids. 1988;8(7):845–61. doi:10.1002/fld.1650080707. [Google Scholar] [CrossRef]

5. Yasir M, Ahmed A, Khan M, Ch MMI, Ayub M. Study on time-dependent Oldroyd-B fluid flow over a convectively heated surface with Cattaneo-Christov theory. Waves Random Complex Medium. 2025;35(1):144–61. doi:10.1080/17455030.2021.2021316. [Google Scholar] [CrossRef]

6. Nayfeh AH. Perturbation methods. Hoboken, NJ, USA: John Wiley & Sons; 2024. [Google Scholar]

7. Coyle DJ. The fluid mechanics of roll coating: steady flows, stability, and rheology. Minneapolis, MN, USA: University of Minnesota; 1984. [Google Scholar]

8. Sullivan TM, Middleman S. Roll coating in the presence of a fixed constraining boundary. Chem Eng Commun. 1979;3(6):469–82. doi:10.1080/00986447908935879. [Google Scholar] [CrossRef]

9. Greener Y, Middleman S. Blade-coating of a viscoelastic fluid. Polym Eng Sci. 1974;14(11):791–6. doi:10.1002/pen.760141110. [Google Scholar] [CrossRef]

10. Hintermaier J, White R. The splitting of a water film between rotating rolls. Tappi J. 1965;48(11):617–25. [Google Scholar]

11. Greener J, Middleman S. Reverse roll coating of viscous and viscoelastic liquids. Ind Eng Chem Fund. 1981;20(1):63–6. doi:10.1021/i100001a012. [Google Scholar] [CrossRef]

12. Owens MS, Vinjamur M, Scriven LE, Macosko CW. Misting of non-Newtonian liquids in forward roll coating. J Non Newton Fluid Mech. 2011;166(19–20):1123–8. doi:10.1016/j.jnnfm.2011.06.008. [Google Scholar] [CrossRef]

13. Lécuyer HA, Mmbaga JP, Hayes RE, Bertrand FH, Tanguy PA. Modelling of forward roll coating flows with a deformable roll: application to non-Newtonian industrial coating formulations. Comput Chem Eng. 2009;33(9):1427–37. doi:10.1016/j.compchemeng.2009.04.001. [Google Scholar] [CrossRef]

14. Ascanio G, Carreau PJ, Tanguy PA. High-speed roll coating with complex rheology fluids. Exp Fluids. 2006;40(1):1–14. doi:10.1007/s00348-005-0025-5. [Google Scholar] [CrossRef]

15. Echendu SOS. Computational and rheological studies for coating flows. Swansea, UK: Swansea University; 2013. [Google Scholar]

16. Raske N, Hewson RW, Kapur N, de Boer GN. A predictive model for discrete cell gravure roll coating. Phys Fluids. 2017;29(6):062101. doi:10.1063/1.4984127. [Google Scholar] [CrossRef]

17. Abdal S, Taha T, Ali L, Zulqarnain RM, Yook SJ. Neural networking-based approach for examining heat transfer and bioconvection in non-Newtonian fluid with chemical reaction over a stretching sheet. Case Stud Therm Eng. 2025;69:106047. doi:10.1016/j.csite.2025.106047. [Google Scholar] [CrossRef]

18. Rajagopal KR. Mechanics of non-Newtonian fluids. In: Recent developments in theoretical fluid mechanics. Boca Raton, FL, USA: Chapman and Hall/CRC; 2023. p. 129–62. doi:10.1201/9781003417026-5. [Google Scholar] [CrossRef]

19. Shahzad MH, Nadeem S, Awan AU, Allahyani SA, Ameer Ahammad N, Eldin SM. On the steady flow of non-Newtonian fluid through multi-stenosed elliptical artery: a theoretical model. Ain Shams Eng J. 2024;15(1):102262. doi:10.1016/j.asej.2023.102262. [Google Scholar] [CrossRef]

20. Khaliq S, Abbas Z. Analysis of calendering process of non-isothermal flow of non-Newtonian fluid: a perturbative and numerical study. J Plast Film Sheeting. 2021;37(3):338–66. doi:10.1177/8756087920979024. [Google Scholar] [CrossRef]

21. Jang J-Y, Chen P-Y. Reverse roll coating flow with non-Newtonian fluids. Commun Comput Phys. 2009;6(3):536. [Google Scholar]

22. Shahzad H, Wang X, Hafeez MB, Shah Z, Alshehri AM. Study of slip effects in reverse roll coating process using non-isothermal couple stress fluid. Coatings. 2021;11(10):1249. doi:10.3390/coatings11101249. [Google Scholar] [CrossRef]

23. Devisetti SK, Bousfield DW. Fluid absorption during forward roll coating of porous webs. Chem Eng Sci. 2010;65(11):3528–37. doi:10.1016/j.ces.2010.02.042. [Google Scholar] [CrossRef]

24. Bhatti MM, Ellahi R, Sait SM. Numerical investigation of magnetized nanofluid flow in a non-Darcian medium with convective heating conditions across a rotating disk. Int J Numer Meth Heat Fluid Flow. 2025;35(5):1736–63. doi:10.1108/hff-03-2025-0135. [Google Scholar] [CrossRef]

25. Farooq A, Rehman S, Ma WX. Comparative study of regression-based data-driven models for thermally stratified Carreau nanofluids with magnesium oxide nanoparticles. Phys Fluids. 2025;37(8):082032. doi:10.1063/5.0274063. [Google Scholar] [CrossRef]

26. Anderson MJ, Whitcomb PJ. RSM simplified: optimizing processes using response surface methods for design of experiments. Boca Raton, FL, USA: Productivity Press; 2016. [Google Scholar]

27. Milani Shirvan K, Mamourian M, Mirzakhanlari S, Ellahi R. Two phase simulation and sensitivity analysis of effective parameters on combined heat transfer and pressure drop in a solar heat exchanger filled with nanofluid by RSM. J Mol Liq. 2016;220:888–901. doi:10.1016/j.molliq.2016.05.031. [Google Scholar] [CrossRef]

28. Greener J, Middleman S. Theoretical and experimental studies of the fluid dynamics of a two-roll coater. Ind Eng Chem Fund. 1979;18(1):35–41. doi:10.1021/i160069a009. [Google Scholar] [CrossRef]

29. Ali F, Hou Y, Feng X, Odeyemi JK, Zahid M, Hussain S. Optimization and sensitivity analysis of heat transfer for Powell-Eyring fluid between rotating rolls with temperature-dependent viscosity: a mathematical modeling approach. Phys Fluids. 2024;36(5):053110. doi:10.1063/5.0211313. [Google Scholar] [CrossRef]

30. Park J, Shin K, Lee C. Roll-to-roll coating technology and its applications: a review. Int J Precis Eng Manuf. 2016;17(4):537–50. doi:10.1007/s12541-016-0067-z. [Google Scholar] [CrossRef]

31. Rehman MIU, Hamid A, Jamshed W, Eid MR. Numerical study of electromagnetically controlled Carreau nanofluid stagnation point flow over a stretchable surface in Darcy-Forchheimer porous medium. Z Angew Math Mech. 2025;105(10):e70243. doi:10.1002/zamm.70243. [Google Scholar] [CrossRef]

32. Ur Rehman MI, Chen H, Hamid A. Multi-physics modeling of magnetohydrodynamic Carreau fluid flow with thermal radiation and Darcy-Forchheimer effects: a study on Soret and Dufour phenomena. J Therm Anal Calorim. 2023;148(24):13883–94. doi:10.1007/s10973-023-12699-9. [Google Scholar] [CrossRef]

33. Montgomery DC. Introduction to statistical quality control. 8th ed. Hoboken, NJ, USA: John Wiley & Sons; 2020. [Google Scholar]