EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5581 ISSN 1307-5543 – ejpam.com Published by New York Business Global Thermal and Concentration Analysis of Mucociliary Johnson-Segalman Fluid Flow in Chronic Obstructive Pulmonary Disease Amel Alaidrous1,∗, Hameed Ashraf2 1 Department of Mathematics, Faculty of Sciences, Umm Al-Qura University, Makkah, Saudi Arabia 2 Department of Mathematics, University of Okara, Okara, Pakistan Abstract. Analyzing mucociliary transport’s thermal and concentration aspects is crucial for com- prehending mucus hypersecretion mechanisms in chronic obstructive pulmonary disease (COPD). The present paper uses a mathematical model of two-layered fluid flow in a finite two-dimensional airway channel. This model uses the hypersecretion scenario, where goblet cells secreted mucus in large volume, leading to complete blockage of the airway channel as the airway ciliary layer (ACL) and peri ciliary layer (PCL). Johnson-Segalman fluid model is used to characterize secreted mucus. The analysis delineated that temperature rises with increasing conductivity ratio, pressure gradient at the airway channel entrance, and Weissenberg number of PCL, while concentration decreases with increasing pressure gradient at the airway channel entrance and Weissenberg num- bers. The analysis also highlights the impact of high Weissenberg numbers on mucus velocity. Johnson-Segalman fluid better characterizes secreted mucus as it improves mucus clearance in the respiratory system, particularly in patients with severe COPD. 2020 Mathematics Subject Classifications: 35Q35, 76A05, 76T06, 92C10 Key Words and Phrases: Chronic obstructive pulmonary disease (COPD), ciliary cells, thermal and concentration, Johnson-Segalman fluid, two-layered model ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5581 Email addresses: aaaidrous@uqu.edu.sa (Amel Alaidrous), hameedashraf09@uo.edu.pk (Hameed Ashraf) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 2 of 30 Nomenclature Dimensional quantities Symbol Interpretation Units a, a1 PCL and ACL mean width mm c Wave speed mm s−1 λ Wavelength mm Û , V̂ Components of velocity that are imparted to the fluid by cilia tips mms−1 ϵ Cilia length mm V̂(k) Velocity vector mm s−1 Û (k), V̂ (k) Velocity component in fixed frame mm s−1 û(k), v̂(k) Velocity component in wave frame mm s−1 P̂ Pressure in fixed frame gm−1s−2 p̂ Pressure in wave frame gm−1s−2 Ŝ(k) Extra stress tensor gmm−1s−2 T0 Temperature at the center of the tube K T1 Temperature at the ciliated wall K T̂ (k) Temperature K C0 Concentration at the center of the tube mol mm−3 C1 Concentration at the ciliated wall mol mm−3 Ĉ(k) Concentration mol mm−3 Ĉk f Specific heat constant Jg−1K−1 χ(k) Thermal conductivity of fluid Wmm−1K−1 D̂k B Coefficient of mass diffusivity mm2s−1 f Body force per unit mass gmm s−2 λ (k) 1 Relaxation time s µ̂(k), η̄(k) Viscosities poise ρ̂k Constant density of fluid g mm−3 X̄0 Reference position of fluid particles mm Dimensionless quantities Symbols Meaning ek Slip parameter We(k) Weissenberg number M Effective viscosity Br(k) Brinkman number SH Schmidt number ST Soret number ζ Conductivity ratio∏ Mass diffusivity ratio ξ Pressure gradient at the airways entrance A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 3 of 30 1. Introduction Secretion of mucus is typically the first line of defense for the human airways. Goblet and ciliary cells are the two main cell types on the mucous membrane of the human airway epithelium. Goblet cells are the specialized secretory cells that secrete a viscoelastic na- ture fluid called mucus. The mucus consists of 95% of water, 2%–3% of macromolecular glycoproteins, 0.1%–0.5% of proteoglycans, and 0.3%–0.5% of lipids, DNA, and proteins. Energy from the hydrolysis of adenosine triphosphate (ATP) opens the pores in a goblet cell’s plasma membrane, triggering the process of exocytosis. The exocytosis, in turn, causes the pores to extend into the airway lumen and secrete mucus quickly. By trapping and clearing inhaled particles, bacteria, and irritants, mucus helps to protect the airways. However, for mucus hypersecretion in chronic obstructive pulmonary disease (COPD), ab- normally high mucus secretion occurs, resulting in excess mucus in the airway lumen. The mucus, in turn, forms the airway surface layer (ASL) and blocks the airway passage. The ASL further comprises two immiscible layers: the upper peri ciliary layer (PCL) and the lower airway ciliary layer (ACL). The mucus buildup may impair lung function, impede mucus clearance, and block airways. This excess of mucus can also trap dust, allergens, and pathogens in the air, increasing respiratory symptoms and raising the risk of infec- tion. To effectively treat airway blockage and enhance respiratory health, it is essential to comprehend the processes driving mucus hypersecretion [13, 15, 28, 29, 31, 34]. The human airways’ mucous membrane epithelium contains ciliary cells found within goblet secretory cells. The structure of ciliary cells is hair-like. A ciliary cell uses the energy from ATP hydrolysis to execute a series of effective and recovery strokes. A cil- iary cell moves the fluid forward during the effective stroke and returns to its initial state in the recovery stroke. A series of effective and recovery strokes cause a ciliary cell’s whip-like swaying motion. A large number of ciliary cells work together to produce a metachronal wave through their collective swaying movements. With varying amplitudes, this wave propagates in the shape of an envelope along an elliptical path at a wave speed [2, 9, 20, 21, 24, 30, 33, 35]. This wave, in turn, makes removing foreign objects like dust, germs, and debris easier, avoiding their buildup and possible damage to the airways. However, ciliary cells may not function properly when there is an increased volume of mucus, leading to impaired mucociliary clearance (MCC) and restricted ciliary function. The excess mucus can hinder the movement of the cilia, causing mucus buildup, airway blockage, and difficulty clearing the mucus [10, 11, 27, 37]. Understanding the role of ciliary cells in hypersecretion is crucial to developing treatments that can improve mucus clearance and alleviate respiratory symptoms. In their investigation, Maqbool et al. [23] considered the effects of mucociliary flow’s concentration and temperature of Phan-Thien-Tanner fluid in an airway channel. The ACL and PCL, two layers of immiscible fluids, were regarded as the components of ASL. The PCL had the highest temperature and velocity, whereas the ACL had the highest concentration. Recently, Ashraf et al. [5] proposed a two-layered mucociliary mathemat- ical model. In their model, the ciliary movement and pressure gradient at the airway channel entrance induced the flow. They characterized the rheological properties of mu- A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 4 of 30 cus in the airways using third-grade fluid. They employed the Adomian decomposition method to solve the formulated nonlinear partial differential equations. They found that the pressure gradient at the channel entrance and the high Deborah number had significant effects. As far as the author knows, previous research has yet to investigate the thermal and concentration aspects of Johnson-Segalman mucociliary fluid flow while considering the hypersecretion scenario of mucus. Further investigation is required to comprehend the flow dynamics related to the mucociliary flow with airway mucus hypersecretion to fill this research void. Most biofluids exhibit viscoelastic properties, such as mucus, vaginal discharge, develop- ing embryos, synovial fluid, and fallopian tubal fluid [8, 14, 19, 26]. Because viscoelastic fluids come in a wide range of varieties, it is still not possible to provide a fluid model that fully captures all their rheological properties. Many models have been developed in the literature to characterize and forecast the behavior and physical attributes of various materials with viscoelastic fluid characteristics. The Johnson-Segalman fluid model is a type of viscoelastic fluid that allows for non-affine deformations. From this model, one can recover both the Maxwell and Newtonian fluids as special cases. The spurt phenomenon can be used to explain this model, which refers to the significant increase in volume until a critical pressure gradient is reached. At this point, there is only a minor increase in the driving pressure gradient [6, 16, 18, 22, 25, 32]. As the mucus moves, it continuously grows, making the Johnson-Segalman fluid model useful for analyzing mucociliary trans- port in the human airways. Mathematics modeling of non-Newtonian fluids produces extremely nonlinear ordinary and partial differential equations. These equations are complex to solve for closed-form solu- tions because of their nonlinear nature. Several analytical techniques have been developed to find approximate series solutions for related issues to overcome this challenge. Differen- tial, ordinary, and partial equations have been solved numerous times with the Adomian decomposition method (ADM). One distinctive feature is that ADM is not constrained by tiny parameters, linearization, or perturbation. Instead, we obtain a convergent solution in a straightforward manner by solving recursive relations. Hosseini and Nasabzadeh [17] provided the criteria for the rapid convergence of the solution obtained through ADM. The first few terms of the series solution can be computed to find an accurate approxima- tion solution. In our case, the mucociliary flow of Johnson-Segalman fluid in the human airways makes finding an exact solution challenging. We employ ADM in this analysis to solve a set of nonlinear partial differential equations [1, 4, 12, 36]. This paper aims to improve the mathematical model of two-layered mucociliary transport in human airways, as proposed by Ashraf et al. [5], within the context of human airway mucus hypersecretion in COPD. This analysis is unique as it focuses on the scenario where the secreted fluid fills the airway channel lumen. The mucus comprises two immiscible layers, the ACL and PCL, which can trap dust, allergens, and pathogens under specific air- way entrance pressure gradient conditions. We characterize the secreted mucus using the Johnson-Segalman fluid model. This analysis focuses on investigating the response of the Johnson-Segalman fluid to concentration and temperature changes in a finite, symmetric airway channel in a mucus hypersecretion context. We will utilize the smaller Reynolds A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 5 of 30 number and long wavelength assumptions to simplify the resultant PDEs. In turn, we will make use of the ADM to solve subsequent PDEs for velocity, pressure gradient, tem- perature, and concentration up to second order. We will show a graphic representation of the temperature, velocity, and concentration distribution to acquire a thorough un- derstanding of mucociliary transport in the case of mucus hypersecretion. By examining the corresponding graphs, we will explore how various parameters influence the ACL and PCL velocities, temperatures, and concentrations. Additionally, we will observe how high Weisenberg numbers affect the velocities of both mucus layers. We will also compare the Johnson-Segalman and Newtonian fluids for velocities, temperatures, and concentrations. 2. Mathematical Modelling 2.1. Formulation of the Problem We consider the two-layered mathematical model of mucociliary two-dimensional Johnson- Segalman fluid flow in a symmetric airway channel of finite length. We improve the math- ematical model proposed formally by Ashraf et al. [5] in the context of mucus hypersecre- tion in chronic obstructive pulmonary disease (COPD). We hypothetically assume that the inner surface of the channel lining the membrane of mucus is populated densely with goblet secretory and ciliary cells. The ciliary cells’ swaying movements, in turn, generate a metachronal wave, which, along with the airway entrance pressure gradient, induces the flow. In addition, we consider a scenario in which the goblet secretory cells poured out a large volume of mucus. The secreted mucus fills the airway channel lumen and forms the airway ciliary layer (ACL) and periciliary layer (PCL). The flow of ACL and PCL occurs in immiscible fluids with differing viscosities, densities, diffusion parameters, and thermal conductivities. Fig. 1 illustrates a schematic representation of the symmetric airway channel. In this rep- resentation, we place the X̄-axis at the midplane of the airways channel, and the Ȳ -axis is perpendicular to it. The mean half-width of the PCL is denoted by a, and the mean half-width of the ACL is denoted by a1. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 6 of 30 Figure 1: Schematic description of the model. This type of flow is described by the velocity vector of the form: V̂(k) = [ Û (k)(X,Y , t), V̂ (k)(X,Y , t), 0 ] , k = A,P. (1) Also, we obtain Ŝ(k) = Ŝ(k)(X,Y , t), (2) T̂ (k) = T̂ (k)(X,Y , t), (3) Ĉ(k) = Ĉ(k)(X,Y , t). (4) Here V̂(k) depict the velocity vector, Ĉ(k) represents the concentration, and T̂ (k) signifies the temperature. We use respectively the k = A and k = P in the superscript to repre- sent the ACL and PCL. The extra stress tensor is represented by Ŝ(k), and the velocity components in the X and Y directions are represented by Û (k) and V̂ (k), respectively. The propagating periciliary metachronal wave H is defined as Y = H = f(X, t) = ± [ a+ aϵ cos( 2π λ (X − ct)) ] , (5) in which ϵ indicates cilia length, λ indicates wavelength, and c denotes the speed of wave. The periciliary metachronal wave and the interface between the periciliary and airway ciliary layers both adhere to the same physical laws, and as a result, the latter exhibits comparable properties. The shape of a wave, a key scientific concept, is determined by the balance of forces acting on a fluid. In a similar manner to the definition of the periciliary metachronal wave, this balance of forces also shapes the interface between the two fluids. Since this metachronal wave travels along the interface, it can be referred to as the airway A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 7 of 30 ciliary metachronal wave H1 [5, 7]: Y1 = H1 = f(X, t) = ± [ a1 + a1ϵ cos( 2π λ (X − ct)) ] . (6) In a symmetric channel, the cilia’s movement follows an elliptical pattern. So, the per- pendicular location of cilia is governed by the following equation: X = g(X, t) = X0 + aϵα sin( 2π λ (X − ct)), (7) here α represents the eccentricity of elliptical path andX0 represents the reference position of the fluid particle. If the no-slip condition is present, the velocities of the fluid particles are determined solely by the cilia tips. Consequently, the components of velocity X and Y can be calculated as follows: Ũ = ∂tX|X=X0 = ∂tg + ∂Xg.∂tX = ∂tg + ∂XgŨ . (8) Ṽ = ∂tY |X=X0 = ∂t̄f + ∂Xf.∂tX = ∂tf + ∂XfŨ . (9) When Eq. (5) and Eq. (7) are used into Eq. (8) and Eq. (9), the following components of velocity imparted by the cilia tip to the fluid particles are obtained: Ũ = −2π λ [ acϵα cos(2πλ (X − ct)) ] 1− 2π λ [ aϵα cos(2πλ (X − ct)) ] . (10) Ṽ = −2π λ [ acϵα sin(2πλ (X − ct)) ] 1− 2π λ [ aϵα cos(2πλ (X − ct)) ] . (11) 2.2. Governing Equations The continuity equation, extra stress, linear momentum conservation, energy equation, and concentration all govern the flow situation under investigation [3, 5, 6, 18, 23, 32]. We introduce and apply these governing equations in the following order: • Equation of continuity: The equation of continuity guarantees the conservation of mass throughout the flow. The continuity equation for an incompressible fluid is stated as follows: ∇ · V̂(k) = 0. (12) The velocity vector (1) identically satisfies the equation of continuity (12). • Extra stress tensor: For the incompressible Johnson-Segalman fluid, we define the extra stress tensor to account for viscoelastic effects and spurt phenomena in the following way: Ŝ(k) = 2µ̂kD̂ (k) + σ̂(k), (13) σ̂(k) + λ (k) 1 [ D̂tσ̂ (k) + σ̂(k) ( Ŵ(k) − ekD̂ (k) ) + ( Ŵ(k) − ekD̂ (k) )Tr σ̂(k) ] = 2η̂kσ̂ (k), (14) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 8 of 30 where the implicit nonlinear viscoelastic behavior and spurt phenomena of the fluid is taken into consideration by σ̂(k), the viscosities are µ̂k and η̂k, the relaxation time is λ (k) 1 , and the slip parameter is ek. The symmetric and skew-symmetric components of the velocity gradient denoted by D̂(k) and Ŵ(k), respectively are defined as D̂(k) = 1 2 [ gradV̂(k) + ( gradV̂(k) )T ] , Ŵ(k) = 1 2 [ gradV̂(k) − ( gradV̂(k) )Tr ] , (15) where we represent the gradient operator with grad and transpose with Tr. Remark :- In Eq. (14), setting ek = 1 and η̂k = 0 recovers the extra stress tensor for Maxwell fluid whereas setting η̂k = λ (k) 1 = 0 yields the extra stress tensor for Newtonian fluid. • Conservation of Linear Momentum: The relationship between the forces acting on the fluid and its momentum change rate allows the momentum balance to describe the behavior of the incompressible Johnson-Segalman fluid. We define the balance of momentum as follows: ρ̂(k)DtV̂ (k) = ∇ · Ŝ(k) −∇P̂ + ρ̂(k)f̂ (k), (16) where f̂ (k) represents body force per unit mass, P̂ represents pressure, D Dt represents material time derivative, and ρ̂(k) represents constant fluid densities. • Energy Equation: The transfer of heat within the human airways is governed by the energy equation, which defines how heat is transferred within the fluid. The energy equation is defined as ρ̂(k)Ĉ (k) f DtT̂ (k) = χ(k)∇2T̂ (k) + tr(Ŝ(k).∇V̂(k)), (17) where the thermal conductivity as χ(k), the specific heat constant as Ĉ (k) f , the constant ratio due to thermal diffusion as DKT , and the coefficient of mass diffusivity as D̂ (k) B . • Concentration Equation: The concentration equation accurately describes the distri- bution and variation of substances within respiratory mucus. We define the concentration equation in the following way DtĈ (k) = D̂ (k) B ∇2Ĉ(k) + DKT T̂0 ∇2T̂ (k). (18) Here, Ĉ (k) f denotes the specific heat constant, DKT represents the constant ratio due to thermal diffusion, and D̂ (k) B is used to denote the coefficient of mass diffusivity. 2.3. Dimensionless Coordinate Transformations We introduce a wave reference frame (x, y) to facilitate the analysis of mucociliary flow. Since this frame propagates at wave speed c, therefore, we can investigate steady behavior A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 9 of 30 of the mucociliary flow. Through the dimensionless transformations the coordinates in the wave reference frame (x, y) and the fixed frame of reference (X, Y , t) are now related to transform the velocity components in the fixed frame of reference (Û (k), V̂ (k)) to the wave reference frame (u(k), v(k)), pressure in the fixed frame of reference (P̂ ) to the wave reference frame (p), and extra stress tensor (Ŝ(k)) from fixed frame of reference to the wave reference frame (S(k)). Following dimensionless transformations we scaled for this purpose: x = X − ct λ , y = Y a , u(k) = Û (k) − c c , v(k) = V̂ (k)λ ca , ρ(k) = ρ̂(k) ρ̂(A) , δ = a1 a , h1 = Ĥ1 a , h = Ĥ a , p = P̂ a2 (µ̂A + η̂A) cλ , S(k) = Ŝ(k)a cµ̂A , A (k) 1 = Â (k) 1 a c , T (k) = T̂ (k) − T0 T1 − T0 , C(k) = Ĉ(k) − C0 C1 − C0 . (19) We can more effectively investigate the mucociliary flow within the reference wave frame and simplify the study owing to these transformations. Eqs. (5), (6) and (12)-(18) follow- ing the application of dimensionless transformations (19), utilizing the smaller Reynolds number (Re → 0) and long wavelength (β → 0) assumptions subsequently yields h(x) = ±[1 + ϵ cosx], (20) h1(x) = ±[δ + δϵ cosx], (21) ∂xu (k) + ∂yv (k) = 0, (22) S(k) xy = 1 M [ ∂yu (k) + ( We(k) )2 ( 1− e2k ) M ( ∂yu (k) )3 1 + ( We(k) )2 ( 1− e2k ) ( ∂yu(k) )2 ] , (23) M∂yS (k) xy = dxp, (24) ∂yyT (k) = −Br(k)S(k) yx ∂yu (k), (25) ∂yyC (k) = −SHST D (k) B ∂yyT (k), (26) in which We(k) = λ (k) 1 c a is used to represent Weissenberg number, M = µ̂A µ̂A + η̂A is used to represent effective viscosity, Br(k) = c2µ̂A χ(k) (T1 − T0) is used to represent Brinkman number, SH = µ̂A D̂ (A) B ρ̂A is used to represent Schmidt number and ST = ρ̂ADKT µ̂A(C1 − C0) is used to represent Soret number. In biological flows, the maximum velocity, temperature, and concentration occur at the A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 10 of 30 midplane of the airway channel. Shear forces, fluid velocities, temperatures, concentra- tions, and temperature and concentration gradients are the same at the ACL-PCL interface [5, 23]. Consequently, we have the following boundary conditions in dimensionless wave frame of reference: S(A) xy (x, 0) = 0, (27) v(A)(x, 0) = 0, (28) ∂yT (A)(x, 0) = 0, (29) ∂yC (A)(x, 0) = 0. (30) S(A) xy (x, h1) = S(P ) xy (x, h1), (31) u(A) xy (x, h1) = u(P ) xy (x, h1), (32) v(A)(x, h1) = v(P )(x, h1), (33) T (A)(x, h1) = T (P )(x, h1), (34) ∂yT (A)(x, h1) = ζ∂yT P (x, h1), (35) ∂yC (A)(x, h1) = Π∂yC (P )(x, h1), (36) C(A)(x, h1) = C(P )(x, h1). (37) T (2)(x, h) = 0, (38) C(2)(x, h) = 0, (39) u(2)(x, h) = −1− 2πϵαβ cos(2πx), (40) v(2)(x, h) = ±2πϵ(sin(2πx) + 2πϵαβ sin(2πx)cos(2πx)). (41) In Eq. (35) and Eq. (36), ζ = χ(P ) χ(A) denotes for conductivity ratio and Π = D (A) B D (P ) B denotes for mass diffusivity ratio. The condition of the pressure gradient at the airway channel entrance is described as [5]: dxp(0) = −ξ. (42) After integrating Eq. (24) concerning ‘y’ and invoking the boundary conditions (27) and (31) into resulting equation, one finally arrives at MS(k) xy = ydxp. (43) Putting Eq. (23) into Eq. (43) and Eq. (25), respectively we have ∂yu (k) = G (k) 1 [( ∂yu (k) )2 ydxp−M ( ∂yu (k) )3 ] + ydxp, (44) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 11 of 30 ∂yyT (k) = −G (k) 2 [( ∂yu (k) )2 +G (k) 3 ( ∂yu (k) )4 −G (k) 4 ( ∂yu (k) )6 ] . (45) Substituting Eq. (45) into Eq. (26), we get ∂yyC (k) = J (k) 1 [( ∂yu (k) )2 +G (k) 3 ( ∂yu (k) )4 −G (k) 4 ( ∂yu (k) )6 ] . (46) In this section, we have derived three nonlinear partial differential equations, namely Eq. (45), Eq. (46), and Eq. (44). The first two are second-order, while the third one is first- order. All these equations seem to be unsolvable by exact methods. In the next section, we will look for upto second-order solution using the ADM while taking into account the boundary conditions (28)-(30) and (32)-(42). 3. Solution of the Problem We can write Eqs. (44)-(46) in operator form by denoting the nonlinear terms: ( ∂yu (k) )2 ,( ∂yu (k) )3 , ( ∂yu (k) )4 and ( ∂yu (k) )6 with N1(u (k)), N2(u (k)), N3(u (k)) and N4(u (k)) respec- tively and defining one fold Ly = ∂y and two fold Lyy = ∂yy linear invertible operators as shown below: Lyu (k)(x, y) = G (k) 1 [ N1(u (k))ydxp−MN2(u (k)) ] + ydxp, (47) LyyT (k)(x, y) = −G (k) 2 [ N1(u (k)) +G (k) 3 N3(u (k))−G (k) 4 N4(u (k)) ] , (48) LyyC (k)(x, y) = J (k) 1 [ N1(u (k)) +G (k) 3 N3(u (k))−G (k) 4 N4(u (k)) ] . (49) We define respectively the inverse operators of one Ly and two Lyy folds operators as follows L−1 y (⋆) = ∫ (⋆)dy, (50) L−1 yy (⋆) = ∫ ∫ (⋆)dydy. (51) Following the application of the inverse operators (50) and (51) respectively to Eq. (47) and Eqs. (48) and (49), one acquires u(k)(x, y) = G (k) 1 [ dxpL −1 y N1(u (k))y −ML−1 y N2(u (k)) ] + y2 2 dxp+D (k) 1 (x), (52) T (k)(x, y) = −G (k) 2 [ L−1 yy N1(u (k)) +G (k) 3 L−1 yy N3(u (k))−G (k) 4 L−1 yy N4(u (k)) ] +D (k) 2 (x)y+D (k) 3 (x), (53) C(k)(x, y) = J (k) 1 [ L−1 yy N1(u (k)) +G (k) 3 L−1 yy N3(u (k))−G (k) 4 L−1 yy N4(u (k)) ] +D (k) 4 (x)y+D (k) 5 (x). (54) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 12 of 30 Here, D (k) 1 (x), D (k) 2 (x), D (k) 3 (x), D (k) 4 (x) and D (k) 5 (x) represent the arbitrary functions. The uk(x, y), vk(x, y), dxp, T k(x, y), and Ck(x, y) are now decomposed as follows: u(k)(x, y) = ∞∑ i=0 u (k) i (x, y), (55) v(k)(x, y) = ∞∑ i=0 v (k) i (x, y), (56) dxp = ∞∑ i=0 dxpi, (57) T (k)(x, y) = ∞∑ i=0 T (k) i (x, y), (58) C(k)(x, y) = ∞∑ i=0 C (k) i (x, y). (59) The nonlinear terms N1(u (k)), N2(u (k)), N3(u (k)), and N4(u (k)) are now expanded with the following form of infinite series of Adomian polynomials: N1(u (k)) = ∞∑ i=0 A (k) i , N2(u (k)) = ∞∑ i=0 B (k) i , N3(u (k)) = ∞∑ i=0 E (k) i , N4(u (k)) = ∞∑ i=0 F (k) i . (60) After inserting the Adomian polynomials (60) and the series of decomposition (55)-(59) into Eqs. (52)-(54), we finally get u (k) 0 (x, y)+ ∞∑ i=1 u (k) i (x, y) = D (k) 1 (x)+ y2 2 dxp0+G (k) 1 [ dxpL −1 y ∞∑ i=0 A (k) i y −ML−1 y ∞∑ i=0 B (k) i ] , (61) T (k) 0 (y) + ∞∑ i=1 T (k) i (x, y) = D (k) 2 (x)y +D (k) 3 (x)−G (k) 2 [ L−1 yy ∞∑ i=0 A (k) i +G (k) 3 L−1 yy ∞∑ i=0 E (k) i −G (k) 4 L−1 yy ∞∑ i=0 F (k) i ] , (62) C (k) 0 (x, y) + ∞∑ i=1 C (k) i (x, y) = D (k) 4 (x)y +D (k) 5 (x) + J (k) 1 [ L−1 yy ∞∑ i=0 A (k) i +G (k) 3 L−1 yy ∞∑ i=0 E (k) i −G (k) 4 L−1 yy ∞∑ i=0 F (k) i ] . (63) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 13 of 30 Eqs. (61)-(63) respectively give rise to the component of zeroth order and a recursive relations of the form: u (k) 0 (x, y) = D (k) 1 (x) + y2 2 dxp0, (64) ∞∑ i=1 u (k) i (x, y) = G (k) 1 [ dxpL −1 y ∞∑ i=0 A (k) i y −ML−1 y ∞∑ i=0 B (k) i ] , (65) T (k) 0 (x, y) = D (k) 2 (x)y +D (k) 3 (x), (66) ∞∑ i=1 T (k) i (x, y) = −G (k) 2 [ L−1 yy ∞∑ i=0 A (k) i +G (k) 3 L−1 yy ∞∑ i=0 E (k) i −G (k) 4 L−1 yy ∞∑ i=0 F (k) i ] , (67) C (k) 0 (x, y) = D (k) 4 (x)y +D (k) 5 (x), (68) ∞∑ i=1 C (k) i (x, y) = J (k) 1 [ L−1 yy ∞∑ i=0 A (k) i +G (k) 3 L−1 yy ∞∑ i=0 E (k) i −G (k) 4 L−1 yy ∞∑ i=0 F (k) i ] . (69) After making use the decomposition series (55)-(59) on Adomian polynomials (60), and in turn using binomial expansion, we get A (k) 0 = ( ∂yu (k) 0 )2 , A (k) 1 = 2 ( ∂yu (k) 0 )( ∂yu (k) 1 ) , (70) B (k) 0 = ( ∂yu (k) 0 )3 , B (k) 1 = 3 ( ∂yu (k) 0 )2 ( ∂yu (k) 1 ) , (71) E (k) 0 = ( ∂yu (k) 0 )4 , E (k) 1 = 4 ( ∂yu (k) 0 )3 ( ∂yu (k) 1 ) , (72) F (k) 0 = ( ∂yu (k) 0 )6 , F (k) 1 = 6 ( ∂yu (k) 0 )5 ( ∂yu (k) 1 ) . (73) Invoking the decomposition series (55)-(59) into the boundary conditions (28)-(30) and (32)-(42) lead to the following results: u (A) 0 (x, h1) = u (P ) 0 (x, h1), u (A) 1 (x, h1) = u (P ) 1 (x, h1), u (A) 2 (x, h1) = u (P ) 2 (x, h1), (74) u (P ) 0 (x, h) = −1− 2πϵαβ cos(2πx), u (P ) 1 (x, h) = 0, u (P ) 2 (x, h) = 0, (75) v (A) 0 (x, 0) = 0, v (A) 1 (x, 0) = 0, v (A) 2 (x, 0) = 0, (76) v (P ) 0 (x, h) = ±2πϵ(sin(2πx) + 2πϵαβ sin(2πx)cos(2πx)), v (P ) 1 (x, h) = 0, v (P ) 2 (x, h) = 0, (77) v (A) 0 (x, h1) = v (P ) 0 (x, h1), v (A) 1 (x, h1) = v (P ) 1 (x, h1), v (A) 2 (x, h1) = v (P ) 2 (x, h1), (78) ∂yT (A) 0 (x, 0) = 0, ∂yT (A) 1 (x, 0) = 0, ∂yT (P ) 2 (x, 0) = 0, (79) T (A) 0 (x, h1) = T (P ) 0 (x, h1), T (A) 1 (x, h1) = T (P ) 1 (x, h1), T (A) 2 (x, h1) = T (P ) 2 (x, h1), (80) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 14 of 30 ∂yT (A) 0 (x, h1) = ζ∂yT (P ) 0 (x, h1), ∂yT (A) 1 (x, h1) = ζ∂yT (P ) 1 (x, h1), ∂yT (A) 2 (x, h1) = ζ∂yT (P ) 2 (x, h1), (81) T (P ) 0 (x, h) = 0, T (P ) 1 (x, h) = 0, T (P ) 2 (x, h) = 0, (82) ∂yC (A) 0 (x, 0) = 0, ∂yC (A) 1 (x, 0) = 0, ∂yC (A) 2 (x, 0) = 0, (83) C (A) 0 (x, h1) = C (P ) 0 (x, h1), C (A) 1 (x, h1) = C (P ) 1 (x, h1), C (A) 2 (x, h1) = C (P ) 2 (x, h1), (84) ∂yC (A) 0 (x, h1) = Π∂yC (P ) 0 (x, h1), ∂yC (A) 1 (x, h1) = Π∂yC (P ) 1 (x, h1), ∂yC (A) 2 (x, h1) = Π∂yC (P ) 2 (x, h1), (85) C (P ) 0 (x, h) = 0, C (P ) 1 (x, h) = 0, C (P ) 2 (x, h) = 0, (86) dxp0(0) = −ξ, dxp1(0) = 0, dxp2(0) = 0. (87) 3.1. Velocity After solving Eq. (64) and Eq. (65) using the concern Adomian polynomials from Eq. (60), in turn utilizing appropriate boundary conditions given in (74) and (75) into the resulting equations, we respectively acquire the zeroth, first, and second order solutions of both ACL and PCL as follows: u (A) 0 (x, y) = −1− 2πϵαβ cos(2πx) + 1 2 dxp0 ( y2 − h2 ) , (88) u (A) 1 (x, y) = 1 4 (dxp0) 3 [Γ1 ( y4 − h41 ) + 4Γ2 ( h41 − h4 )] + 1 2 dxp1 ( y2 − h2 ) , (89) u (A) 2 (x, y) = 1 6 (dxp0) 4 (2dxp1 − 3Mdxp0) [ Γ3 ( y6 − h61 ) + Γ4 ( h61 − h6 )] + 1 4 (dxp0)(dxp1) × (2dxp1 + 3dxp0) [ G (A) 1 ( y4 − h41 ) +G (P ) 1 ( h41 − h4 )] + 1 2 dxp2 ( y2 − h2 ) . (90) u (P ) 0 (x, y) = −1− 2πϵαβ cos(2πx) + 1 2 dxp0 ( y2 − h2 ) , (91) u (P ) 1 (x, y) = Γ2 4 (dxp0) 3 (y4 − h4 ) + 1 2 dxp1 ( y2 − h2 ) , (92) u (P ) 2 (x, y) = Γ4 6 (dxp0) 4 (2dxp1 − 3Mdxp0) ( y6 − h6 ) + 1 4 (dxp0)(dxp1) (2dxp1 + 3dxp0) ( y4 − h4 ) + 1 2 dxp2 ( y2 − h2 ) . (93) Substitution of Eqs. (88)-(90) into Eq. (55) for k = A and Eqs. (91)-(93) into Eq. (55) for k = P yield respectively the ACL and PCL velocities. Using decomposition series (56) into Eq. (22), in turn the resulting equation for k = A upon using Eqs. (88)-(90) and for k = P upon using Eqs. (91)-(93) along with associated A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 15 of 30 boundary conditions (76) and (77), we respectively obtain the zeroth, first, and second order solutions of both ACL and PCL as follows: v (A) 0 (x, y) = −4π2ϵαβ sin(2πx)y − 1 6 dxxp0 ( y3 − 3h2y ) + dxp0hdxh, (94) v (A) 1 (x, y) = − 3 20 (dxp0) 2 dxxp0 [ Γ1 ( y5 − 5h41y ) + 4Γ2 ( h41 − h4 ) y ] + (dxp0) 3 [Γ1h 3 1dxh1 − 4Γ2 ( h31dxh1 − h3dxh )] y − 1 6 dxxp1 ( y3 − 3h2y ) + dxp1hdxh, (95) v (A) 2 (x, y) = − 2 21 (dxp0) 3 ((2dxp1 − 3Mdxp0) dxxp0 + 4dxp0 (2dxxp1 − 3Mdxxp0)) [ Γ3 ( y7 − 7h61y ) + 7Γ4 ( h61 − h6 ) y ] + (dxp0) 4 (2dxp1 − 3Mdxp0) [ Γ3h 5 1dxh1 − Γ4 ( h51dxh1 − h5dxh )] y − 1 4 ( 4dxp0dxp1dxxp1 + 2dxxp0 (dxp1) 2 + 6dxp0dxxp0dxp1 + 3 (dxp0) 2 dxxp1 ) [ G (A) 1 × ( y5 − 5h41y ) + 5G (P ) 1 ( h41 − h4 ) y ] + ( 2dxp0 (dxp1) 2 + 3 (dxp0) 2 dxp1 ) [ G (A) 1 h31dxh1 −4G (P ) 1 ( h31dxh1 − h3dxh )] y − 1 6 dxxp2 ( y3 − 3h2y ) + dxp2hdxh. (96) v (P ) 0 (x, y) = ±2πϵ(sin(2πx) + 2πϵαβ sin(2πx)cos(2πx))− 4π2ϵαβ sin(2πx)(y − h) − 1 6 dxxp0 ( y3 − 3h2y ) + dxp0hdxh, (97) v (P ) 1 (x, y) = −Γ2 20 [ 3 (dxp0) 2 dxxp0 ( y5 − 5h4y ) + 20 (dxp0) 3 h3dxhy ] − 1 6 dxxp1 ( y3 − 3h2y ) + dxp1hdxh, (98) v (P ) 2 (x, y) = −Γ4 42 [ (dxp0) 3 ((2dxp1 − 3Mdxp0) dxxp0 + 4dxp0 (2dxxp1 − 3Mdxxp0)) ( y7 − 7h6y ) − 42 (dxp0) 4 (2dxp1 − 3Mdxp0)h 5dxhy ] − 1 20 [( 4dxp0dxp1dxxp1 + 2dxxp0 (dxp1) 2 + 6dxp0dxxp0dxp1 + 3 (dxp0) 2 dxxp1 ) ( y5 − 5h4y ) − 20 ( 2dxp0 (dxp1) 2 + 3 (dxp0) 2 dxp1 ) × h3dxhy ] − 1 6 dxxp2 ( y3 − 3h2y ) + dxp2hdxh. (99) When we make use of Eqs. (94)-(96) into Eq. (56) for k = A and Eqs. (97)-(99) into the Eq. (56) for k = P respectively acquire the ACL and PCL velocities. 3.2. Pressure Gradient The following components of zeroth, first, and second-order pressure gradients are derived by solving Eqs. (94)-(99) and applying the related boundary conditions (78): dxp0 = − 6(1 + 2πϵαβcos(2πx))h+ 3 ( 2ϵcos(2πx) + πϵ2αβcos(4πx) ) + 6 (Γ5 − Γ6) 3 ( h31 − h21 + h2h1 ) + 2h3 , (100) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 16 of 30 dxp1 = 6Λ (dxp0) 3 20(3h31 − 5h3) [ 4G (A) 1 ( Γ9 − h51 ) +G (P ) 1 ( h51 − h4h1 − Γ8 ) +G (P ) 1 ( Γ7 − h51 + 4h5 − 5h1h 2 )] , (101) dxp2 = 6 3h21 − 3h2 − 7h31 + 3h2h1 + 2h3 [ Λ 42 (dxp0) 4 (2dxp1 − 3Mdxp0) {( G (P ) 1 )2 ( h71 − 7h6h1 + 13h7 − 7h61h− Γ13 − Γ8 ) − 7 ( G (A) 1 )2 ( h71 + Γ10 )} + dxp0 (dxp1) 2 { G (P ) 1 ( h51 − 5h4h1 − 4h5 − 4h41 + 4h4 ) − 20G (A) 1 ( h51 − Γ7 )} − dxp0dxp1 G (P ) 1 20 { 2dxp1 ( h51 − 5h5h1 − 4h5 ) × (3dxp0 − 2dxp1) Γ5}] . (102) By inserting the pressure gradient components from Eqs. (100)-(102) into Eq. (57), we can get the pressure gradient dxp. 3.3. Temperature Substitution of Adomian polynomials (70), (72) and (73) into Eqs. (66) and (67). After applying the correct boundary conditions (79)-(82) to the resultant equations, the com- ponents of zeroth, first, and second-order are derived when we set k = P followed by k = A: T (A) 0 (x, y) = 0, (103) T (A) 1 (x, y) = G (A) 2 840 [ 15G (A) 4 (dxp0) 6 (y8 − h81 ) − 28G (A) 3 (dxp0) 4 (y6 − h61 ) − 70 (dxp0) 2 (y4 − h41 )] + T (P ) 1 (x, h1) , (104) T (A) 2 (x, y) = −G (A) 2 840 [ 56G (A) 1 Λ (dxp0) 4 (y6 − h61 ) + 140dxp0dxp1 ( y4 − h4 ) + 4G (A) 3 { 15G (A) 1 (dxp0) 6 × Λ ( y8 − h81 ) + 28 (dxp0) 3 dxp1 ( y6 − h61 )} − 7G (A) 4 { 8G (A) 1 Λ (dxp0) 8 (y10 − h101 ) + 15 (dxp0) 5 dxp1 ( y8 − h81 )}] + T (2) 2 (x, h1). (105) T (P ) 0 (x, y) = 0, (106) T (P ) 1 (x, y) = G (P ) 2 840 [ 15G (P ) 4 (dxp0) 6 (y7 − h7 ) − 28G (P ) 3 (dxp0) 4 (y5 − h5 ) − 70 (dxp0) 2 (y4 − h4 )] + G (A) 2 105ζ [ 15G (A) 4 (dxp0) 6 h71 − 21G (A) 3 (dxp0) 4 h51 − 35 (dxp0) 2 h31 ] (y − h)− G (P ) 2 105 [ 15G (P ) 4 × (dxp0) 6 h71 + 21G (P ) 3 (dxp0) 4 h51 + 35 (dxp0) 2 h31 ] (y − h), (107) T (P ) 2 (x, y) = G (P ) 2 840 [ 56G (P ) 1 (dxp0) 2 Λ ( y6 − h6 ) + 140dxp0dxp1 ( y4 − h4 ) + 8G (P ) 3 { 15G (P ) 1 (dxp0) 6 A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 17 of 30 × Λ ( y8 − h8 ) + 28 (dxp0) 3 dxp1 ( y6 − h6 )} − 7G (P ) 4 { 8G (P ) 1 (dxp0) 8 Λ ( y10 − h10 ) + 15 (dxp0) 5 dxp1 ( y8 − h8 )}] − G (P ) 2 105 [ 42G (P ) 1 (dxp0) 4 h51Λ + 35dxp0dxp1h 3 1 + G (P ) 3 { 20G (P ) 1 (dxp0) 6 h71Λ + 28 (dxp0) 3 dxp1h 5 1 } −G (P ) 4 { 2G (P ) 1 (dxp0) 8 h41Λ + 5 (dxp0) 5 dxp1h 7 1 }] (y − h) + G (A) 2 105ζ [ 42G (A) 1 (dxp0) 4 h51Λ + 70dxp0dxp1h 3 1 + 3G (A) 3 { 20G (A) 1 (dxp0) 6 h71Λ + 28 (dxp0) 3 dxp1h 6 1 } − 35G (A) 4 { 2G (A) 1 (dxp0) 8 h41Λ + 3 (dxp0) 5 dxp1h 7 1 }] (y − h). (108) When we make use of temperature components calculated in Eqs, (103)-(105) into Eq. (58) for k = A and Eqs. (106)-(108) into Eq. (58) for k = P , we get the final ACL and PCL temperature expressions. 3.4. Concentration Using Adomian polynomials (70), (72), and (73) is the first step towards solving Eqs. (68) and (69). By applying the boundary conditions of concern (83)-(87) to the resultant equations, we obtain the concentrations components of zeroth, first, and second-order when we set respectively the k = P followed by k = A: C (A) 0 (x, y) = 0, (109) C (A) 1 (x, y) = − J (A) 1 840G (A) 5 [ 15G (A) 4 (dxp0) 6 (y8 − h81 ) − 28G (A) 3 (dxp0) 4 (y6 − h61 ) − 70 (dxp0) 2 (y4 − h41 )] + C (P ) 1 (x, h1) , (110) C (A) 2 (x, y) = J (A) 1 840G (A) 5 [ 56G (A) 1 Λ (dxp0) 4 (y6 − h61 ) + 140dxp0dxp1 ( y4 − h4 ) + 4G (A) 3 { 15G (A) 1 (dxp0) 6 × (1−M) ( y8 − h81 ) + 28 (dxp0) 3 dxp1 ( y6 − h61 )} − 7G (A) 4 { 8G (A) 1 (1−M) (dxp0) 8 × ( y10 − h101 ) + 15 (dxp0) 5 dxp1 ( y8 − h81 )}] + C (P ) 2 (x, h1). (111) C (P ) 0 (x, y) = 0, (112) C (P ) 1 (x, y) = − J (P ) 1 840G (P ) 5 [ 15G (P ) 4 (dxp0) 6 (y7 − h7 ) − 28G (P ) 3 (dxp0) 4 (y5 − h5 ) − 70 (dxp0) 2 (y4 − h4 )] − J (A) 1 105ζG (A) 5 [ 15G (A) 4 (dxp0) 6 h71 − 21G (A) 3 (dxp0) 4 h51 − 35 (dxp0) 2 h31 ] (y − h) + J (P ) 1 105G (P ) 5 × [ 15G (P ) 4 (dxp0) 6 h71 + 21G (P ) 3 (dxp0) 4 h51 + 35 (dxp0) 2 h31 ] (y − h), (113) A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 18 of 30 C (P ) 2 (x, y) = − J (P ) 1 840G (P ) 5 [ 56G (P ) 1 (dxp0) 2 Λ ( y6 − h6 ) + 140dxp0dxp1 ( y4 − h4 ) + 8G (P ) 3 { 15G (P ) 1 (dxp0) 6 × Λ ( y8 − h8 ) + 28 (dxp0) 3 dxp1 ( y6 − h6 )} − 7G (P ) 4 { 8G (P ) 1 (dxp0) 8 Λ ( y10 − h10 ) + 15 (dxp0) 5 dxp1 ( y8 − h8 )}] + J (P ) 1 105G (P ) 5 [ 42G (P ) 1 (dxp0) 4 h51Λ + 35dxp0dxp1h 3 1 + G (P ) 3 { 20G (P ) 1 (dxp0) 6 h71Λ + 28 (dxp0) 3 dxp1h 5 1 } −G (P ) 4 { 2G (P ) 1 (dxp0) 8 h41Λ + 5 (dxp0) 5 dxp1h 7 1 }] (y − h)− J (A) 1 105ζG (A) 5 [ 42G (A) 1 (dxp0) 4 h51Λ + 70dxp0dxp1h 3 1 + 3G (A) 3 { 20G (A) 1 (dxp0) 6 h71Λ + 28 (dxp0) 3 dxp1h 6 1 } − 35G (A) 4 { 2G (A) 1 (dxp0) 8 h41Λ + 3 (dxp0) 5 dxp1h 7 1 }] (y − h). (114) Putting the concentration components from Eqs. (109)-(111) into the Eq. (59) for k = A and from Eqs. (112)-(114) into the (59) for k = P , one gets the final expression for con- centration. Remark : By putting the We = 0 into Eqs. (88)-(114), we can obtain the expressions for Newtonian mucus velocity components, pressure gradient, temperatures, and concentra- tions. 4. Results and Discussion A mathematical model incorporating a finite narrow two-dimensional channel comprising two layers (ACL and PCL), where the Johnson-Segalman fluid model characterizes the mucus rheological properties secreted within the airways, was considered [5]. The ADM and the theory of lubrication approximation have been employed to obtain the series form up to the second-order solution of the resultant simplified set of non-linear PDEs. The pri- mary focus of this section is to provide an estimation of the quantitative variation impact of slip parameters (eA and eP ), Weissenberg numbers ((We)(A) and (We)(P )), Brinkmann numbers (Br(A) and Br(P )), pressure gradient at the airway channel entrance ξ, conduc- tivity ratio ζ, mass diffusivity ratio Π, Schmidt number SH , and Soret number ST involved in the present analysis. Graphs are displayed and discussion is made to seek the impacts of aforementioned involved parameters on the velocity (u versus y plots), temperature (T versus y plots) and concentration (C versus y plots) (see Figs. 2-9). Variation impact of high Weissenberg numbers ((We)(A) and (We)(P )) on velocity is demonstrated and analyzed (see Table 1 and Table 2). 4.1. Velocity We constructed the graphs shown in Figs. 2 and 3 to examine the effects of variation on eA, eP , ξ, (We)(A), and (We)(P ) on mucociliary velocity. These graphs show plots of A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 19 of 30 velocity against normal distance. The velocity rises as normal distance does. The effects of eA and eP on the two layers’ velocities are shown in Figs. 2(a) and (b). The velocity and both eA and eP have a direct relationship. In other words, a rise in eA and eP causes the ACL and PCL to undergo a velocity increase. (We)(A) has very little influence on velocity, whereas (We)(P ) has a significant effect. For both layers, prolonged relaxation times are correlated with an increase in (We)(A) and (We)(P ) as the Weissenberg number is directly proportional to the fluid’s relaxation time and inversely proportional to the observation time. The fluid behaves viscoelastically, with the elastic effects becoming more apparent as the relaxation time increases. Furthermore, we can state that both layers’ velocities rise with increased ACL and PCL relaxation times. The variation impacts of ξ on velocity are presented in Fig. 3(a). These graphs show that there is rise in the velocity as the values of ξ increment. The influence of ξ on velocity is apparent. That said, we add that the both layers velocity tend to increase with increasing pressure gradient at the airway entrance. This figure shows that both layers’ velocities increase as ξ values do. When there is a rise in the pressure gradient at the entrance, the fluid flows more swiftly. The Johnson-Segalman and Newtonian fluid velocity comparison is shown in Fig. 3(b). We demonstrated that when mucus is characterized using the Johnson-Segalman fluid, it flows more quickly in the ACL and PCL than the Newtonian fluid. Figure 2: The variation impacts of (a) eA, (b) eP , (c) (We)(A) and (d) (We)(P ) on velocity u. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 20 of 30 Figure 3: The variation impacts of (a) ξ on velocity u and (b) comparison between the Johnson-Segalman fluid and Newtonian fluid for velocity u. When the Weissenberg number of a Johnson-Segalman fluid is low, it exhibits more vis- cous effects, and when it is high, it demonstrates more elastic effects. A high Weissenberg number causes the fluid to flow easily with more velocity, while a low Weissenberg number results in more viscous behavior. This value offers a quantitative analysis that differen- tiates fluids that are viscous and viscoelastic. The effects of a high Weissenberg number, (We)(A), on the ACL and PCL velocities of Johnson-Segalman fluid are displayed numeri- cally in Table 1. These statistics clearly show that when there is an increment in the values from 1.0 to 10.0, the mucociliary velocity reduces and the flow direction inverts. It may assist the fluid in maintaining its structure by exhibiting resistance to deformation. On the other hand, Table 2 demonstrates that as the values of (We)(A) increase from 1.0 to 10.0, velocity increases rapidly, potentially causing the fluid to break down its structure or offer lesser resistance to deformation, behaving more like a liquid. The effects of (We)(A) and (We)(P ) are vice versa on ACL and PCL velocity. Different effects of high Weissenberg numbers are shown in ACL and PCL velocities; (We)(A) influences the direction of flow, whereas (We)(P ) forces the fluid to move more easily and quickly in the forward direction. Table 1: Variation impact of high Weissenberg number (We)(A) on velocity. (We)(A) = 1.00 (We)(A) = 5.00 (We)(A) = 10.00 y u(A)(x, y) u(P )(x, y) u(A)(x, y) u(P )(x, y) u(A)(x, y) u(P )(x, y) 0.0 2.31568 2.39972 11.4225 9.25325 -29.7104 -64.4501 0.2 2.17309 2.25658 10.9828 8.81956 -27.2559 -61.747 0.4 1.75015 1.82635 9.01264 7.51723 -30.3327 -53.6409 0.6 1.05404 1.10653 -0.76465 5.34247 -126.173 -40.141 0.8 0.07215 0.09247 -43.2254 2.2885 -646.53 -21.2634 1.0 -1.26441 -1.22322 -183.468 -1.65501 -2447.23 2.96882 A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 21 of 30 Table 2: Variation impact of high Weissenberg number (We)(P ) on velocity. (We)(P ) = 1.00 (We)(P ) = 5.00 (We)(P ) = 10.00 y u(A)(x, y) u(P )(x, y) u(A)(x, y) u(P )(x, y) u(A)(x, y) u(P )(x, y) 0.0 418.307 418.899 3.96301× 107 3.96334× 107 9.45422× 109 9.45443× 109 0.2 400.517 401.108 3.79475× 107 3.79507× 107 9.05244× 109 9.05265× 109 0.4 347.148 347.712 3.28998× 107 3.29016× 107 7.84711× 109 7.84722× 109 0.6 258.195 258.627 2.44870× 107 2.44826× 107 5.83822× 109 5.83792× 109 0.8 133.653 133.670 1.27089× 107 1.26874× 107 3.02577× 109 3.02435× 109 1.0 -26.488 -27.514 −2.43432× 106 −2.49344× 107 −5.90229× 108 −5.94113× 108 4.2. Temperature Figs. 4 and 5 display the impact of several parameters, namely eA, eP , (We)(A), (We)(P ), Br(A), Br(P ), Π, and ξ, on the ACL and PCL temperature. Figs. 4(a) and 4(b) specif- ically show the impact of eA and eP on temperature. Significantly, these characteristics display contrasting impacts on temperature. An increased value in eA causes a drop in temperature, whereas an elevation in eP leads to a boost in temperature. Similarly, the variations in (We)(A) and (We)(P ) also have different impacts on temperature, as seen in Figs. 4(c) and 4(d). An increase in the (We)(A) leads to a fall in temperature, whereas an increase in the (We)(P ) causes the temperature to increment. The temperature variations can be ascribed to the relaxation time. Changes in the ACL’s relaxation time impact tem- perature, with prolonged relaxation times resulting in lower temperatures. Conversely, prolonged relaxation times of the PCL tend to cause an increase in temperature. Fig. 5(a) and 5(b) depict the opposite effects of Br(A) and Br(P ) on temperature. With the increment in Br(A), the temperature decreases, while with the increment in Br(P ), the temperature rises. This temperature variation is a result of the interplay between thermal diffusion and viscous dissipation. The viscous dissipation rate rises as Br(A) increases, but the thermal diffusion rate falls. Consequently, the temperature rises by the increase in thermal diffusion and falls with a decrease in viscous dissipation. Conversely, as Br(P ) rises, the rate of viscous dissipation rises and the rate of thermal diffusion falls. Because of viscous dissipation, the temperature rises. Conversely, an increment in thermal diffu- sion causes the temperature to decrease. The temperatures of both fluid layers increase when the parameter Π increases, as seen in Fig. 5(c). Fig. 5(d) illustrates the correlation between temperature and ξ, showing that a higher value of ξ correlates to the rise in tem- perature. Thus, there is of directly proportional relation between the temperature and pressure gradient condition. Fig. 6 presents a comparison between the Johnson-Segalman fluid and Newtonian fluid in terms of temperature. This figure suggests that mucus char- acterized by the Johnson-Segalman fluid has a higher temperature than the Newtonian fluid. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 22 of 30 Figure 4: The variation impacts of (a) eA, (b) eP , (c) (We)(A) and (d) (We)(P ) on temperature T . Figure 5: The variation impacts of (a) Br(A), (b) Br(P ), (c) ζ and (d) ξ on temperature T . A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 23 of 30 Figure 6: The comparison between the Johnson-Segalman fluid and Newtonian fluid for temperature T . 4.3. Concentration The variation impacts of eA, eP , (We)(A) and (We)(P ), Br(A), Br(P ), ξ, Π, SH and ST on the ACL and PCL concentrations are plotted in Figs. 7–9. The concentration decreases as the values of eA and eP tend to increment, as delineated in Figs. 7(a) and 7(b). Figs. 7(c) and 7(d) elucidate how (We)(A) and (We)(P ) affect the concentrations of ACL and PCL, respectively. These graphs show that the concentration reduces with increasing (We)(A), while the concentration increases with increasing (We)(P ). The behavior of (We)(A) and (We)(P ) can be associated with the relaxation time of ACL and PCL, respectively. As the ACL relaxation time increases, the concentration reduces. Conversely, PCL concentration increases with increasing relaxing time. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 24 of 30 Figure 7: The variation impacts of (a) eA, (b) eP , (c) (We)(A) and (d) (We)(P ) on concentration C. Figs. 8(a) and 8(b) are plotted to seek the impacts of Br(A) and Br(P ) on concentration. It is delineated that by the increase of Br(A), the concentration increases minutely, whereas by the increase of Br(P ), the concentration decreases. Since Br(k) is directly proportional to viscous dissipation and inversely proportional to thermal diffusion rate. At this, we can add that with the increment in viscous dissipation of ACL, the concentration increases, and with the increment in thermal diffusion rate, the concentration decreases. The recip- rocal impacts of thermal diffusion rate and viscous dissipation rate of PCL are observed. Fig. 8(c) illustrates the effects of ξ on concentration. From this figure, one perceived that the concentration decreases with an increase in ξ—moreover, the more the pressure gradient value at the airway entrance, the lesser the concentration. Fig. 8(d) illustrates that the concentration also rises as the value of Π increases. As Π represents the quotient of D (A) B and D (P ) B . It can be said that Π is directly proportional to D (A) B and inversely proportional to D (P ) B . The decrease in diffusivity of the fluid is the cause of the increase in concentration. Decreased diffusivity generally leads to reduced mixing and slower diffu- sion, resulting in a less complete mixing of the material and increased concentration. The concentration tends to decrease when the values of SH and ST increase, as observed in Figs. 9(a) and 9(b). An increased value of SH elucidates that the mass diffusivity is lower than the viscosity, leading to higher solute concentration and reduced diffusion within the mucus. A significant ST suggests that temperature gradients promote solute trans- port, which affects mucus concentration. These parameters depict the distribution and movement of solutes throughout the mucus layers, ultimately influencing the effectiveness A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 25 of 30 of mucociliary clearance. Fig. 9(c) compares the concentration levels between Johnson- Segalman and Newtonian fluids. It indicates that when the mucus is characterized as a Johnson-Segalman fluid, the concentration exceeds that of the Newtonian fluid. Figure 8: The variation impacts of (a) Br(A), (b) Br(P ), (c) ξ and (d) Π on concentration C. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 26 of 30 Figure 9: The variation impacts of (a) SH and (b) ST on concentration C and (c) comparison between the Johnson-Segalman fluid and Newtonian fluid for concentration C. 5. Concluding Remarks The present prospective analysis has been conducted to theoretically investigate the ther- mal and concentration aspects of human airways using the mathematical model of two- layered mucociliary transport. The model has used the scenario of mucus hypersecretion in COPD. The Johnson-Segalman fluid model has characterized the viscoelastic nature of secreted mucus within the airway channel lumen. The subsequent set of nonlinear PDEs has been first simplified using the theory of lubrication approximation. The ADM was then used to solve subsequent PDEs in series form up to the second-order. Velocity, tem- perature, and concentration have been used as flow variables to characterize mucociliary transport. The variation impact of the parameters of interest like slip parameters (eA and eP ), Weissenberg numbers ((We)(A) and (We)(P )), Brinkmann numbers (Br(A) and Br(P )), pressure gradient at the airway channel entrance ξ, conductivity ratio ζ, mass diffusivity ratio Π, Schmidt number SH , and Soret number ST have been observed. We have concisely outlined the significant results of the present mathematical analysis in the following manner: • When the values of eA, eP , (We)(A), (We)(P ), and ξ increase, the ACL and PCL’s velocities tend to increase. Higher Weissenberg numbers have distinct influences on the velocities of ACL and PCL. Higher values of (We)(P ) accelerate fluid flow in a favorable direction, whereas higher values of (We)(A) reverse the direction of flow. • The ACL and PCL’s temperature elevates by the increase of eP , (We)(P ), Br(P ), ζ and A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 27 of 30 ξ. The temperature of the ACL and PCL layers rises because of viscous dissipation, but thermal diffusion reduces the temperature of these layers when Br(P ) increases. On the other hand, it decreases when eA, (We)(A), and Br(A) increment. An icremet in Br(A) results in an increase in temperature in the ACL and PCL layers because of thermal dif- fusion. However, the temperature within these layers is reduced by viscous dissipation. • The ACL and PCL’s concentration decreases with increasing eA, eP , (We)(A), (We)(P ), ξ, SH and ST while it increases with increasing (We)(P ), Br(A) and Π. The concentration of ACL increases with viscous dissipation, while it decreases with thermal diffusion rate, while PCL’s viscous dissipation and thermal diffusion rate differ. The pressure gradient at the airway channel entrance reduces concentration due to decreased fluid diffusivity. • The comparison between Johnson-Segalman and Newtonian fluids for velocity, temper- ature, and concentration reveals that Johnson-Segalman fluid flows with more velocity, higher temperature, and more concentration. The findings of this analysis suggest that when we characterize secreted mucus by the Johnson-Segalman fluid, mucus clearance in the respiratory system improves, particularly in patients of severe COPD. It’s essential to grasp the rheological properties of mucus and consider environmental influences to develop successful treatments for COPD, particularly in scenarios involving hypersecretion of mucus in the airways. Acknowledgements The authors wish to express their very sincere thanks to the editors and reviewers for their valuable suggestions and comments. Conflict of Interest The authors have no conflicts of interest to declare. Data availability The data used to support the findings of this study are available from the corresponding author upon request. References [1] G Adomian. Analytical solutions for ordinary and partial differential equations. In Differential Equations and Mathematical Physics: Proceedings of an Interna- tional Conference held in Birmingham, Alabama, USA March 3–8, 1986, pages 1–15. Springer, 2006. [2] Noreen Sher Akbar, Salman Akhtar, Ehnber N Maraj, Ali E Anqi, and Raad Z Homod. Heat transfer analysis of mhd viscous fluid in a ciliated tube with entropy generation. Mathematical Methods in the Applied Sciences, 46(10):11495–11508, 2023. [3] Noreen Sher Akbar, Abbasali Abouei Mehrizi, Maimona Rafiq, M Bilal Habib, and Taseer Muhammad. Peristaltic flow analysis of thermal engineering nano model with effective thermal conductivity of different shape nanomaterials assessing variable fluid properties. Alexandria Engineering Journal, 81:395–404, 2023. A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 28 of 30 [4] Amel A Alaidrous and Mohamed R Eid. 3-d electromagnetic radiative non-newtonian nanofluid flow with joule heating and higher-order reactions in porous materials. Scientific Reports, 10(1):14513, 2020. [5] H Ashraf, Tariq Ali, Hamood Ur Rehman, Nehad Ali Shah, Sidra Irshad, and Bander Almutairi. Heat and concentration analysis of two-layered muco-ciliary third grade fluid flow in human airways. Case Studies in Thermal Engineering, 59:104512, 2024. [6] H Ashraf, AM Siddiqui, and MA Rana. Analysis of the peristaltic-ciliary flow of johnson–segalman fluid induced by peristalsis-cilia of the human fallopian tube. Math- ematical biosciences, 300:64–75, 2018. [7] Hameed Ashraf, Imran Siddique, Ayesha Siddiqa, Ferdous MO Tawfiq, Fairouz Tchier, Rana Muhammad Zulqarnain, Hamood Ur Rehman, Shahzad Bhatti, and Abida Rehman. Analysis of two layered peristaltic-ciliary transport of jeffrey fluid and in vitro preimplantation embryo development. Scientific Reports, 14(1):1469, 2024. [8] Rama Bansil and Bradley S Turner. The biology of mucus: Composition, synthesis and organization. Advanced drug delivery reviews, 124:3–15, 2018. [9] John Blake. Mucus flows. Mathematical Biosciences, 17(3-4):301–313, 1973. [10] Ronald G Crystal, Scott H Randell, John F Engelhardt, Judith Voynow, and Mary E Sunday. Airway epithelial cells: current concepts and challenges. Proceedings of the American Thoracic Society, 5(7):772–777, 2008. [11] Henry Danahay and Alan D Jackson. Epithelial mucus-hypersecretion and respiratory disease. Current Drug Targets-Inflammation & Allergy, 4(6):651–664, 2005. [12] MR Eid, SM Abdel-Gaied, and AA Idarous. On effectiveness chemical reaction on viscous flow of a non-darcy nanofluid over a non-linearly stretching sheet in a porous medium. Jokull J, 65(12):76–92, 2015. [13] Christopher M Evans, Kyubo Kim, Michael J Tuvim, and Burton F Dickey. Mucus hypersecretion in asthma: causes and effects. Current opinion in pulmonary medicine, 15(1):4–11, 2009. [14] Yuan-Cheng Fung and Yuan-Cheng Fung. Bioviscoelastic fluids. Biomechanics: Me- chanical Properties of Living Tissues, pages 220–241, 1993. [15] David B Hill, Brian Button, Michael Rubinstein, and Richard C Boucher. Physiology and pathophysiology of human airway mucus. Physiological Reviews, 102(4):1757– 1836, 2022. [16] S Hina, T Hayat, and A Alsaedi. Heat and mass transfer effects on the peristaltic flow of johnson–segalman fluid in a curved channel with compliant walls. International Journal of Heat and Mass Transfer, 55(13-14):3511–3521, 2012. [17] Mohammad Mahdi Hosseini and Hamideh Nasabzadeh. On the convergence of ado- mian decomposition method. Applied mathematics and computation, 182(1):536–543, 2006. [18] MW Johnson Jr and D Segalman. A model for viscoelastic fluid behavior which allows non-affine deformation. Journal of Non-Newtonian fluid mechanics, 2(3):255– 270, 1977. [19] M King. Viscoelastic properties of airway mucus. In Federation proceedings, vol- A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 29 of 30 ume 39, pages 3080–3085, 1980. [20] WL Lee, PG Jayathilake, Zhijun Tan, DV Le, HP Lee, and BC Khoo. Muco-ciliary transport: effect of mucus viscosity, cilia beat frequency and cilia density. Computers & Fluids, 49(1):214–221, 2011. [21] Marcus A Mall. Role of cilia, mucus, and airway surface liquid in mucociliary dys- function: lessons from mouse models. Journal of aerosol medicine and pulmonary drug delivery, 21(1):13–24, 2008. [22] N Manzoor, O Anwar Bég, K Maqbool, and S Shaheen. Mathematical modelling of ciliary propulsion of an electrically-conducting johnson-segalman physiological fluid in a channel with slip. Computer methods in biomechanics and biomedical engineering, 22(7):685–695, 2019. [23] K Maqbool, S Shaheen, E Bobescu, and R Ellahi. Thermal and concentration analysis of phan–thien–tanner fluid flow due to ciliary movement in a peripheral layer. J Cent South Univ, 28(11):3327–3339, 2021. [24] Janna C Nawroth, Anne M Van Der Does, Amy Ryan, and Eva Kanso. Multiscale mechanics of mucociliary clearance in the lung. Philosophical Transactions of the Royal Society B, 375(1792):20190160, 2020. [25] Kh Nowar, EM Abo-Eldahab, and EI Barakat. Peristaltic pumping of johnson- segalman fluid in an asymmetric channel under the effect of hall and ion slip currents. J. Appl. Comput. Math, 1(102):1–6, 2012. [26] Rahul R Rajendran and Arindam Banerjee. Effect of non-newtonian dynamics on the clearance of mucus from bifurcating lung airway models. Journal of Biomechanical Engineering, 143(2):021011, 2021. [27] Samriddha Ray and Jeffrey A Whitsett. Airway mucus and mucociliary system. 1998. [28] Duncan F Rogers. Airway mucus hypersecretion in asthma: an undervalued pathol- ogy? Current opinion in pharmacology, 4(3):241–250, 2004. [29] Duncan F Rogers. Airway mucus hypersecretion: Rationales for pharmacotherapy. Journal of Organ Dysfunction, 2(3):183–191, 2006. [30] SM Ross and S Corrsin. Results of an analytical model of mucociliary pumping. Journal of Applied Physiology, 37(3):333–340, 1974. [31] Bruce K Rubin. Secretion properties, clearance, and therapy in airway disease. Trans- lational respiratory medicine, 2:1–7, 2014. [32] Mohammed R Salman. The ciliary propulsion of an electrically conducting johnson- segalman physiological fluid through a porous medium in an inclined symmetric a channel with slip. In Journal of Physics: Conference Series, volume 1818, page 012192. IOP Publishing, 2021. [33] Mohammad Hadi Sedaghat, Uduak Z George, and Omid Abouali. A nonlinear vis- coelastic model of mucociliary clearance. Rheologica Acta, 60(6):371–384, 2021. [34] Binay Kumar Shah, Bivek Singh, Yukun Wang, Shuanshuan Xie, and Changhui Wang. Mucus hypersecretion in chronic obstructive pulmonary disease and its treat- ment. Mediators of Inflammation, 2023(1):8840594, 2023. [35] Noreen Sher Akbar and Salman Akhtar. Metachronal wave form analysis on cilia- driven flow of non-newtonain phan–thien–tanner fluid model: A physiological math- A. Alaidrous, H. Ashraf / Eur. J. Pure Appl. Math, 18 (1) (2025), 5581 30 of 30 ematical model. Proceedings of the Institution of Mechanical Engineers, Part E: Journal of Process Mechanical Engineering, 237(6):2567–2573, 2023. [36] Abdul-Majid Wazwaz. Partial differential equations and solitary waves theory. Springer Science & Business Media, 2010. [37] A Yaghi and MB Dolovich. Airway epithelial cell cilia and obstructive lung disease. cells, 5 (4), 40, 2016. Appendix Λ = 1−M , G (k) 1 = ( (We)(k) )2 ( 1− e2k ) , G (k) 2 = Br(k) M , G (k) 3 = − ( (We)(k) )2 ( 1− e2k ) Λ, G (k) 4 = ( (We)(k) )4 ( 1− e2k ) , G (k) 5 = SHST D (k) B , J (k) 1 = G (k) 2 G (k) 5 , Γ1 = G (A) 1 Λ, Γ2 = G (P ) 1 Λ, Γ3 = (G (A) 1 )2Λ, Γ4 = (G (P ) 1 )2Λ, Γ5 = ξ 6 [ 3 (δ − δϵ)2 (δ + δϵ+ 1) + (1 + ϵ) (3δ + 3δϵ+ 2) ] , Γ6 = 1 + 2πϵαβ + 1 2 ( 2ϵ+ πϵ2αβ ) , Γ7 = (δ + δϵ)5 + ( 4 (1 + ϵ)− 5δ (1 + ϵ) (1 + ϵ)4 ) , Γ8 = δ (1 + ϵ) ( (δ + δϵ)4 − (1 + ϵ) ) , Γ9 = (δ + δϵ)5, Γ10 = (δ + δϵ)7,