Numerical Solution of the Newtonian Plane Couette Flow with Linear Dynamic Wall Slip

Abstract

An efficient numerical approach based on weighted-average finite differences is used to solve the Newtonian plane Couette flow with wall slip, obeying a dynamic slip law that generalizes the Navier slip law with the inclusion of a relaxation term. Slip is exhibited only along the fixed lower plate, and the motion is triggered by the motion of the upper plate. Three different cases are considered for the motion of the moving plate, i.e., constant speed, oscillating speed, and a single-period sinusoidal speed. The velocity and the volumetric flow rate are calculated in all cases and comparisons are made with the results of other methods and available results in the literature. The numerical outcomes confirm the damping with time and the lagging effects arising from the Navier and dynamic wall slip conditions and demonstrate the hysteretic behavior of the slip velocity in following the harmonic boundary motion.

1. Introduction

In the past few decades, there has been a growing interest in the study of Newtonian and non-Newtonian viscometric flows in the presence of static or dynamic wall slip conditions for their importance in rheometry and in industrial applications [1,2,3]. Reviews of the slip conditions prevailing at the fluid–structure interface for various media of practical importance have been reported by Hatzikiriakos [2,3]. Malkin and Patlazan [4] also reviewed wall slip in complex fluids of different types, focusing on fluid–wall interaction and shear-induced fluid-to-solid transitions.

The simplest dynamic wall slip equation in the case of unidirectional one-dimensional Newtonian flow, such that the x-velocity component is $v=v(y,t)$ and the plane $y=0$ represents a wall, reads as follows [1]:

$$v_s(0,t)+\lambda {\displaystyle \frac{\partial v_s}{\partial t}}(0,t)={\displaystyle \frac{\mu }{\beta }}\left|{\displaystyle \frac{\partial v}{\partial y}}(0,t)\right|,$$

(1)

where $v_s$ is the slip velocity, defined as the relative velocity of the fluid particles adjacent to the wall with respect to that of the wall, $\mu$ is the constant fluid viscosity, $\beta$ is the slip coefficient, and $\lambda$ is the slip relaxation parameter. When the latter parameter vanishes, Equation (1) is reduced to the Navier slip condition [5]:

$$v_s(0,t)={\displaystyle \frac{\mu }{\beta }}\left|{\displaystyle \frac{\partial v}{\partial y}}(0,t)\right|.$$

(2)

If the wall is fixed, then $v_s$ = v and Equation (1) can be written as follows:

$$v(0,t)+{\lambda }_s{\displaystyle \frac{\partial v}{\partial t}}(0,t)={\displaystyle \frac{\mu }{\beta }}\left|{\displaystyle \frac{\partial v}{\partial y}}(0,t)\right|,$$

Abbatiello et al. [6] presented the mathematical analysis of Navier–Stokes-like problems involving a boundary where dynamic slip applies. Ferrás et al. [7] presented analytical and semi-analytical solutions to some linear and nonlinear problems for Couette and Poiseuille flows for Newtonian and non-Newtonian media with slip boundary conditions at different walls. Thalakkottor and Mohseni [8] used molecular dynamic simulations to study slip at the fluid–solid boundary in an unsteady flow based on Stokes’ second problem, when the wall undergoes an oscillatory motion. They showed the existence of dynamic wall slip and discussed the resulting hysteresis phenomena. Hysteresis was attributed to the unsteady inertial forces of the fluid. Kaoullas and Georgiou [9] derived analytical solutions to some Poiseuille and Couette problems including dynamic wall slip and discussed its effects on the flow development. More recent work on the flow of power law fluids in circular annuli and analytical approximate solutions was presented by Deterre et al. [10]. Ali et al. [11] solved the axial, annular Couette flow of a Newtonian viscous fluid of constant density, taking into account both Navier and dynamic slip boundary conditions, using the Laplace transform technique and inversion by Laguerre polynomials. Farragui et al. [12] used separation of variables to derive analytical solutions to the problem of cessation of annular Poiseuille and Couette flows of a Newtonian fluid with dynamic wall slip.

The analytical solution of flows with dynamic wall slip is possible only for linear problems, e.g., for Newtonian flows with a linear slip equation, such as Equation (1). However, the visualisation of the solution and the analysis of the flow still require numerical calculations, which may not be trivial (see, e.g., [13]). Moreover, the slip equation may be nonlinear [2,3]. Hence, it is necessary to use numerical methods to approximate the solution of this model. These approximation techniques require great effort. In the literature, there exist many methods used to numerically solve partial differential equations: spectral methods [14]; finite element [15]; finite difference, e.g., the weighted average finite difference method (WAFDM) [16,17]; and the collocation method [18,19].

