Analytical Investigation of Time-Dependent Two-Dimensional Non-Newtonian Boundary Layer Equations

Abstract

In this study, five different time-dependent incompressible non-Newtonian boundary layer models in two dimensions are investigated with the self-similar Ansatz, including external magnetic field effects. The power-law, the Casson fluid, the Oldroyd-B model, the Walter fluid B model, and the Williamson fluid are analyzed. For the first two models, analytical results are given for the velocity and pressure distributions, which can be expressed by different types of hypergeometric functions. Depending on the parameters involved in the analytical solutions of the nonlinear ordinary differential equation obtained by the similarity transformation, a vast range of solution types is presented. It turned out that the last three models lack self-similar symmetry; therefore, no analytic solutions can be derived.

1. Introduction

The scientific field of classical fluid mechanics is vast and cannot be summarized in a few finite arbitrary books. If we restrict ourselves to non-Newtonian fluids or boundary layer flows, the relevant literature is still remarkably large. The fundamental physics of such fluid flows can be found in introductory textbooks such as [1,2,3,4,5]. In our previous publication on time-dependent self-similar solutions of compressible and incompressible heated boundary layer equations [6], we collected the most relevant literature in the field, which we omit here. Some analytic results are also available for stationary non-Newtonian boundary layers, for example, [7,8,9].

Today, there is a growing need to study the motion of non-Newtonian fluids, which are often found in industrial applications and in modeling many manufacturing processes (heat exchangers, pharmaceuticals, food, and paper). All fluids for which Newton’s law cannot describe the shear stress rate relation are called non-Newtonian fluids. They can be characterized by formulas of a highly diverse nature. Without completeness, we will investigate five different non-Newtonian fluid models in our study; these are the following: The first is the non-Newtonian power-law model, which describes the viscosity of time-independent flow behavior for both shear thinning (or pseudoplastic) and shear thickening (or dilatant) cases (see [1,4]). This model is used in vast engineering applications, for instance, lubricants, polymers, and slurries [5]; this model is well suited for numerical and analytical studies of fluid flows [8].

The second can be directly derived from the usual Newtonian fluids. This is the widely used viscosity model for viscoplastic non-Newtonian behavior and is named the Casson fluid model. Fluids that behave like a solid when the shear stress is less than the yield stress applied to a liquid start to deform when the shear stress is greater than the yield stress [10]. Numerical studies regarding the Casson fluid one may find in [11]. The system of unsteady MHD fluid that flows through a porous medium is discussed in [12].

Our third model is the non-Newtonian Oldroyd-B fluid, which is used in the literature for the rheological characterization of a type of viscoelastic fluid. Oldroyd applied the principle that stresses in a continuous medium can only arise from deformations and cannot change if the material is only rotated [13].

The fourth is the Walters-B model, which was proposed by Walters for viscoelastic fluids [14]. It is a generalization of the Oldroyd-B model and describes the behavior of fluids mainly used in food and food processing technologies.

The final model is the Williamson fluid model, which can also describe the viscoelastic shear-thinning properties of non-Newtonian fluids [15]. In Williamson’s fluid model, the effective viscosity should decrease infinitely as the shear rate increases, which is effectively infinite viscosity at rest (no fluid motion) and zero viscosity as the shear rate approaches infinity.

The actual existing literature on non-Newtonian boundary layers is enormous. Therefore, we cannot give a general overview; we can mention some relevant recent publications in the right place in this study.

To analyze the variation in flow characteristics of a fluid in a magnetic field, an additional magnetic term is added to the momentum equation of the boundary layer equation system so that the system can describe the magnetic effects.

The aim of this paper is twofold: firstly, to derive analytic results for the non-Newtonian boundary layer equations and secondly, to investigate how an additional stationary magnetic field changes these effects.

We examine the five types of non-Newtonian behavior given above. A similarity transformation transforms the governing equations into a system of ordinary differential equations. We analyze the cases in which the fluid flow rate functions can be analytically given by the method used. In these two cases, we give the velocity distributions in two space dimensions and point out the effect of each parameter.

2. Theory and Results