The above literature review clearly shows the importance of using numerical schemes for the solution of Couette and Poiseuille flows of power-law fluids, side by side with the possibility of obtaining exact or approximate analytical solutions.

In the present work, we use an efficient WAFDM scheme to solve the Newtonian planar Couette flow when static or dynamic wall slip applies along the fixed plate and the other plate moves either at constant or oscillatory speed. The evolution of velocity and the volumetric flow rate are calculated for three cases: constant, oscillating, and single-period sinusoidal plate velocity. For the last two cases, the hysteretic behavior of the fluid in following the motion of the wall is put in evidence, as this is implies energy dissipation.

2. Governing Equations

Consider the time-dependent Couette flow of a viscous fluid between infinite parallel walls placed a distance d apart, as illustrated in Figure 1. The fluid is initially at rest and suddenly the upper wall starts moving with velocity $V_{0}f\left(t\right)$ in the x direction. The x-momentum equation takes the form

$$\rho {\displaystyle \frac{\partial v}{\partial t}}=\mu {\displaystyle \frac{{\partial }^{2}v}{\partial y^2}},t>0,0\le y\le d,$$

(3)

where $\rho$ is the constant fluid density. Regarding the boundary conditions, the dynamic slip law applies at the fixed wall (y = 0) and no slip is assumed at the upper wall (y = d). Hence, the boundary conditions of the problem read

$$v(0,t)+{\lambda }_s{\displaystyle \frac{\partial v}{\partial t}}(0,t)={\displaystyle \frac{\mu }{\beta }}\left|{\displaystyle \frac{\partial v}{\partial y}}(0,t)\right|,t>0$$

(4)

and

$$v(d,t)=V_{0}f\left(t\right),t>0.$$

(5)

The initial condition is

$$v(y,0)=0,0\le y\le d.$$

(6)

We work with the dimensionless equations scaling lengths by d, time by $d^2\rho /\mu$, and the velocity by the characteristic velocity $V_0$ of the moving wall. For the sake of simplicity, we keep the same symbols for the dimensionless variables. Hence, the dimensionless form of the initial boundary value problem defined by Equations (3)–(6) is

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial v}{\partial t}}={\displaystyle \frac{{\partial }^{2}v}{\partial y^2}},t>0,0\le y\le 1,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & v(0,t)+{\Lambda }_s{\displaystyle \frac{\partial v}{\partial t}}(0,t)={\displaystyle \frac{1}{B}}{\displaystyle \frac{\partial v}{\partial y}}(0,t),t>0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & v(1,t)=f\left(t\right),t>0\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & v(y,0)=0,0\le y\le 1\hfill \end{array}$$

(10)

where

$$B\equiv {\displaystyle \frac{\beta d}{\mu }}$$

(11)

is the dimensionless slip number, and

$$\Lambda \equiv {\displaystyle \frac{\lambda \mu }{\rho }}$$

(12)

is the dimensionless slip relaxation number. It should be noted that when ${\Lambda }_s=0$, the slip equation is reduced to the static Navier slip equation. Moreover, when $B\to \infty$, the no-slip condition is recovered.

3. Weighted-Average Finite Difference Method

Here, we use a WAFDM scheme [20] to simulate and study the behavior of solutions for problems (7)–(10) of planar Couette flow with Navier and dynamic slip conditions prevailing at the fixed wall only. The motion is triggered by assigning a given motion at the other boundary in its own plane.

The WAFDM is a widely used numerical technique for solving differential equations by approximating the derivatives of a function at discrete points in space and time using finite differences. The discretized solution is eventually obtained by solving a system of algebraic equations [16]. The WAFDM is relatively simple to implement, computationally efficient, and can be applied to a wide range of problems [20]. Its formulation for the problem of Equations (7)–(10) is outlined below. The domain $[0,1]\times [0,T]$ in the $(y,t)$-plane is discretized by a uniform grid with steps $h=\Delta y$ and $k=\Delta t$ so that

$$h={\displaystyle \frac{1}{N}},k={\displaystyle \frac{T}{M}},$$

(13)

where N and M are the numbers of subintervals used for y and t, respectively. Hence, the coordinates of the grid points are

$$y_n=n\Delta y,n=0,1,2,\dots ,N,t_m=m\Delta t,m=0,1,2,\dots ,M.$$