In the following, we present our various non-Newtonian boundary layer models with an additional magnetic field after some introductory remarks. This gives readers a hint of how the models were derived and introduced. As the main point, we investigate boundary layers. Therefore, our starting point is the pioneering work of Prandtl [16], who used scaling arguments and derived that specific terms of the Navier–Stokes equations are negligible in boundary layer flows. Some years later, Blasius [17] gave the solutions of the steady-state, incompressible, two-dimensional, laminar boundary layer equation forms on a semi-infinite flat surface parallel to a constant unidirectional flow. Here, we study the uniform two-dimensional flow of an incompressible, viscous, electrically conducting fluid in a tension plate placed in a uniform velocity ambient fluid. Assuming a magnetic field $B_0$ perpendicular to the plate for small magnetic Reynolds number, the induced magnetic field justified for MHD flow is neglected, i.e., $R{e}_m={\mu }_0\sigma VL\ll 1$, where ${\mu }_0$ is the magnetic permeability, $\sigma$ is the electrical conductivity of the fluid, V is the characteristic velocity, L is the characteristic length scale of the fluid, and M is the characteristic length scale. It is also assumed that the external electric field is zero and that the electric field polarization is negligible. The next step is introducing the non-Newtonian viscous terms in the impulse equation and modifying the term in the remaining Navier–Stokes equation after simplifying the boundary layer. This is performed by hand, just writing down the power-law, Oldroyd-B types, the Walters-T type, or even the Williamson type non-Newtonian viscosity. Finally, an additional term is needed if an external magnetic field is present. Hartmann first derived this in 1937 [18] for the theory of laminar flow on an electrically conductive liquid in a homogeneous magnetic field. In recent years, the same term was used by Waqas [19] and Lone [20].

Secondly, we have a few statements about the boundary conditions of our future solutions. In this study, we generally obtain our results as solutions of second-order ordinary differential equations, meaning we obtained three free real integral constants: $c_1$, $C_1$, and $C_2$. The solutions are a linear combination of the regular and the irregular solutions. When the boundary conditions are given at the boundary of the sheet and far from the boundary at the boundary layer edge, we can obtain the similarity solutions for the velocities and pressure as a function of the space and time in analytic forms.

Unfortunately, with this self-similar Ansatz, not all the mathematically and physically possible boundary conditions can be described. The good news is that physically relevant solutions with proper asymptotic temporal and spatial decay can be defined and studied. This is the shortcomings of the method.

2.1. Power-Law Viscosity

Under our assumptions, the MHD equations for the steady two-dimensional flow in the boundary layer for a fluid with power-law viscosity read as follows:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}& =& 0,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \rho \left({\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}\right)& =& -{\displaystyle \frac{\partial p}{\partial x}}+m{\displaystyle \frac{\partial }{\partial y}}\left({\left|{\displaystyle \frac{\partial u}{\partial y}}\right|}^{n-1}{\displaystyle \frac{\partial u}{\partial y}}\right)+\sigma B_0^{2}u,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial p}{\partial y}}& =& 0,\hfill \end{array}$$

(3)

where the dynamical variables are the velocity components $u(x,y,t),v(x,y,t)$, the fluid pressure $p(x,y,t)$, $B_0=const$ is the magnetic induction, and $\sigma$ is the electrical conductivity. Additional physical parameters are $\rho ,m,n$, the fluid density, consistency parameter, and the power-law index, respectively. Consider Figure 1 to fix our system’s geometrical relations.

The boundary conditions for the fluid flow with constant velocity $U_e$ and using a continuously moving permeable surface are as follows:

$$u(x,y=0,t)=U_w,v(x,y=0,t)=V_w,u(x,y\to \infty ,t)=U_e.$$

(4)

Here, we take $U_e$, $U_w$, and $V_w$ as constants.

The temporal decay of the strong solutions of the power-law type of non-Newtonian fluids for variable power-law index was mathematically proven by Ko [21]. Flow reversal effects in an expanding channel for this type of viscous fluid were analyzed by [22]. A non-iterative transformation method for an extended Blasius problem describing a 2D laminar boundary layer with power-law viscosity for non-Newtonian fluids was developed by Fazio [23]. Pater et al. [24] solved the stationary power-law equations by the one-parameter deductive group theory technique. We mention at this point that the stationary boundary layer equations can be reduced to the third-order non-linear differential equation, which is called the Blasius equation and has extensive literature; see [4,25].

To obtain the solutions to the system (1)–(3), we apply the following self-similar Ansatz for $u,v$ and p [25]:

$$\begin{array}{c}\hfill u(x,y,t)=t^{-\alpha }f\left(\eta \right),v(x,y,t)=t^{-\delta }g\left(\eta \right),p(x,y,t)=t^{-\gamma }h\left(\eta \right),\end{array}$$

(5)

with the argument $\eta ={\displaystyle \frac{x+y}{t^{\beta }}}$ of the shape functions. All the exponents $\alpha ,\beta ,\gamma ,\delta$ are real numbers. Solutions with integer exponents are called self-similar solutions of the first kind; non-integer exponents generate self-similar solutions of the second kind [25].

It is necessary to remark that if we add the regular heat conduction mechanism to (1)–(3), the boundary layer equations exclude the self-similarity; therefore, using the Ansatz of (5), one obtains a contradiction among self-similar exponents.

The shape functions $f,g$, and h could be any arbitrary continuous functions with existing first and second continuous derivatives and will be derived later on. The physical and geometrical interpretation of the Ansatz was exhaustively analyzed in former publications [26]; therefore, we skip it here.

The main points are that $\alpha ,\delta ,\gamma$ are responsible for the decay rate, and $\beta$ is for the spreading rate of the corresponding dynamical variable for positive exponents. Negative exponents mean physically irrelevant cases, blowing up and contracting solutions. The numerical values of the exponents are considered as follows:

$$\alpha =\delta =n/2,\beta =1-\alpha =1-n/2,\gamma =n.$$

(6)

Exponents with numerical values of one-half mean the regular Fourier heat conduction (or Fick’s diffusion) process. One-half values for the exponent of the velocity components and unit value exponent for the pressure decay are usual for the incompressible Navier–Stokes equation [26]. Note that the value of $\gamma$ is responsible for the decay of the pressure field.

Applying (6), the derived ordinary differential equation (ODE) system reads as

$$f^{\prime }+g^{\prime }=0,$$

(7)

$$\rho \left[-{\displaystyle \frac{n}{2}}f-\left(1-{\displaystyle \frac{n}{2}}\right)\eta f^{\prime }+f{f}^{\prime }+g{f}^{\prime }\right]=-h^{\prime }+mn{f}^{″}{\left|f^{\prime }\right|}^{n-1}+\sigma B_0^{2}f,$$

(8)

$$h^{\prime }=0,$$

(9)

where the prime means derivation concerning variable $\eta$. Equations (7) and (9) can be integrated to obtain $f+g=c_1$ and $h=c_2$. However, the dynamic variable under consideration is the velocity component u, which is f. From (8), the derived ODE is the following:

$${\displaystyle \frac{mn}{\rho }}{f}^{″}{\left|f^{\prime }\right|}^{n-1}+\left[\left(1-{\displaystyle \frac{n}{2}}\right)\eta -c_1\right]f^{\prime }+Kf=0,$$

(10)

where $K=(\rho n/2+\sigma B_0^2).$ The constant $c_1$ can be determined from possible boundary conditions applied to the fluid flow problem.

Here, we give analytic and numeric solutions to the flow Equations (7)–(9):

(i)

First, we consider the case $B_0=0$. Equation (10) has analytic solutions only for the cases of $n=1,0$, and $-1$.

For the Newtonian fluid case, when $n=1$, and other parameters ($c_1,m,\rho$) are general, the solution reads as

$$f\left(\eta \right)=C_{1}M\left[{\displaystyle \frac{\rho }{2}},{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {(\eta -2c_1)}^2}{4m}}\right]+C_{2}U\left[{\displaystyle \frac{\rho }{2}},{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {(\eta -2c_1)}^2}{4m}}\right],$$

(11)

where $M(,,)$ and $U(,,)$ are the Kummer’s M and Kummer’s U functions [27], and $C_1$ and $C_2$ stand for the usual integration constants. The parameter $c_1$ shifts the maxima of the shape function, and m, the consistency parameter (which is positive), just changes the full width at half the function’s maximum. The larger the numerical value of m, the broader the shape function. The third parameter, the density of the fluid (which must also be positive), is the most relevant parameter of the flow, which, of course, meets our physical intuition. Figure 2 shows four different velocity shape functions for different fluid densities. It is essential to mention that for density $0<\rho$, the shape function has a global maxima in the origin and a decay to zero (larger densities mean quicker decay). However, for $1<\rho$, the shape function also shows additional oscillations. The final velocity field has the asymptotic form of $u(x,y=0,t)=t^{-{\displaystyle \frac{1}{2}}}f\left({\displaystyle \frac{x}{t^{1/2}}}\right)$. For completeness, we give the entire formula for u as well:

$$u(x,y,t)={\displaystyle \frac{1}{t^{1/2}}}\left(C_{1}M\left\{{\displaystyle \frac{\rho }{2}},{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {\left[\frac{(x+y)}{t^{1/2}}-2c_1\right]}^2}{4m}}\right\}+C_{2}U\left\{{\displaystyle \frac{\rho }{2}},{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {\left[\frac{(x+y)}{t^{1/2}}-2c_1\right]}^2}{4m}}\right\}\right),$$

(12)

where $c_1$, $C_1$, and $C_2$ are subjected to the boundary conditions. We remark that this case was investigated in [6]. Figure 3 shows the velocity field of $v(x,y=0,t)$ for a parameter set. The slight oscillation and the quick decay are clear to see. It is interesting to mention here that during our investigations using the self-similar Ansatz, we regularly obtain solutions that contain Kummer’s M and Kummer’s U functions for hydrodynamic processes like the Bénard-Rayleigh convection [26] or even for regular diffusion equations [28].

For $n=1$, the velocity u for long times has a power-law decay $u\sim {\displaystyle \frac{1}{t^{\alpha }}}\sim {\displaystyle \frac{1}{t^{1/2}}}$.