The numerical values of the variable v and the function f at the general grid point $(y_n,t_m)$ are denoted, respectively, by $v_n^m$ and $f_n^m$. The following difference approximations are used for the time and spatial derivatives of the problem:

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{\partial v}{\partial t}}\right)}_n^m& ={\displaystyle \frac{v_n^{m+1}-v_n^m}{k}}+O\left(k\right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{\partial v}{\partial y}}\right)}_n^m& ={\displaystyle \frac{v_{n+1}^m-v_n^m}{h}}+O\left(h\right),\hfill \end{array}$$

(15)

and

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{{\partial }^{2}v}{\partial y^2}}\right)}_n^m& ={\displaystyle \frac{v_{n+1}^m-2v_n^m+v_{n-1}^m}{h^2}}+O\left(h^2\right).\hfill \end{array}$$

(16)

With the WAFDM, the spatial derivatives of v in the RHS of Equations (7) and (8) at time step $m+1$ are approximated as linear combinations of their corresponding values at time steps m and $m+1$ by means of a weight parameter $\theta$, $0\le \theta \le 1$. Substituting Equations (14)–(16) into the governing Equations (7)–(10) leads to a linear system of equations for the unknowns $v_n^m,n=0,1,2,\cdots ,N,m=0,1,2,\cdots ,M$:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{v_n^{m+1}-v_n^m}{k}}=\theta {\displaystyle \frac{v_{n+1}^m-2v_n^m+v_{n-1}^m}{h^2}}+(1-\theta ){\displaystyle \frac{v_{n+1}^{m+1}-2v_n^{m+1}+v_{n-1}^{m+1}}{h^2}},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & \theta v_0^m+(1-\theta )v_0^{m+1}+{\Lambda }_s{\displaystyle \frac{v_0^{m+1}-v_0^m}{k}}={\displaystyle \frac{\theta }{B}}{\displaystyle \frac{v_1^m-v_0^m}{h}}+{\displaystyle \frac{1-\theta }{B}}{\displaystyle \frac{v_1^{m+1}-v_0^{m+1}}{h}},m=0,1,2,\cdots ,M\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & v_N^m=f\left(mk\right),m=1,2,3,\cdots ,M\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & v_0^n=0,n=0,1,2,\cdots ,N.\hfill \end{array}$$

(20)

Depending on the value of the weight factor $\theta$, the method can be explicit ($\theta$ = 1; easy and simple for coding and can be used when the function being approximated is relatively smooth and well behaved) or implicit ($\theta$ = 0; more accurate, has a larger stability region, and can be used when the function being approximated is not smooth). When $\theta =0.5$, the Crank–Nicolson implicit scheme is recovered [20]. More information on the stability and convergence of the numerical method can be found in Ref. [20].

After some manipulations and simplifications, the system of Equations (17)–(19) may be cast in matrix form as follows:

$$\mathbf{A}{\mathbf{V}}^{m+1}=\mathbf{B}{\mathbf{V}}^m+{\mathbf{F}}^m,$$

(21)

where ${\mathbf{V}}^{m+1}$ is the vector of unknown function values at time $m+1$,

$$\mathbf{A}={\left(\begin{array}{ccccccc}(1-\theta )(1+{\displaystyle \frac{1}{Bh}})+{\displaystyle \frac{{\Lambda }_s}{k}}& -{\displaystyle \frac{1-\theta }{Bh}}& 0& 0& \cdots & 0& 0\\ & & & & & & \\ a& b& a& 0& \cdots & 0& 0\\ 0& a& b& a& \cdots & 0& 0\\ 0& 0& a& b& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & b& a\\ 0& 0& 0& 0& \cdots & 0& \alpha \end{array}\right)}_{N+1},$$

(22)

$$\mathbf{B}={\left(\begin{array}{ccccccc}-\theta (1+{\displaystyle \frac{1}{Bh}})+{\displaystyle \frac{{\Lambda }_s}{k}}& {\displaystyle \frac{\theta }{Bh}}& 0& 0& \cdots & 0& 0\\ & & & & & & \\ a^{\prime }& b^{\prime }& a^{\prime }& 0& \cdots & 0& 0\\ 0& a^{\prime }& b^{\prime }& a^{\prime }& \cdots & 0& 0\\ 0& 0& a^{\prime }& b^{\prime }& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & b^{\prime }& a^{\prime }\\ 0& 0& 0& 0& \cdots & 0& 0\end{array}\right)}_{N+1},$$

(23)

and

$${\mathbf{F}}^m={(0,0,0,0,\dots ,f_n^m)}_{N+1}^T,$$

(24)

where $f_n^m\equiv f(x_n,t_m)$,