In the case $n=K=0$, the solution is a trivial constant $f=c_2$.

There are analytic solutions available for $n=-1$; unfortunately, not for the general case where all the other three parameters $(c_1,\rho ,m)$ are entirely free. The constant $c_1$ should equal zero, and the ratio of $\rho /m$ should have some special values. The existing possibilities are not so many.

In the following, we present some solutions that came from our numerical experimental experiences:

For $\rho =1$ and arbitrary integer or rational m, we obtain the following implicit relation for the solutions:

$$ln\left(\eta \right)-ln\left[f\left(\eta \right)\right]-ln\left\{C_{1}U\left[-1,{\displaystyle \frac{3}{2}},{\displaystyle \frac{1}{4m}}f{\left(\eta \right)}^2\right]+M\left[-1,{\displaystyle \frac{3}{2}},{\displaystyle \frac{1}{4m}}f{\left(\eta \right)}^2\right]\right\}-C_2=0.$$

(13)

If the first parameters of both Kummer’s functions are non-positive integers, the series becomes finite:

$$\begin{array}{c}\hfill M\left[-1,{\displaystyle \frac{3}{2}},{\displaystyle \frac{1}{4m}}f{\left(\eta \right)}^2\right]=1-{\displaystyle \frac{2}{3m}}f{\left(\eta \right)}^2,\\ \hfill U\left[-1,{\displaystyle \frac{3}{2}},{\displaystyle \frac{1}{4m}}f{\left(\eta \right)}^2\right]=-{\displaystyle \frac{3}{2}}+{\displaystyle \frac{1}{m}}f{\left(\eta \right)}^2,\end{array}$$

(14)

considering that $C_2>0$ and ${\tilde{C}}_2=ln\left(C_2\right)$ after some trivial algebraic steps, we obtain

$$\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)f{\left(\eta \right)}^3+\left({\tilde{C}}_2-{\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}\right)f\left(\eta \right)-\eta =0.$$

(15)

Three different velocity shape functions are presented in Figure 4 for various parameter sets to show the global and general properties of the solution Equation (15).

One can see in Equation (15) that if $\eta =0$, then $f(\eta =0)=0$ is a solution of the equation. Further possibilities give the following rearrangement:

$$\left[\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)f{\left(\eta \right)}^2+\left({\tilde{C}}_2-{\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}\right)\right]f\left(\eta \right)=0.$$

(16)

This yields

$$f(\eta =0)=\pm \sqrt{\left({\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}-{\tilde{C}}_2\right)/\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)}.$$

(17)

It shows that for $\eta =0$, there are one or three corresponding values to function $f\left(\eta \right)$. In the following, we try to find the constraints on the parameters that separate the two cases. If we consider the inverse function $f\left(\eta \right)=y$, we have

$$\eta =\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)y^3+\left({\tilde{C}}_2-{\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}\right)y.$$

(18)

The derivative of this function is

$${\displaystyle \frac{d\eta }{dy}}=3\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)y^2+\left({\tilde{C}}_2-{\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}\right).$$

(19)

If this expression has always the same sign, the function is monotonous, and then there is just one single root $\eta (y=0)=0$ (i.e., $f(\eta =0)=0$). If this derivative may change the sign, then multiple roots are possible. The inverse function is not monotonous and may have three roots if there are real roots of the equation of the derivative $d\eta /dy=0$. This means the following condition:

$${\displaystyle \frac{\left(\frac{3{\tilde{C}}_2{C}_1}{2}-{\tilde{C}}_2\right)}{\left(\frac{{\tilde{C}}_2{C}_1}{m}-\frac{2{\tilde{C}}_2}{m}\right)}}>0.$$

(20)

Figure 4 shows two implicit solutions of Equation (15) for three different parameter sets. Note that Equation (15) is a third-order equation. We may obtain multi-valued solutions. If the parameter m is much smaller than the two integration constants, then single-valued solutions emerge as well. Whether the multi-valued solutions have any physical significance is not yet clear, but this property may indicate the existence of finite oscillations or eddies. Figure 5 shows the projected velocity fields (this means $y=0$) for single- and double-valued cases. The decrease of the velocity field with time is visible. For the sake of completeness, we give the final form of the velocity field, which reads as

$$\left({\displaystyle \frac{{\tilde{C}}_2{C}_1}{m}}-{\displaystyle \frac{2{\tilde{C}}_2}{m}}\right)t^{3/2}f{\left({\displaystyle \frac{x+y}{t^{3/2}}}\right)}^3+\left({\tilde{C}}_2-{\displaystyle \frac{3{\tilde{C}}_2{C}_1}{2}}\right)t^{1/2}f\left({\displaystyle \frac{x+y}{t^{3/2}}}\right)-{\displaystyle \frac{x+y}{t^{3/2}}}=0.$$