$$a=-(1-\theta )\xi ,b=1+2(1-\theta )\xi ,\xi ={\displaystyle \frac{k}{h}}$$

(25)

and

$$a^{\prime }=\theta \xi ,b^{\prime }=1-2\theta \xi .$$

(26)

The algorithmic steps of the proposed method are summarized in Algorithm 1. The system (21) is easily solved. The local truncation error of the scheme is of order $O(h+k)$. This scheme is conditionally stable when $\theta >0.5$, and it is unconditionally stable when $\theta \le 0.5$ [20].

Algorithm 1. Algorithmic steps of the proposed method

Initialize parameters and finite difference grid

Calculate initial condition $V^0$

Build A

% Time marching loop

for m = 1, …, M

% Node loop

for n = 1, …, N − 1

Build $B{V}^{m-1}+F^{m-1}$

end

Solve Equation (21) for $V^m$

end

Visualize the results

4. Results and Discussion

In this section, we use the introduced numerical scheme (17)–(20) to simulate the approximation solution of (7)–(10). It will be shown that the introduced WAFDM provides good approximations for the solution of (7)–(10). It is applicable and efficient for solving the given system of governing equations with the accompanying boundary and initial conditions. Different test examples are carried out to approximate the solutions and compare between the solutions using different techniques. We compared the computational time of the WAFDM and the Legendre spectral collocation method and found that the WAFDM is more efficient.

The results assess the effect of the various parameters in the initial and boundary conditions to which the solution is subjected, demonstrate the efficiency of the introduced technique, and justify the accuracy in comparison with other methods. In our numerical experiments we choose different values for the weight factor and full agreement is reached with the theoretical stability condition.

4.1. Case 1: Plate Moving at a Constant Speed

In order to assess the efficiency of the proposed numerical scheme, we have compared our results (WAFDM) for the special case when $\theta =0$ and $f\left(t\right)=1$ with those obtained using the MATLAB2021b pdepe toolbox. This is one of the cases treated in [21], where one plate is kept fixed, while the other one is suddenly set to motion from rest with constant speed. In Figure 2, the predictions of pdepe and the WAFDM for ${\Lambda }_s$ = 0 (Navier slip) and three values of the slip number B, i.e., B = 0.1 (strong), 1 (moderate), and 10 (weak slip), are compared. The agreement is excellent and this is perhaps more visible when comparing the slip velocities $v(0,t)$ in Figure 3. It should be noted that full agreement is also reached with the results presented in [21] based on a Fourier expansion and on the use of one-sided Fourier transform. The main feature of the solution is the reduction in the slip velocity at the boundary $y=0$ as parameter B increases. The efficiency on the computational time of the WAFDM is demonstrated by comparing it with that of the Legendre spectral collocation method depending on the Legendre polynomial of degree M in

Table 1. Comparison of the WAFDM with the spectral collocation method.

WAFDM		Spectral Method
N	CPU	M	CPU
500	5.9 s	5	963.6 s
5000	18.3 s	12	7579.9 s
10,000	44.8 s	15	88,937.5 s
50,000	789.9 s	20	219,357.6 s

.

The volumetric flow rate, defined by

$$Q\left(t\right)={\int }_0^{1}v(y,t)dy,$$

(27)

is shown in Figure 4 for $\theta =0$, B = 0.1 (strong slip) and 1 (moderate slip), and  ${\Lambda }_s$ = 0 (Navier slip), 1, and 10. The damping effect of the slip relaxation parameter in reaching a fully developed flow is clearly shown. As expected, the initial value of the volumetric flow rate is independent of ${\Lambda }_s$ and decreases with the slip parameter B (since wall slip is reduced).

4.2. Case 2: Oscillating Plate

Consider now the case where the moving plate is performing continuous harmonic oscillations with period equal to unity, i.e.,

$$f\left(t\right)=sin4\pi t.$$

(28)

A similar problem was addressed in [8] for unsteady flow to estimate the increment in slip at the boundary due to wall acceleration and uncover the hysteretic behavior of the slip velocity for this harmonic motion of the boundary. Here, it is required to evaluate the effect of the dynamic slip parameter ${\Lambda }_s$ on the flow in general, and on the hysteretic behavior of the slip velocity due to oscillating wall motion.

Figure 5 illustrates the distribution of the solution in space and time with different values of ${\Lambda }_s$, as time is taken to run along two complete periods of the wall oscillations. The boundary condition at the moving wall is clearly satisfied. As ${\Lambda }_s$ increases from $0.1$ to 10, a damping of the amplitude of the oscillations with time takes place. The differences become more noticeable as one approaches the fixed wall at $y=0$, where the velocity profile becomes more flattened.