(21)

There are analytic solutions available for $\rho =3$ and for arbitrary real m and $c_1$ in the implicit form of

$${\displaystyle \frac{\eta }{f\left(\eta \right)}}+{\displaystyle \frac{C_1{e}^{\frac{9f{\left(\eta \right)}^2}{4m}}}{3f\left(\eta \right)}}-{\displaystyle \frac{C_1·\mathrm{erf}\left[\frac{3}{2}\sqrt{-\frac{1}{m}}f\left(\eta \right)\right]}{2m\sqrt{-\frac{1}{m\pi }}}}-{\displaystyle \frac{2c_1}{3f\left(\eta \right)}}-C_2=0,$$

(22)

where erf() stands for the usual error function. For more information please consult the basic handbook of [27].

Equation (22) clearly shows that for $C_1=0$, we obtain the trivial linear function as a solution. Figure 6 presents the solutions of Equation (22) for three different parameter sets.

The implicitplot command of Maple 12 was used to evaluate the presented functions in Figure 6. As we can see, multivalued functions are presented. This is because the implicit function is second-order in f(eta).

The result functions increase strongly (or decrease in the case of a negative branch) very close to the origin and for larger arguments on a nearly horizontal plateau. The general structure of the solution is relatively stable; even a tenfold change in the parameters does not change the overall feature of the derived curve. To have a better overview of the properties of the solution, Figure 7 presents the usual velocity projection of $u(x,y=0,t)$ for a given parameter set from an unusual point of view. The physically relevant rapid decay over time is again clearly visible.

(ii)

In the second part of our analysis, we investigate the solutions when the magnetic induction is not equal to zero ($B_0\ne 0$). Even now four different cases exist: $n=2,1,0$, and $-1$. We have to examine them one by one.

First, we consider the dilatant case $n=2$. Equation (10) is reduced to

$${\displaystyle \frac{2m}{\rho }}{f}^{″}{f}^{\prime }+Kf=0,$$

(23)

multiplying by $f^{\prime }$, we obtain a total derivative, which can be integrated, giving us the ODE of

$${\displaystyle \frac{2m}{3\rho }}{f}^{\prime 3}+{\displaystyle \frac{1}{2}}(\rho +\sigma B_0^2)f^2+c_0=0.$$

(24)

First considering the $c_0=0$ simpler case, the ODE can be directly integrated and gives multiple solutions: there is a trivial one of $f=0$, there are two complex conjugated solutions that we skip, and a relevant real one. After some algebraic manipulation, it reads as

$$f=-{\displaystyle \frac{6\rho (\rho +\sigma B_0^2){(\eta -c_1)}^3}{216m}},$$

(25)

which is a third-order parabola in the variable of $\eta$. The final velocity distribution is

$$u(x,y,t)=-{\displaystyle \frac{6\rho (\rho +\sigma B_0^2){(x+y-c_1)}^3}{216mt}},$$

(26)

which is divergent at $t=0$ for large spatial coordinates and has quick $t^{-1}$ decay in time in general. Figure 8 shows the graph of Equation (26) for a given parameter set.

For $c_0\ne 0$, there is a formal implicit solution containing an integral:

$$\eta -{\int }^{f\left(\eta \right)}{\displaystyle \frac{2m}{{[-6\rho m^2\left(a^2\rho +a^2\sigma B_0^2+2c_0\right)]}^{1/3}}}da-c_1=0.$$

(27)

Unfortunately, it cannot be evaluated for general arbitrary parameters $(\rho ,m,\sigma ,B_0^2,c_0)$.

The second case is $n=1,$ the Newtonian fluid case:

$$\begin{array}{c}\hfill f\left(\eta \right)=C_1·M\left[{\displaystyle \frac{\rho }{2}}+\sigma B_0^2,{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {(\eta -2c_1)}^2}{4m}}\right]+C_2·U\left[{\displaystyle \frac{\rho }{2}}+\sigma B_0^2,{\displaystyle \frac{1}{2}},-{\displaystyle \frac{\rho {(\eta -2c_1)}^2}{4m}}\right].\end{array}$$

(28)

Note the interesting feature that both the density and the magnetic induction give contributions to the first parameter of Kummer’s functions. The density $\rho$ and the square of the magnetic induction $B_0^2$ are always positive; the electric conductivity $\sigma$ is also positive for regular materials. (With the need for completeness, we have to note that so-called metamaterials can have negative electrical conductivity, but we will skip that case now. More on materials can be found in [29]. For such media, the first parameter of Kummer’s function would be negative, which would mean divergent velocity fields, which contradicts energy conservation.) The external magnetic field enhances the first positive parameter of Kummer’s functions for regular materials, which means more oscillations and quicker decay. The shape function, or the final velocity distribution $u(x,y,t)$, is very similar to what was presented in Figure 2 and Figure 3. If the first parameter of Kummer’s M function $(\rho /2+\sigma B_0^2)$ is larger than one, the larger the magnetic induction, the larger the number of oscillations. This is true for the density and for the electric conductivity as well. Such shape functions are presented in Figure 9.