Figure 6 illustrates the solutions at $y=0$, i.e., the slip velocity, as functions of time for different values of ${\Lambda }_s$. The figure clearly shows time damping and lagging effects at the boundary $y=0$. The amplitudes of oscillations in the slip velocity show a reduction of approximately $33\%$ as the value of ${\Lambda }_s$ increases from $0.1$ to 10.

In order to put in evidence the crucial role played by the weight factor $\theta$ in the present numerical study, we have considered one case with a new value of this parameter to test the stability of the numerical scheme. As mentioned above, the scheme is stable provided that $\theta <0.5$ [20]. The instability occuring when $\theta =0.7>0.5$ is illustrated in the 3D plot in Figure 7.

The 3D plots in Figure 8 show how the solutions change in space and time when the slip parameter B assumes the values 0.1, 1, and 10. Figure 9 illustrates the dependence of the slip velocity on parameter B. In both plots, we have fixed the value ${\Lambda }_s=0.5$. The effect of parameter B becomes more pronounced as the boundary $y=0$ is approached. Here, again, one notes the presence of lag in following the boundary motion, and time damping of the peaks at the boundary $y=0$ as the value of B increases.

In order to visualize the hysteretic behavior of the slip velocity and the lagging following the applied motion on one wall noticed in [8], plots have been provided in Figure 10 of the slip velocity against the boundary motion $sin\omega t$ for three values of the dynamic slip parameter ${\Lambda }_s$. It is noticed that the hysteresis loops become narrower as the value of the dynamic slip parameter increases. In [8], the width of the hysteresis loop is related to the loss of energy transfer from the wall to the fluid.

The volumetric flow rate is shown in Figure 11, where it is seen that this is oscillatory and is damped by the slip relaxation parameter. However, for sufficiently large values of ${\Lambda }_s$, it is seen that the amplitude of oscillations in $Q\left(t\right)$ will be less sensitive to any further increase in ${\Lambda }_s$. This phenomenon becomes more striking for larger values of the slip parameter B. Notice that after two complete oscillations, the volumetric flow rate has not vanished due to retardation.

4.3. Case 3: Single-Plate Oscillation

In this section, we assume that the upper plate oscillates only once and then comes to rest, i.e.,

$$f\left(t\right)=\left\{\begin{array}{cc}sin4\pi t,\hfill & t\le 0.5\hfill \\ 0,\hfill & t>0.5\hfill \end{array}\right..$$

(29)

Figure 12 illustrates the effect of the dynamic slip parameter ${\Lambda }_s$ on the solution, where the value of the slip parameter was set to $B=1$.

Figure 13 shows the effect of the relaxation parameter ${\Lambda }_s$ on the slip velocity $v(0,t)$ for $B=1$ (moderate slip). As for the case of continuous harmonic boundary motion, the effect of the dynamic slip parameter becomes stronger as the boundary $y=0$ is approached as can be seen in Figure 12. The same remains valid for the effect of slip parameter B on the fluid motion as shown in Figure 14 and Figure 15, where the value of the dynamic slip parameter was set ${\Lambda }_s=0.5$.

Figure 16 represents the hysteretic behavior of the slip velocity in following the motion of the wall in Case 3 for $B=1$ and three values of the dynamic slip parameter ${\Lambda }_s$. Here, again, it is noticed that the hysteresis loop becomes narrower as the value of ${\Lambda }_s$ increases.

Finally, we show in Figure 17 the volumetric flow rate in Case 3 with $\theta =0$ for two values of B and three values of ${\Lambda }_s$. The volumetric flow rate is almost independent of these two slip parameters.

5. Conclusions

A weighted-average finite difference scheme has been used to solve the time-dependent (start-up) plane Couette flow with dynamic slip along the fixed plate. Three different cases have been considered for the motion of the moving plate, i.e., constant-speed, sinusoidal, and single oscillation. The numerical solutions compare well with available analytical and numerical solutions. Both the slip and relaxation parameters appear to decelerate the evolution of the flow. For the two cases of accelerated (sinusoidal and one-period sinusoidal) boundary motion, the numerical results have clearly demonstrated a hysteretic behavior of the slip velocity that will be responsible for time lag and loss of energy transfer by the moving wall to the fluid.

Future work will be devoted to extending the numerical scheme for solving nonlinear extensions of the problem studied here. These include the non-Newtonian (e.g., power-law) flow under both Navier and dynamic slip at the fixed wall and the Newtonian flow with slip obeying a nonlinear dynamic slip law.