The third case is $n=0$, which is the simplest one. The original ODE of Equation (10) is reduced to

$$\eta f^{\prime }+\sigma B_0^{2}f=0,$$

(29)

with the trivial solution of

$$f\left(\eta \right)=C_1{\eta }^{-\sigma B_0^2}.$$

(30)

The velocity distribution has the form of

$$u(x,y,t)=C_1{\left({\displaystyle \frac{x+y}{t}}\right)}^{-\sigma B_0^2}.$$

(31)

The pressure field is a pure constant as well. These are simple power-laws. Therefore, we skip to present additional figures.

The last case, $n=-1$, is again the most complicated one. There is no analytic solution available when all parameters are arbitrary. There is, however, a special case when the last term in Equation (10) is zero, which defines the constraint of $\rho =2\sigma B_0^2$. If we additionally fix $c_1=0$, we obtain the following solution:

$$f\left(\eta \right)=\pm {\displaystyle \frac{1}{3}}\sqrt{{\displaystyle \frac{6m}{3\rho }}}\arctan \left({\displaystyle \frac{\sqrt{6m\rho }\eta }{\sqrt{-6m\rho {\eta }^2+2m{C}_1}}}\right)+C_2,$$

(32)

where arctan() is the usual inverse trigonometric arcus tangent function. It can be easily shown that the solution has a compact support, and the function is only defined in the region of $-\sqrt{{\displaystyle \frac{C_1}{3\rho }}}\le \eta \le \sqrt{{\displaystyle \frac{C_1}{3\rho }}}$. Note that the smaller the fluid density, the larger the acceptable velocity range if $c_1$ remains the same. The larger the integral constant $c_1$ parameter, the larger the available velocity range. The second integral constant $C_2$ just shifts the solution parallel to the y-axis. Figure 10 shows the solution for three different parameter sets. Figure 11 presents the final typical velocity distribution $v(x,y=0,t)\sim t^{{\displaystyle \frac{1}{2}}}f\left({\displaystyle \frac{x}{t^{3/2}}}\right)$. The quick temporal decay is clear to see. Note that the prefactor $t^{{\displaystyle \frac{1}{2}}}$ makes the time backpropagation impossible for negative values.

For the sake of completeness, we provide solutions for the pressure as well. The ODE of the shape function is trivial with the solution of

$$h^{\prime }=0,h=c_2.$$

(33)

Therefore, the final pressure distribution reads as

$$p(x,y,t)=t^{-\gamma }h(x,y,t)={\displaystyle \frac{c_2}{t^n}},$$

(34)

which means that the pressure is constant in the entire space at a given time, and the time decay can be different from the velocity field.

2.2. Casson Fluid

We investigated additional fluid models to perform an even more comprehensive analysis and study how the transients happen in non-Newtonian boundary layers.

The first in our line is the simplest one, the so-called Casson fluid model; now the boundary layer equation has the form of

$$\rho \left({\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}\right)=-{\displaystyle \frac{\partial p}{\partial x}}+\left(1+{\displaystyle \frac{1}{\lambda }}\right)\mu {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\sigma B_0^{2}u,$$

(35)

where $\lambda$ is the Casson parameter. Note that for $\lambda \to \infty$ the term goes over the regular Newtonian viscous term. (All five physical parameters $\rho ,\lambda ,\mu ,B_0$, and $\sigma$ should have positive real values, and the $\lambda \ne 0$ is also evident.) Sochi analyzed the variational approach for the flow of Casson fluids in pipes [30]. The blood flow through a stenotic tube was modeled by the Casson fluid by Tandon et al. [31].

We still consider the self-similar Ansatz of (5) for the three dynamical variables. After the usual algebraic manipulations, we arrive at the non-linear second-order ODE for the shape function of the horizontal velocity component $f\left(\eta \right)$ in the next form:

$$\rho \left(-{\displaystyle \frac{f}{2}}-{\displaystyle \frac{\eta f^{\prime }}{2}}\right)=\left(1+{\displaystyle \frac{1}{\lambda }}\right)\mu f^{″}+\sigma B_0^{2}f,$$

(36)

In this model, all four self-similar exponents have fixed values: $\alpha =\beta =\delta =1/2,\gamma =1$, which is usual for the Newtonian viscous fluid equations [26]. The derived solutions are mainly different from the first model; therefore, we have to give a detailed analysis in the following. With our usual mathematical program package, Maple 12, we can easily derive the solution in closed form:

$$\begin{array}{c}\hfill f\left(\eta \right)=e^{-{\displaystyle \frac{\rho {\eta }^2}{4\mu \left(\frac{\lambda +1}{\lambda }\right)}}}·\eta ·\left\{C_{1}M\left[-{\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }},{\displaystyle \frac{3}{2}},{\displaystyle \frac{\rho {\eta }^2}{4\mu \left(\frac{\lambda +1}{\lambda }\right)}}\right]\right.+\\ \hfill \left.C_{2}U\left[-{\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }},{\displaystyle \frac{3}{2}},{\displaystyle \frac{\rho {\eta }^2}{4\mu \left(\frac{\lambda +1}{\lambda }\right)}}\right]\right\}.\end{array}$$

(37)

It is important to note, at this point, that the derived formula shows some similarities to the results of the regular diffusion equations [28], which were exhaustively discussed in our former studies. The situation is, however, a bit different here. Thanks to the positivity of all parameters, the exponent sign is always negative (a real Gaussian function), which dictates a rapid decay for large arguments $\eta$ and for any additional parameter sets. The crucial parameter that qualitatively classifies the solutions is the numerical value of the first parameter of Kummer’s functions. We can distinguish three cases:

• $-\left({\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }}\right)<0$, which is equivalent to $2\sigma B_0^2>\rho$ (with the stipulation of $\rho \ne 0$). The solutions have oscillations. If this parameter is a negative integer, the infinite series of Kummer’s function will break down into a finite-order polynomial. The smaller the parameter, the larger the number of oscillations.
• $-\left({\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }}\right)=0$, which is equivalent to $2\sigma B_0^2=\rho$ (with the stipulation of $\rho \ne 0$). Now both Kummer’s M and Kummer’s U functions are unity, and the solution is reduced to the Gaussian function times $\eta$ function, which has odd symmetry. This is the limiting solution between the oscillating and the non-oscillating solutions.
• $-\left({\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }}\right)>0$, which is equivalent to $2\sigma B_0^2<\rho$ (with the stipulation of $\rho \ne 0$). The larger the parameter, the smaller the peak value of the solution and the quicker the decay.

Figure 12 presents five shape functions of Equation (37) for five different fluid densities; for a clear comparison, all the other parameters remain the same. The larger the density $\rho$, the quicker the decay and the smaller the global maximum of the shape functions.

In general, the $\rho /[4\mu (\lambda +1)/\lambda ]$ factor is responsible for the extent (or the full width at half maximum (FWHM) value of the solution). The larger this parameter is, the quicker the decay will occur. For changing the Casson parameter in the physically relevant positive range $\lambda >0$, we found no remarkable change in the solutions; this is because the ${\displaystyle \frac{\lambda +1}{\lambda }}$ (is almost unit) factor appears in the Gaussian together with two other parameters;

The role of magnetic induction is also important. Figure 13 shows the effect of the magnetic induction $B_0^2$, where all the other parameters are unchanged. If the first parameter of Kummer’s M function is smaller than unity, the solutions have a global maximum and a quick decay; if this parameter is larger than the unity, then the larger the magnetic induction the larger the oscillations of the solutions.

Most of these solutions have odd symmetry. It can be easily seen because the Gaussian is an even function, $\eta$ is an odd function, and the series of Kummer’s M and U functions, which have a quadratic argument, are again even infinite or finite series, which, together, give odd symmetry. However, using the series expansion, it can be shown that if the first parameter of Kummer’s U (only for U) function $-\left({\displaystyle \frac{2\sigma B_0^2-\rho }{2\rho }}\right)$ is a negative half-integer, then the solutions obtain even symmetry. The properties of such kinds of solutions are exhaustively discussed in [28].

To complete our analysis, we will also provide solutions for the pressure. The ODE of the shape function is trivial with the solution of

$$h^{\prime }=0,h=c_2.$$

(38)

Therefore, the final pressure distribution is

$$p(x,y,t)=t^{-1}h(x,y,t)={\displaystyle \frac{c_2}{t^1}},$$

(39)

which means that the pressure is constant in the entire space at a given time, and the temporal decay follows the simple inverse law.

2.3. Oldroyd-B Model

The next fluid flow model is the Oldroyd-B model [13], where the momentum equation is the following:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{{\lambda }_1}{\rho }}\left(u^2{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+2uv{\displaystyle \frac{{\partial }^{2}u}{\partial x\partial y}}+v^2{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{\nu }{\rho }}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\\ \hfill +{\displaystyle \frac{\nu {\lambda }_2}{\rho }}\left(u{\displaystyle \frac{{\partial }^{3}u}{\partial x\partial y^2}}+v{\displaystyle \frac{{\partial }^{3}u}{\partial y^3}}-{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}-{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{{\partial }^{2}v}{\partial y^2}}\right)+\sigma B_0^{2}u.\end{array}$$

(40)

Even without the external term with the magnetic induction, the self-similar Ansatz of (5) leads to a contradiction among the four self-similar exponents, indicating that the system has no self-similar symmetry with power-law time decay. The main cause of the lack of symmetry is the appearance of the third spatial derivatives together with the second ones. When terms with second spatial derivatives and terms with third spatial derivatives appear together in one equation, then the exponents of the explicit time dependence cause contradictions that cannot be solved (e.g., $t^{-\alpha -2\beta }$ should be equal to $t^{-\alpha -3\beta }$ ), and on the other side, an additional equation—usually the continuity equation—dictates the $\beta =1/2$ condition. In such cases, the dynamical system has no self-symmetry. So, one cannot see if a PDE system has self-symmetry in such cases.

The Oldroyd-B model was exhaustively investigated in the last decades from many points of view; optimal time decay rates for the higher order spatial derivatives of solutions were analyzed by Wang [32]. The global well-posedness of the model was investigated by Elgindi and Liu [33]. The global existence results of some Oldroyd-B models were proven by [34]. Finally, we mention the review paper of [35], which summarizes most of the known mathematical results.

2.4. Walters’ Liquid B Model

The third possible non-Newtonian candidate is the so-called Walters’ Liquid B model [14], where the momentum equation has the form of

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{\mu }{\rho }}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+\\ \hfill {\displaystyle \frac{k_0}{\rho }}\left(u{\displaystyle \frac{{\partial }^{3}u}{\partial x\partial y^2}}+v{\displaystyle \frac{{\partial }^{3}u}{\partial y^3}}-{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{{\partial }^{2}u}{\partial x\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)+\sigma B_0^{2}u.\end{array}$$

(41)

Practically, we observed the same property as in the previous case: the existence of the terms with the third spatial partial derivatives destroyed the self-similar symmetry, causing controversy beyond the self-similar exponents. Therefore, no solutions can be derived with this Ansatz.

Flow and heat transfer of Walter’s liquid model with stretching walls were applied to hemodynamics by Misra, Shit, and Rath [36]. Khrisna [37] studied the Hall and ion slip effects on MHD laminar flow in a Walter’s-B fluid.

2.5. Williamson Fluid

As the last model, we took the Williamson fluid [15] with the equation of motion of

$${\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{\partial }{\partial y}}\left[{\displaystyle \frac{\mu }{\rho }}{\displaystyle \frac{\partial u}{\partial y}}+\mu {\displaystyle \frac{\Gamma }{2}}{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2\right]+\sigma B_0^{2}u.$$

(42)

We cannot present analytical results even if this model (again due to the third-order partial spatial derivative term) lacks the property of time-dependent self-similar symmetry. Malik et al. [38], however, numerically studied the effects of variable thermal conductivity and heat generation/absorption on Williamson fluid flow and heat transfer. The swimming effects of microbes in the blood flow of nano-bioconvective Williamson’s fluid were investigated by Rana et al. [39]. We found these recent results, which are worth mentioning.

3. Summary and Outlook

In our study, we investigated five different time-dependent non-Newtonian boundary layer problems using the self-similar approach. We show that analytical results exist and give them. For the non-Newtonian power-law fluid, we found different types of solutions depending on the value of the power-law exponent, which can be expressed by Kummer’s functions or other implicit functions. The effects of density and other parameters in the nonlinear ordinary differential equation have been investigated. The effects of the magnetic field can be described by including an additional magnetic term in the equation. Analytical solutions have also been found for such cases. In some cases, only implicit solutions can be derived, which is unusual in our experience. For non-Newtonian Casson fluid flow, the effect of the Casson parameter is also analyzed. For the other three non-Newtonian cases, the Oldroyd-B model, Walter’s Liquid B model, and Williamson fluid, the self-similar transformation (5) cannot be performed; contradictions arose among the self-similar exponents. This is not usual for dissipative systems and is due to the third-order spatial derivative term. All five models might be investigated with the traveling wave Ansatz in the far future because the spatial and temporal symmetry shift is always present in these equations. The equations have no explicit temporal or spatial dependence. As an additional not-so-well-known reduction function, we may mention the traveling profile [40] Ansatz, which interpolates between the self-similar and the traveling wave properties.

4. Conclusions

This paper investigates the time-dependent non-Newtonian boundary layer equations with and without the influence of a magnetic field. In the absence of a magnetic field, the horizontal velocity’s shape function can be represented using Kummer functions, whose properties significantly influence the system’s dynamics. When a magnetic field is introduced, the velocity component can typically still be expressed using Kummer functions, although with parameters modified by the field’s presence. In certain cases—depending on the specific non-Newtonian characteristics—simpler functional forms may emerge. Other dynamical variables are either directly related to this velocity component or possess more straightforward expressions.