Unified Theory of Unsteady Planar Laminar Flow in the Presence of Arbitrary Pressure Gradients and Boundary Movement

Abstract

A general unified solution of the plane Couette–Poiseuille–Stokes–Womersley incompressible linear fluid flow in a slit in the presence of oscillatory pressure gradients with periodic synchronous vibrating boundaries is presented. Oscillatory flow remains stable and laminar with no-slip boundary conditions applied. Eigenfunction expansion method is used to obtain the exact analytical solution of the general linear inhomogeneous boundary value problem. Fourier expansion of arbitrary harmonic pressure gradient and non-harmonic wall oscillations was used to calculate arbitrary driving of the fluid. In-house developed optimized computational fluid dynamics marching-in-time finite-volume method was used to test and verify all analytical results. A number of particular transients, steady-state and combined flows were obtained from the general analytical result. Generalized Stokes and Womersley flows were solved using the analytical computations and numerical experiments. The combined effects of periodic non-harmonic wall movements with oscillatory pressure gradients offers rich and interesting flow patterns even for a linear Newtonian fluid and may be particularly interesting for pumping-assist microfluidic devices. The main motivation for developing a unified solution of the unsteady laminar planar Couette–Stokes–Poiseuille–Womersley flow originates in a need for, but is not limited to, in-depth exploration of flow patterns in hemodynamic and microfluidic pumping applications.

1. Introduction

The solutions of the classical fully-developed steady laminar flows, such as the planar Couette (CT) and the planar Poiseuille (PO) flows in an infinite slit, have been known for a long time [1,2,3,4,5,6,7,8,9]. The Couette flow originates in a fluid motion between two coaxial cylinders and is shear-driven (wall-driven). If the radii of the cylinders are large, then planar CT flow is approximated. On the other hand, PO flow is solely pressure-driven. A planar steady-state Couette–Poiseuille (CP) flow is a linear combination of the CT flow with one boundary moving at uniform speed and the PO flow between two fixed infinite plates caused by the steady-state (proverse or adverse) pressure gradients (PG). These classic solutions of the formidable nonlinear Navier–Stokes (N-S) equations also paved the way to early investigation of the flow stability and transition to turbulence [2,5,10]. The Womersley (WO) flow that can be seen as a generalization of the PO flow is pressure-driven caused by oscillating PGs with fixed conduit boundaries. Such flows have also been known for quite some time [4,5,7,11,12,13,14,15] and often under different names too. WO-type flow has important applications in pulsatile arterial hemodynamics and many biological fluid systems [15,16,17].

Unsteady shear-driven flows with moving boundaries were first tackled with, so-called, Stokes’ (ST) 1st and 2nd problem [1,4,5]. Stokes’ 1st problem (ST1) treats flow transients after a single infinite flat plate is suddenly set in uniform motion. A fluid then asymptotically settles into the steady-state (Couette-like) flow after decay of initial transients. Stokes’ 2nd problem (ST2) treats the harmonically vibrating infinite flat plate problem that sets still semi-infinite fluid in oscillatory motion. The quasi-steady-state (QSS) flow oscillations take place after the initial transients are attenuated. A fluid closer to the moving boundary will be more affected by the oscillating wall due to viscous diffusivity, which is chiefly responsible for the decaying transfer of the linear momentum deeper into the fluid mass [17,18,19]. Stokes’ 2nd problem (ST2) is in many respects similar to Kelvin’s problem of penetration of harmonically induced temperature waves into the semi-infinite (soil) medium. Carslaw and Jaeger [20] provide detailed analytical solutions and discussion of this important problem. Thus, ST (ST1 and ST2) and WO flows are by nature unsteady, while CT and PO flows are steady flows.

Numerous experiments have confirmed the validity and theoretical soundness of these particular basic flows. The general validity of the N-S equations in continuum hydrodynamic approximation has also been confirmed in numerous experiments. The question of the molecular, apparent, or effective slip at the liquid–solid boundary is still open and a topic of intense current research activities. Although based on an entirely different concept, the Lattice Boltzmann Method (LBM) originating in mesoscopic consideration of the particle density function on symmetry-preserving lattice including streaming and collision operators have also confirmed the general validity of the incompressible N-S equations in a low-Mach number LBM limit [21,22,23,24]. The continuum fluid behavior has also been confirmed by the Molecular Dynamic (MD) and Metropolis Monte Carlo (MMC) simulations on a true microscopic (particle) scale involving Hamiltonian of atomic and molecular interactions [25,26,27,28,29]. At very small length scales, microscales for gases and nanoscales for liquids, we do however expect departure from the continuum N-S equations [30,31].

The basic flow models mentioned above have been already used on many occasions in modeling hemodynamic flow and cardiovascular circulation. Early measurements by McDonald [32] and Hale et al. [33] were supported by theoretical considerations of oscillatory pressure-driven flows in elastic tubes by Womersley [34,35]. A good review of blood flow phenomena in arteries was presented in [36]. Daidzic [17] successfully used the Womersley flow model to smooth out velocity profiles and recover oscillatory shear stress information lost during multi-gate Doppler ultrasound measurements of blood mimicking fluid in blood-vessel-type conduits. Additionally, very recently, Giusti and Mainardi [37] studied viscoelastic fluid effects in fluid-filled elastic tubes. Although it was suspected for a long time, it was actually only quite lately shown experimentally that the longitudinal motion of blood vessel walls exists due to wall shear [38] on the blood vessel interface.

The main goal of this research article is to present a unified theory of laminar Newtonian fluid flows with moving boundaries (shear-driven) and unsteady PGs (pressure-driven) for the combined planar CT, PO, ST, and WO, or as termed here, a planar Couette–Poiseuille–Stokes–Womersley (CPSW) flow. Flow is generated in planar (slit) geometry between two infinite rigid impenetrable flat plates. Since the walls are parallel and impermeable only the axial velocity component remains. Interaction of forced wall oscillations and applied unsteady PGs results in interesting transient and QSS flow patterns. Specifically, non-harmonic shear and/or pressure forces for linear fluid are discussed for the first time in a consistent and in-depth manner to the best of our knowledge. In many instances the hemodynamic flows can be approximated with a Newtonian fluid and so this theoretical numerical analysis is directly applicable. One of the important applications of the findings from this study is in the utilization of shear-driven flows and the design of micropumps.

2. Mathematical Model

Within the realm of continuum hydrodynamics, the N-S equations will be used as a starting point to describe an unsteady motion of the Newtonian fluid. The flow is incompressible and divergence-free. The oscillatory non-harmonic synchro-phased periodic motion of two solid interfaces in a slit geometry will be modeled using the Fourier expansion of the arbitrary periodic excitation waveforms. The PGs are modeled by assuming an infinite pressure disturbances group speed and will thus be time-dependent only [8].

We made few important assumptions: fluid behavior is incompressible, linear, and time-independent (no internal dynamics) where viscosity is shear- and temperature-independent, flow remains laminar and isothermal at all times, entrance effects are neglected, boundary walls are rigid and impenetrable, Dirichlet no-slip BCs exist, periodic oscillations of boundary movement and PGs and hyperbolic effects are neglected. The IBVP discussed here is well-posed [39,40].

2.1. Fluid Flow Model

The general 3D unsteady conservation equations for mass, momentum, and energy for a linear compressible fluid in conventional tensor and vector notation and Cartesian coordinates [2,5] yields:

$$\begin{array}{l}\frac{1}{\rho }\frac{D\rho }{Dt}=-\frac{\partial u_j}{\partial x_j}\hspace{1em}\left(=-div\hspace{0.17em}\overrightarrow{u}=-\nabla \hspace{0.17em}\overrightarrow{u}\right)\\ \rho \frac{D{u}_i}{Dt}=f_i+\frac{\partial {\sigma }_{ij}}{\partial x_j}\\ \rho c_p\frac{DT}{Dt}=\frac{\partial }{\partial x_j}\left[k\left(T\right)\frac{\partial T}{\partial x_j}\right]+\Phi \hspace{0.17em}+\beta T\frac{Dp}{Dt}\\ p=p\left(\rho ,T\right)\\ x_i=(x,y,z)\hspace{1em}\hspace{1em}{u}_i=(w,v,u)\end{array}$$

(1)

We can regard Equation (1) as a sort of “standard model” of continuum (low-Knudsen number) fluid mechanics. The N-S equations and the thermal energy equation are essentially phenomenological as the linear fluid stress–strain rate and the heat conduction relationships are semi-empirical. Generally, for Newtonian (linear) fluids, we may write:

$$\begin{array}{l}{\sigma }_{ij}=-p{\delta }_{ij}+{\tau }_{ij}\hspace{1em}\hspace{1em}\hspace{1em}{\tau }_{ij}=2\mu D_{ij}+\lambda \frac{\partial u_k}{\partial x_k}{\delta }_{ij}=2\mu \left(D_{ij}-\frac{1}{3}{\delta }_{ij}\frac{\partial u_k}{\partial x_k}\right)+\xi {\delta }_{ij}\frac{\partial u_k}{\partial x_k}\\ D_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)=\frac{1}{2}\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)=\frac{1}{2}\left(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T\right)\\ \xi =\left(\lambda +\frac{2}{3}\mu \right)=0\hspace{1em}\Rightarrow \hspace{1em}\lambda =-\frac{2}{3}\mu \hspace{1em}\left(\mathrm{Stokes}’ \mathrm{Hypothesis}\right)\end{array}$$

(2)

where for conservative forces:

$$f_i=\nabla \psi =\frac{\partial \psi }{\partial x_i}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{\partial p}{\partial x_i}+f_i=-\frac{\partial }{\partial x_i}\left(p+\psi \right)$$

and:

$$\beta =-\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial T}\right)}_{p=const}\hspace{1em}\hspace{1em}\mu =\mathrm{const}.,\hspace{0.17em}{c}_p=\mathrm{const}.,\hspace{0.17em}k=\mathrm{const}.$$

The viscous dissipation function is a non-negative expression of the irreversible flow energy conversion into heat due to deviatoric part of the stress tensor [2,5]:

$$\Phi ={\tau }_{ij}\frac{\partial u_i}{\partial x_j}=\left(2\mu D_{ij}+\lambda \frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)\frac{\partial u_i}{\partial x_j}$$

(3)

The vorticity, which is the curl of vorticity and twice the angular rotation speed, is [2,5]:

$${\omega }_k=-{\epsilon }_{ijk}\frac{\partial u_i}{\partial x_j}=-{\epsilon }_{ijk}{\partial }_j{u}_i$$

(4)

Typically, it is assumed that the bulk viscosity $\xi$ vanishes (Stokes hypothesis). Additionally, if the flow is incompressible, the second-viscosity term $\lambda$ vanishes identically and all is left from the deviatoric part of the stress tensor is the deformation tensor $D_{ij}$.

A Newtonian incompressible fluid in an unsteady laminar fully developed flow $\left(w,v,u\right)$ with oscillatory PG components and both walls vibrating periodically in plane slit geometry $\left(x,y,z\right)$ is investigated. Both the PGs and the wall motion may oscillate in an arbitrary periodic fashion (ST2). A schematic of the flow geometry is illustrated in Figure 1. Other than gravity, no other body forces are assumed to be acting.

The unsteady fully developed isothermal laminar flow in an infinite slit geometry with no-slip BCs and incompressible Newtonian (linear) fluid can thus be reduced to the following IBVP:

$$\begin{array}{l}\frac{\partial u}{\partial z}=0\hspace{1em}\Rightarrow \hspace{1em}u=u\left(y,t\right)\\ \rho \frac{\partial u}{\partial t}=\mu \frac{{\partial }^{2}u}{\partial y^2}-\frac{\partial p}{\partial z}\left(t\right)+\rho \cdot g\cdot \cos \gamma \\ 0=-\frac{\partial p}{\partial y}+\rho \cdot g\cdot \sin \gamma \end{array}$$

(5)

The auxiliary conditions are:

$$\begin{array}{l}\mathrm{BC}:\hspace{1em}u\left(+h,t\right)=G_{ss}+g\left(t\right)=G\left(t\right)={\epsilon }_{ss}^u{U}_w+{\displaystyle \sum _{k=-\infty }^{\infty }{\epsilon }_{{}_k}^u{U}_w\exp \left(ik{\Omega }_0^{u}t\right)}\\ \mathrm{BC}:\hspace{1em}u\left(-h,t\right)=R_{ss}+r\left(t\right)=R\left(t\right)={\epsilon }_{ss}^l{U}_w+{\displaystyle \sum _{m=-\infty }^{\infty }{\epsilon }_m^l{U}_w\exp \left(im{\Omega }_0^{l}t\right)}\\ \mathrm{IC}:\hspace{1em}u\left(y,0\right)=z\left(y\right)\\ T=\mathrm{const}.\hspace{1em}\hspace{1em}\rho =\mathrm{const}.\hspace{1em}\hspace{1em}\mu =\mathrm{const}.\hspace{1em}\hspace{1em}{\epsilon }_{ss}^u,{\epsilon }_{ss}^l=\mathrm{const}.\hspace{0.17em}\left(\le \ge 0\right)\end{array}$$

The BCs are split into the steady-state part due to uniform wall motion (CT component) and an arbitrary periodic wall oscillation (ST component) described by the appropriate complex Fourier series. Similarly, the PGs consist of the steady-state (PO) component superposed onto a linear combination of many (WO) harmonic components. The phasors of the PGs are complex numbers which contain amplitude (modulus) and phase (argument) information:

$$-\frac{\partial p}{\partial z}\left(t\right)+\rho \cdot g\cdot \cos \gamma ={\left(-\frac{\partial p}{\partial z}\right)}_{ss}+{\displaystyle \sum _{k=-\infty }^{\infty }{P}_k\exp \left(ik{\omega }_{0}t\right)}\hspace{1em}\hspace{1em}{P}_k=\left|P_k\right|\cdot e^{i{\phi }_k}\in C\hspace{1em}\hspace{1em}i=\sqrt{-1}$$

(6)

The steady-state (static) part of the PG includes gravity component:

$${\left(-\frac{\partial p}{\partial z}\right)}_{ss}={\mathrm{\Pi }}_{ss}+\rho \cdot g\cdot \cos \gamma =P_{ss}=\mathrm{const}.\hspace{1em}\hspace{1em}{\mathrm{\Pi }}_{ss}\le \ge 0\hspace{1em}\hspace{1em}{P}_{ss}\le \ge 0\hspace{1em}\hspace{1em}\gamma \in \left[0,\pi \right]$$

Although the real-valued steady-state externally applied PG and the DC-value, or average-value, zero-order Fourier coefficient have different origins they can be lumped together.

Expressions for oscillating PGs and wall movement are conveniently given in a complex form greatly simplifying theoretical analysis. All periodic functions for the PGs and the BCs are real-valued, thus possessing conjugate symmetry, ($C_n^{*}=C_{-n}$) [39,40,41,42]. The complex Fourier series can be converted into the classical or the phase form for an arbitrary function:

$$\begin{array}{l}f\left(t\right)={\displaystyle \sum _{n=-\infty }^{\infty }{C}_n{e}^{in\omega \hspace{0.17em}t}}=\frac{A_0}{2}+{\displaystyle \sum _{n=1}^{\infty }\left[A_n\cos \left(n\omega \hspace{0.17em}t\right)+B_n\sin \left(n\omega \hspace{0.17em}t\right)\right]}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{A_0}{2}+{\displaystyle \sum _{n=1}^{\infty }{D}_n\cos \left(n\omega \hspace{0.17em}t+{\phi }_n\right)}=\frac{A_0}{2}+{\displaystyle \sum _{n=1}^{\infty }{K}_n\sin \left(n\omega \hspace{0.17em}t+{\varphi }_n\right)}\end{array}$$

(7)

where:

$$D_n=K_n=\sqrt{A_n^2+B_n^2}=2\sqrt{C_n{C}_{-n}}=2\left|C_n\right|\hspace{1em}\hspace{1em}{\phi }_n=-{\tan }^{-1}\left(\frac{B_n}{A_n}\right)\hspace{1em}\hspace{1em}{\varphi }_n={\tan }^{-1}\left(\frac{A_n}{B_n}\right)$$

The complex Fourier coefficients are calculated using the inherent orthogonal property of conjugate complex exponentials:

$$\begin{array}{l}{C}_n=\frac{\langle f,{\phi }_n^{*}\rangle }{{\Vert {\phi }_n\Vert }^2}=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}f\left(t\right)\cdot e^{-in\omega \hspace{0.17em}t}dt}\\ {\Vert {\phi }_n\Vert }^2=\langle {\phi }_n,{\phi }_n^{*}\rangle =\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}{e}^{in\omega \hspace{0.17em}t}\cdot e^{-in\omega \hspace{0.17em}t}dt}=T\hspace{1em}\hspace{1em}n=0,\pm 1,\pm 2,\dots \end{array}$$

(8)

It is convenient to express the complex Fourier coefficients as phasors where the modulus is $\left|C_{-n}\right|=\left|C_n^{*}\right|=\left|C_n\right|$, and the phase satisfies ${\varphi }_{-n}=-{\varphi }_n$. The complex Fourier series now yields:

$$f\left(t\right)={\displaystyle \sum _{n=-\infty }^{\infty }{C}_n\hspace{0.17em}{\phi }_n}={\displaystyle \sum _{n=-\infty }^{\infty }\frac{\langle f,{\phi }_n^{*}\rangle }{{\Vert {\phi }_n\Vert }^2}\hspace{0.17em}{\phi }_n}=C_0+{\displaystyle \sum _{n=1}^{\infty }\left[C_n{e}^{in\omega \hspace{0.17em}t}+C_{-n}{e}^{-in\omega \hspace{0.17em}t}\right]}=C_0+{\displaystyle \sum _{n=1}^{\infty }2\left|C_n\right|\cos \left(n\omega \hspace{0.17em}t+{\varphi }_n\right)}$$

(9)

One must be careful when calculating the phase constant. All quadrants must be inspected as the inverse tangent is not a single-valued function. The modulus of the complex Fourier coefficient must be a real non-negative number. The complex Fourier coefficients carry the information from both, the real and the imaginary, part if the function is neither even nor odd ($A_n=2\cdot \Re \left\{C_n\right\},\hspace{0.17em}{B}_n=-2\cdot \Im \left\{C_n\right\}$). For an even real-valued function, $C_n=C_{-n}=C_n^{*}=A_n/2$ ($A_n=A_{-n}$), while for an odd real-valued function, $C_{-n}=-C_n=B_n/2$, $B_{-n}=-B_n$, and, $C_0=A_0=0$. Any function can be represented as the sum of respective even and odd functions. In the case of non-harmonic periodic wall oscillations, the Fourier series will necessarily contain both, real and imaginary, components.

Due to the mass conservation the only velocity component surviving is the longitudinal (axial) along the conduit axis. Due to the divergence-free flow, the second-viscosity term disappears identically from the general strain rate expression (bulk viscosity is zero) without the need to invoke Stokes’ hypothesis. Although we calculate dissipation function, the energy equation is not solved explicitly as the flow is enforced to remain isothermal. Alteration of physical properties due to temperature changes are thus neglected and difficult nonlinear coupled problem is avoided. As the slit thickness is small, the gravity plays no significant role in the lateral direction and the pressure is uniformly distributed in each cross-section. However, the gravitational potential may be important for driving flows in non-horizontal slits.

2.2. Dimensional Form of the IBVP

The main IVBP in compact dimensional form yields:

$$\begin{array}{l}\rho \frac{\partial u}{\partial t}=\mu \frac{{\partial }^{2}u}{\partial y^2}+\Gamma \left(y,t\right)\\ \mathrm{BC}:\hspace{1em}u\left(+h,t\right)=G\left(t\right)\\ \mathrm{BC}:\hspace{1em}u\left(-h,t\right)=R\left(t\right)\\ \mathrm{IC}:\hspace{1em}u\left(y,0\right)=f\left(y\right)\end{array}$$

(10)

The expression for PG source function (only time dependent here but we keep spatial dependence for generality), is:

$$\Gamma \left(y,t\right)=P_{ss}+{\displaystyle \sum _{m=-\infty }^{\infty }{P}_m\exp \left(i{\omega }_{m}t\right)}=P_{ss}+P_0+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|\cos \left(m{\omega }_{0}t+{\varphi }_m\right)}$$

(11)

Additionally, for BCs:

$$\begin{array}{l}G\left(t\right)=G_{ss}+g\left(t\right)={\epsilon }_{ss}^u{U}_w+{\displaystyle \sum _{k=-\infty }^{\infty }{\epsilon }_k^u{U}_w\exp \left(i{\Omega }_k^{u}t\right)}={\epsilon }_{ss}^u{U}_w+{\epsilon }_0^u{U}_w+U_w{\displaystyle \sum _{k=1}^{\infty }2\left|{\epsilon }_k^u\right|\cos \left[\left(k{\Omega }_0^{u}t+{\varphi }_k\right)\right]}\\ R\left(t\right)=R_{ss}+r\left(t\right)={\epsilon }_{ss}^l{U}_w+{\displaystyle \sum _{j=-\infty }^{\infty }{\epsilon }_j^l{U}_w\exp \left(i{\Omega }_j^{l}t\right)={\epsilon }_{ss}^l{U}_w+{\epsilon }_0^l{U}_w+U_w{\displaystyle \sum _{j=1}^{\infty }2\left|{\epsilon }_j^l\right|\cos \left[\left(j{\Omega }_0^{l}t+{\psi }_j\right)\right]}}\end{array}$$

(12)

In the case of the periodic wall motion, there will be no Couette velocity components (${\epsilon }_{ss}^l={\epsilon }_{ss}^u=0$). In the case of an odd wall-driving function, the (Stokes’ component) DC-value will be zero. In the case of an even function, a non-vanishing velocity DC-component can be regarded as a standalone steady-state value. The ICs are general and an arbitrary velocity distribution at time zero can be specified.

2.3. Nondimensional form of the IBVP

Transforming general CPSW flow into dimensionless form is difficult. There are several choices of the characteristic velocity [5,7]. In the case of the classic CT flow, the plate’s uniform velocity is the characteristic velocity and could never be exceeded in the flow field. For a planar PO flow, the maximum velocity in the center of the parabolic profile which is a function of the uniform PG can be used as a characteristic velocity. For the oscillatory PG (WO-component) we could use amplitude and frequency of any harmonic as a characteristic velocity. It is recommended to use the one with the largest magnitude [7] to ensure range [0, 1], but that is not always possible. Let us now define new dimensionless variables:

$$\begin{array}{l}{u}^{*}=\frac{u}{U_0}\hspace{1em}\hspace{1em}{y}^{*}=\frac{y}{h}\hspace{1em}\hspace{1em}{t}^{*}=\frac{t}{t_D}=t\frac{\nu }{h^2}=\frac{t}{h^2/\nu }\hspace{1em}\hspace{1em}{\omega }^{*}=\omega \frac{h^2}{\nu }\hspace{1em}\hspace{1em}{\Omega }^{*}=\Omega \frac{h^2}{\nu }\hspace{1em}\hspace{1em}{t}_D=\frac{h^2}{\nu }\\ P_{ss}^{*}=\frac{P_{ss}\hspace{0.17em}{h}^2}{2\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}{P}_0^{*}=\frac{P_0\hspace{0.17em}{h}^2}{2\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}{\left|P_m\right|}^{*}=\frac{\left|P_m\right|\hspace{0.17em}{h}^2}{2\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}\hspace{1em}{\Gamma }^{*}=\frac{\Gamma \hspace{0.17em}{h}^2}{\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}\hspace{1em}{\tau }^{*}=\frac{\tau \hspace{0.17em}h}{\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}\hspace{1em}{\Phi }^{*}=\frac{\Phi \hspace{0.17em}{h}^2}{\mu \hspace{0.17em}{U}_{ch}^2}\end{array}$$

(13)

The velocity scaling will change depending on the particular flow situation, i.e., presence of the wall movement or steady and/or dynamic PGs. The characteristic velocity could be uniform wall velocity, maximum plane PO velocity in the slit centerline or the mth-harmonic WO velocity caused by respective dynamic PG oscillations:

$$U_{ch}=U_w\hspace{1em}\mathrm{or}\hspace{1em}{U}_{ch}=\frac{P_{ss}\hspace{0.17em}{h}^2}{2\mu }\hspace{1em}\mathrm{or}\hspace{1em}{U}_{ch}=\frac{P_0\hspace{0.17em}{h}^2}{2\mu }\hspace{1em}\mathrm{or}\hspace{1em}{U}_{ch}=\frac{\left|P_m\right|}{\rho \hspace{0.17em}{\omega }_m}$$

(14)

One can now define a Womersley (WO) number which is often also referred to as the “Witzig” number [37] or the “kinetic Reynolds” number (square of the WO number) [2,5,7]:

$${\alpha }_{\omega }=h\sqrt{\frac{\omega }{\nu }}\hspace{1em}\hspace{1em}{\alpha }_{\Omega }=h\sqrt{\frac{\Omega }{\nu }}\hspace{1em}\hspace{1em}{\omega }^{*}={\alpha }_{\omega }^2=h^2\frac{\omega }{\nu }\hspace{1em}\hspace{1em}{\Omega }^{*}={\alpha }_{\Omega }^2=h^2\frac{\Omega }{\nu }$$

(15)

The WO-number, ${\alpha }_n={\alpha }_0\sqrt{n}$, has a very specific physical meaning and can be thought of as the ratio of the half channel height versus Stokes’ or viscous penetration layer. For example, due to wall oscillations the momentum (or more appropriately—vorticity) is being diffused. The ensuing pseudo “shear-waves” are purely dissipative (and dispersive) and attenuate exponentially with distance from the solid vorticity-generating surfaces. Alternatively, we may think of the WO number as the ratio of times due to diffusion and oscillation or as the ratio of velocities due to oscillations and viscous momentum diffusion:

$${\alpha }_{\Omega }=\frac{h}{\sqrt{\nu /\Omega }}\cong \frac{h}{\delta }\hspace{1em}\hspace{1em}{\Omega }^{*}={\alpha }_{\Omega }^2=\frac{h^2/\nu }{1/\Omega }=\frac{t_{diffusion}}{t_{oscillation}}\hspace{1em}\hspace{1em}{\Omega }^{*}={\alpha }_{\Omega }^2=\frac{\Omega \hspace{0.17em}h}{\nu /h}=\frac{u_{oscillation}}{u_{diffusion}}$$

A very high WO number (>10) would in this case mean that the lateral channel dimension is much larger than the zone of influence (radius of action) of the moving surface and that the bulk of the fluid will not “feel” its presence:

$${\alpha }_{\Omega }=\sqrt{{\Omega }^{*}}>>1\hspace{1em}\Rightarrow \hspace{1em}h>>\delta ,\hspace{1em}{t}_{oscillation}<<t_{diffusion},\hspace{1em}{u}_{oscillation}>>u_{diffusion}$$

For small WO numbers, the viscous penetration layer extends practically instantly over the entire lateral flow domain and diffusion is rapid compared to slow wall oscillations. In this case, we have a sort of quasi-steady-state planar PO flow where the fluid follows the slow boundary motion with vanishing phase lag. Hence, we arrive at the quasi-equilibrium solutions:

$${\alpha }_{\Omega }=\sqrt{{\Omega }^{*}}<<1\hspace{1em}\Rightarrow \hspace{1em}h<<\delta ,\hspace{1em}{t}_{oscillation}>>t_{diffusion},\hspace{1em}{u}_{oscillation}<<u_{diffusion}$$

By using the dimensionless variables, defined in Equation (13), and substituting them into Equations ((10)–(12)), we obtain the non-dimensional IBVP:

$$\begin{array}{l}\frac{\partial u^{*}}{\partial t^{*}}=\frac{{\partial }^2{u}^{*}}{\partial y^{*}{}^2}+{\Gamma }^{*}\left(y^{*},t^{*}\right)\\ \mathrm{BC}:\hspace{1em}{u}^{*}\left(+1,t^{*}\right)=G_0^{*}+g^{*}\left(t\right)=G^{*}\left(t^{*}\right)\\ \mathrm{BC}:\hspace{1em}{u}^{*}\left(-1,t^{*}\right)=R_0^{*}+r^{*}\left(t\right)=R^{*}\left(t^{*}\right)\\ \mathrm{IC}:\hspace{1em}{u}^{*}\left(y^{*},0\right)=f^{*}\left(y^{*}\right)\end{array}$$

(16)

with

$$\begin{array}{l}{\Gamma }^{*}\left(y^{*},t^{*}\right)=2P_{ss}^{*}+{\displaystyle \sum _{n=-\infty }^{\infty }2\left|P_n^{*}\right|\exp \left[i\left(n{\omega }_0^{*}{t}^{*}+{\phi }_n\right)\right]}=2P_{ss}^{*}+q^{*}\left(y^{*},t^{*}\right)\\ G^{*}\left(t^{*}\right)={\epsilon }_{ss}^{u*}+{\displaystyle \sum _{k=-\infty }^{\infty }{\epsilon }_k^{u*}\exp \left(ik{\Omega }_0^{*}{t}^{*}\right)}\\ R^{*}\left(t^{*}\right)={\epsilon }_{ss}^{l*}+{\displaystyle \sum _{m=-\infty }^{\infty }{\epsilon }_m^{l*}\exp \left(im{\mathrm{\Theta }}_0^{*}{t}^{*}\right)}\end{array}$$

(17)

where:

$${\left({\epsilon }_{ss}^u\right)}^{*}={\epsilon }_{ss}^u\frac{U_w}{U_{ch}}\hspace{1em}\hspace{1em}{\left({\epsilon }_k^u\right)}^{*}={\epsilon }_k^u\frac{U_w}{U_{ch}}\hspace{1em}\hspace{1em}{P}_{ss}^{*}=\frac{P_{ss}{h}^2}{2\mu \hspace{0.17em}{U}_{ch}}\hspace{1em}\hspace{1em}\left|P_n^{*}\right|=\frac{\left|P_n\right|h^2}{2\mu \hspace{0.17em}{U}_{ch}}$$

The asterisk here refers to dimensionless and not the complex conjugate values. In order to avoid cluttering the symbols and equations we are going to intentionally omit asterisk superscript for dimensionless quantities in subsequent derivations.

2.4. Boundary Conditions and Stability of Pulsatile Flows

The question of Dirichlet’s no-slip BCs for velocity in micro-fluidic flows of similar geometry has already been thoroughly discussed in [18]. Additionally, the general mixed BCs often referred to as Robin BCs for fluid velocity can be expressed as in the case of partial slip:

$$\Delta u_{wall}=u\left(0^{+},t\right)-u_{wall}=L_s\hspace{0.17em}\hspace{0.17em}{\left(\frac{\partial \hspace{0.17em}u}{\partial y}\right)}_{y\to 0^{+}}$$

(18)

where L_s is Navier–Maxwell slip length. Here, we have assumed full no-slip condition in both forward and return wall motion and the Robin BC boils down to Dirichlet BC.

The classical linear stability theory [2,4,5,6,10] predicts absolute stability for steady planar CT and PO flows in rigid tubes, i.e., entirely contrary to experimental evidence. On the other hand, plane PO flow has critical Reynolds number of about 5750 obtained by solving the Orr-Sommerfeld equation [2,4,5,6,10]. Experimentally, the first Tollmien–Schlichting (TS) wave appears, and the onset of flow transition occurs at Reynolds numbers as low as 360. The stability of pulsatile flows in cylindrical and plane geometries is of equal importance and is gaining more importance for microfluidic applications. The science of the flow stability is vast and cannot be treated here in any detail. It suffices to say that the maximum Reynolds number occurring in our analysis stays well below the critical and the assumption of an all-time laminar flow thus seems appropriate.

3. Methods and Materials

To obtain the closed-form analytical solution of this linear non-homogeneous IBVP problem we use the powerful method of eigenfunction expansion (EEM) [39,40,41,42,43]. The entire linear non-homogeneous IBVP problem can be treated directly by EEM. However, in order to arrive at a general solution in a physically more transparent way we will decompose the velocity field into a steady-state, transient and quasi-steady-state (QSS). This is permissible as the governing PDE is linear and the general principle of linear superposition applies. Keeping in mind that all variables are dimensionless, we can write for the velocity field:

$$u(y,t)=u_{ss}\left(y\right)+v(y,t)+w(y,t)$$

(19)

The first component on the RHS is the steady CP flow that may also contain DC PG components from the respective Fourier expansions. The 2nd term represents transient response to ICs and will vanish after the finite amount of time (rapidly in most cases). The 3rd velocity component represents the QSS velocity field (ST-WO component) as a response to arbitrary non-harmonic boundary motion and unsteady PG. However, the third component may also contain some transient responses due to phase information contained in the Fourier expansions.

To solve resulting IBVPs by EEM we invoke the formidable and well-developed mathematical apparatus of self-adjoint linear operators in infinite-dimensional separable Hilbert function space [44,45,46,47,48] (or L²-metric space in particular). It is the same theoretical basis that lies at the mathematical foundation of the quantum mechanics. The eigenfunctions of the self-adjoint (Hermitian) linear differential operator in Hilbert space form a closed and complete basis set and are used in generalized Fourier series expansions.

3.1. Solution of the Steady-State BVP

The first BVP defines steady CP flow and is described in a dimensionless form as:

$$\begin{array}{l}0=\frac{d^2{u}_{ss}}{d{y}^2}+2P_{ss}\\ \mathrm{BCs}:v\left(+1\right)=G_{ss}={\epsilon }_{ss}^u\hspace{1em}\hspace{1em}v\left(-1\right)=R_{ss}={\epsilon }_{ss}^l\hspace{1em}\hspace{1em}{P}_{ss}=\mathrm{const}.\end{array}$$

(20)

The solutions for the dimensionless velocity field, flow rate, velocity gradient, tangential stress, vorticity, and dissipation ($P_{ss}\le \ge 0,\hspace{0.17em}-1\le y\le +1$) are:

$$\begin{array}{l}{u}_{ss}\left(y\right)=\left[P_{ss}\left(1-y^2\right)+\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{2}+\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)y\right]\hspace{1em}\hspace{1em}\hspace{1em}{Q}_{ss}=\left({\epsilon }_{ss}^u+{\epsilon }_{ss}^l\right)+\frac{4P_{ss}}{3}=2\cdot {\overline{u}}_{ss}\\ \frac{\partial u\left(y,t\right)}{\partial y}=\frac{d{u}_{ss}\left(y\right)}{dy}=\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)-2\cdot P_{ss}\cdot y\hspace{1em}\hspace{1em}\hspace{1em}{\tau }_{zy}=\hspace{0.17em}\frac{d{u}_{ss}\left(y\right)}{dy}=\hspace{0.17em}\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)-2\cdot P_{ss}\cdot y\\ {\omega }_z=-\frac{d{u}_{ss}\left(y\right)}{dy}=2\cdot P_{ss}\cdot y-\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)\hspace{1em}\hspace{1em}\hspace{1em}\varphi =\hspace{0.17em}{\left(\hspace{0.17em}\frac{d{u}_{ss}}{dy}\right)}^2=\hspace{0.17em}\left[{\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)}^2-2\cdot \left({\epsilon }_{ss}^u-{\epsilon }_{ss}^l\right)\hspace{0.17em}{P}_{ss}\cdot y+4\cdot P_{ss}^2\cdot y^2\right]\end{array}$$

(21)

In general CP flow, the uniform wall velocity (of a faster moving wall) is used to obtain the dimensionless form and the non-dimensional PG can assume any value ($P_{ss}\le \ge 0$). Clearly, a plethora of uniform wall movements in combination with a choice of static (proverse or adverse) PGs are possible.

If bounding infinite walls are immovable, we obtain the classical plane PO flow, $u_{ss}\left(y\right)=P_{ss}\left(1-y^2\right)$, with, $P_{ss}=1$. In the plane PO flow, the wall velocity is zero and cannot be used for scaling purposes. Therefore, the characteristic velocity is the maximum velocity in the center and dimensionless PG must be $P_{ss}=1$, in accordance with Equations (13) and (14).

3.2. Solution of the Homogeneous IBVP

An infinite-dimensional Hilbert space is an inner-product separable normed metric space with definite algebra and topology. Completeness of eigenfunctions in L² or Hilbert space is guaranteed by, often called, Parseval’s equality [39,40,41,42,43] which forms the equivalence part of the Bessel’s inequality and represents energy conservation principle (completeness of eigenfunctions). The L² or mean-square convergence is automatically assured for complete basis sets only requiring that the function is square-integrable in the Lebesgue integral (absolute) sense or is Lebesgue-measurable [44,45,46,47,48,49,50]:

$${\Vert f\Vert }^2={\displaystyle \underset{a}{\overset{b}{\int }}{\left|f\left(x\right)\right|}^{2}dx\hspace{0.17em}}\hspace{0.17em}<\infty$$

(22)

The mathematical apparatus of functional analysis and L² space used here is based on the celebrated Riesz–Fisher theorem [44,45,46,47] which establishes the isomorphism of vector space l² and function space L² as Hilbert spaces. Assuming a complete basis set of eigenfunctions, a Lebesgue square-integrable function (vast majority also Riemann-integrable) can be represented as a linear combination of orthonormal eigenfunctions:

$$f={\displaystyle \sum _{n=1}^{\infty }{C}_n{Y}_n}\hspace{1em}\hspace{1em}\hspace{1em}f={\displaystyle \sum _{n=1}^{\infty }\langle f,Y_n\rangle Y_n}$$

where the generalized Fourier coefficients $C_n=\langle f,Y_n\rangle$ follow from orthogonality:

$$\langle f,Y_m\rangle ={\displaystyle \sum _{n=1}^{\infty }{C}_n\langle Y_n,Y_m\rangle }={\displaystyle \sum _{n=1}^{\infty }{C}_n{\delta }_{nm}}=C_m\hspace{1em}\hspace{1em}\langle Y_n,Y_n\rangle ={\Vert Y_n\Vert }^2=1$$

(23)

The second problem is a simple homogeneous IBVP with modified ICs originating due to steady-state CP flow:

$$\begin{array}{l}{v}_t=v_{yy}\\ \mathrm{BCs}:v\left(+1,t\right)=v\left(-1,t\right)=0\hspace{1em}\\ \mathrm{IC}:v\left(y,0\right)=f\left(y\right)-u_{ss}\left(y\right)=F\left(y\right)\end{array}$$

(24)

This is a parabolic PDE problem that can be solved using the SOV [31,32,33,34] method (sometimes also called Fourier method). The solution of the simple Sturm–Liouville ($-Y_n″={\lambda }_n{Y}_n$) BVP is:

$$v\left(y,t\right)={\displaystyle \sum _{n=1}^{\infty }{A}_n\cdot T_n\left(t\right)\cdot Y_n\left(y\right)}$$

(25)

The cosine eigenfunction would be a solution of the BVP if the flow is symmetric (BCs are identical or periodic). The eigenvalues (zeros) of the cosine eigenfunctions are:

$${\lambda }_n=k_n^2={\left[\left(2n-1\right)\frac{\pi }{2}\right]}^2\hspace{1em}\hspace{1em}n=1,2,3,\dots$$

(26)

However, when problems are asymmetric, we need both cosine and sine eigenfunctions appearing alternatively and with changing signs. The general asymmetric solution can be thus represented by the following eigenvalues and eigenfunctions:

$${\lambda }_n=k_n^2={\left(n\pi \right)}^2\hspace{1em}\hspace{1em}{Y}_n\left(y\right)=\sin \left[k_n\left(\frac{y+1}{2}\right)\right]\hspace{1em}\hspace{1em}n=1,2,3,\dots$$

(27)

This eigenfunction is basically a sum of cosine and sine eigenfunctions appearing alternatively in infinite sums. All eigenvalues here must be non-negative real numbers as required by the regular Sturm–Liouville (S-L) theory [39,40,41,42,43,51,52,53,54]. Classical Fourier series are only one of the infinitely many possible complete eigenfunction basis sets.

We mostly treat symmetric flows here with cosine eigenfunctions, describing a complete basis set of normalized and orthogonal (orthonormal) eigenfunctions:

$$\langle Y_n,Y_m\rangle ={\displaystyle \underset{-1}{\overset{+1}{\int }}{Y}_n\left(y\right)\cdot Y_m\left(y\right)\hspace{0.17em}dy}={\displaystyle \underset{-1}{\overset{+1}{\int }}\cos \left(k_{n}y\right)\cdot \cos \left(k_{m}y\right)}\hspace{0.17em}dy={\delta }_{nm}\hspace{1em}\hspace{1em}\langle Y_n,Y_n\rangle ={\Vert Y_n\Vert }^2=1$$

We have not yet stated ICs explicitly. For example, a PO parabolic velocity profile in a slit could be used for studying stopping (shutdown) flow dynamics. For simplicity and to obtain fast closed-form solutions, we assume uniform velocity profiles only (say, zero or one):

$$\mathrm{IC}:\hspace{0.17em}{\epsilon }_i^0={\epsilon }_i\frac{U_0}{U_{ch}}\hspace{1em}\hspace{1em}0\le {\epsilon }_i^0\le 1$$

The Fourier coefficients of the 2nd IBVP are now:

$$A_n=\frac{\langle F,Y_n\rangle }{{\Vert Y_n\Vert }^2}=\langle F,Y_n\rangle ={\left(-1\right)}^{n+1}\frac{2}{k_n}\left[{\epsilon }_i^0-\frac{2P_{ss}}{k_n^2}-\left(\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{2}\right)\right]\hspace{1em}\hspace{1em}{k}_n=\frac{\left(2n-1\right)\pi }{2}\hspace{1em}n=1,2,3,\dots$$

(28)

The complete solution of the symmetric (decaying) transient flow is a linear superposition involving several ICs:

$$v\left(y,t\right)={\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\frac{2}{k_n}\left[{\epsilon }_i^0-\frac{2P_{ss}}{k_n^2}-\left(\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{2}\right)\right]\cdot e^{-k_n^{2}t}\cdot \cos \left(k_{n}y\right)}$$

(29)

The first term in the above Fourier expansion describes the effect of the initial velocity distribution. The second term gives transients due to sudden application of the constant PG, while the third component is the response to the sudden instant finite uniform movement (ST1 problem) of the upper and/or lower infinite plate at time zero. Negative exponential factor ensures decaying flow with higher flow harmonics diminishing more rapidly. Starting and stopping flow transients can be deduced from the general solution here.

3.3. Solution of the Non-Homogeneous IBVP for Symmetric Flows

The third and final IBVP problem which has both inhomogeneous PDE and BCs.

$$\begin{array}{l}{w}_t=w_{yy}+q\left(y,t\right)\\ \mathrm{BCs}:w\left(+1,t\right)=g\left(t\right)\hspace{1em}\hspace{1em}w\left(-1,t\right)=r\left(t\right)\\ \mathrm{IC}:w\left(y,0\right)=0\\ q\left(y,t\right)={\displaystyle \sum _{m=-\infty }^{\infty }2\left|P_m\right|\exp \left[i\left(m{\omega }_{0}t+{\phi }_m\right)\right]}\\ g\left(t\right)={\displaystyle \sum _{k=-\infty }^{\infty }\left|{\epsilon }_k^u\right|\exp \left[i\left(k{\Omega }_0^{u}t+{\varphi }_k\right)\right]}\hspace{1em}\hspace{1em}\hspace{1em}r\left(t\right)={\displaystyle \sum _{j=-\infty }^{\infty }\left|{\epsilon }_j^l\right|\exp \left[i\left(j{\Omega }_0^{l}t+{\psi }_j\right)\right]}\end{array}$$

(30)

Both arbitrary, yet periodic, BCs and the pressure source term are given in terms of complex Fourier coefficients. We can always find a suitable (shifting) function that will satisfy the non-homogeneous BCs [19,31,32,33,34,35]. In this way we obtain homogeneous BCs, but the new PDE often ends up with additional source terms. While this technique is perfectly good and often results in solutions containing uniformly convergent series in the entire domain, it is time consuming and problem specific. More general, elegant, and faster solutions are obtained by the straightforward use of EEM although we may get discontinuity at the boundaries thus invalidating uniform convergence in the entire domain. However, mean-square convergence is automatically assured due to the completeness of eigenfunctions. Other analytical methods, such as Laplace transform, Finite Fourier Transform, Green’s functions, etc., could be used to solve bounded non-homogenous problem [20,39,40,41,42,43].

The essence of the EEM is to tailor the inhomogeneous solution in terms of the orthonormal eigenfunctions φ_n(y) of the corresponding (associated) homogeneous BVP (CHBVP):

$$\begin{array}{l}w\left(y,t\right)={\displaystyle \sum _{n=1}^{\infty }{T}_n\left(t\right)\cdot {\phi }_n\left(y\right)}={\displaystyle \sum _{n=1}^{\infty }{T}_n\left(t\right)\cdot \cos \left(k_{n}y\right)}\\ q\left(y,t\right)={\displaystyle \sum _{n=1}^{\infty }{q}_n\left(t\right)\cdot {\phi }_n\left(y\right)}={\displaystyle \sum _{n=1}^{\infty }{q}_n\left(t\right)\cdot \cos \left(k_{n}y\right)}\\ w\left(y,0\right)={\displaystyle \sum _{n=1}^{\infty }{T}_n\left(0\right)\cdot {\phi }_n\left(y\right)}={\displaystyle \sum _{n=1}^{\infty }{T}_n\left(0\right)\cdot \cos \left(k_{n}y\right)}\end{array}$$

(31)

By utilizing the orthogonal property of the eigenfunctions and the norm, we obtain:

$$\begin{array}{l}{T}_n\left(t\right)=\langle w,Y_n\rangle ={\displaystyle \underset{-1}{\overset{+1}{\int }}w\left(y,t\right)\cdot \cos \left(k_{n}y\right)\hspace{0.17em}dy\hspace{1em}\hspace{1em}n=1,2,3\dots }\\ q_n\left(t\right)=\langle q,Y_n\rangle ={\displaystyle \underset{-1}{\overset{+1}{\int }}q\left(y,t\right)\cdot \cos \left(k_{n}y\right)\hspace{0.17em}dy}\\ T_n\left(0\right)=\langle 0,Y_n\rangle ={\displaystyle \underset{-1}{\overset{+1}{\int }}w\left(y,0\right)\cdot \cos \left(k_{n}y\right)\hspace{0.17em}dy}=0\end{array}$$

(32)

Using these relationships, the decomposed oscillating PG source term is:

$$q_n\left(t\right)={\left(-1\right)}^{n+1}\frac{4}{k_n}\left\{{\displaystyle \sum _{m=-\infty }^{\infty }\left|P_m\right|\exp \left[i\left(m{\omega }_{0}t+{\varphi }_m\right)\right]}\right\}$$

(33)

One may think that all IC transients were taken care in the homogeneous IBVP (2nd problem), but that is not the case. Although the ICs are explicitly zero here, the initial information is also hidden in phasors. Therefore, we are dealing with the QSS response caused by periodic forcing functions and few transients caused by non-zero phases. Any transient due to sudden application of the oscillating PG will be taken care of in the solution of this non-homogeneous IBVP. With compatibility condition satisfied [42], the resulting Fourier series is uniformly convergent in the entire domain including boundaries.

Since the term-by-term differentiation is prohibited for the Fourier series that does not converge uniformly, when utilizing Green’s 2nd identity [39,40,41,42,43], one obtains:

$$\frac{d\hspace{0.17em}{T}_n\left(t\right)}{d\hspace{0.17em}t}+k_n^2{T}_n\left(t\right)={\left(-1\right)}^{n+1}{k}_n\left[g\left(t\right)+r\left(t\right)\right]+q_n\left(t\right)=Q_n(t)\hspace{1em}\hspace{1em}{T}_n\left(0\right)=0\hspace{1em}\hspace{1em}n=1,2,3,\dots$$

(34)

Implementing the general solution of the inhomogeneous linear first-order ODE, we obtain:

$$\begin{array}{l}{T}_n\left(t\right)={\left(-1\right)}^{n+1}{k}_n{\displaystyle \sum _{j=-\infty }^{\infty }\left\{\frac{\left|{\epsilon }_j^u\right|\exp \left[i\left({\Omega }_j^{u}t+{\varphi }_j\right)\right]}{\left(k_n^2+i{\Omega }_j^u\right)}-\frac{\left|{\epsilon }_j^u\right|\exp \left(i{\varphi }_j\right)}{\left(k_n^2+i{\Omega }_j^u\right)}\exp \left(-k_n^{2}t\right)\right\}}\\ +{\left(-1\right)}^{n+1}{k}_n{\displaystyle \sum _{p=-\infty }^{\infty }\left\{\frac{\left|{\epsilon }_p^l\right|\exp \left[i\left({\Omega }_p^{l}t+{\psi }_p\right)\right]}{\left(k_n^2+i{\Omega }_p^l\right)}-\frac{\left|{\epsilon }_p^l\right|\exp \left(i{\psi }_p\right)}{\left(k_n^2+i{\Omega }_p^l\right)}\exp \left(-k_n^{2}t\right)\right\}}\\ +{\left(-1\right)}^{n+1}\frac{4}{k_n}{\displaystyle \sum _{m=-\infty }^{\infty }\left\{\frac{\left|P_m\right|\exp \left[i\left({\omega }_{m}t+{\phi }_m\right)\right]}{\left(k_n^2+i{\omega }_m\right)}-\frac{\left|P_m\right|\exp \left(i{\phi }_m\right)}{\left(k_n^2+i{\omega }_m\right)}\exp \left(-k_n^{2}t\right)\right\}}\end{array}$$

(35)

The solution of the 3rd non-homogenous IBVP problem is thus:

$$\begin{array}{l}w\left(y,t\right)={\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}{k}_n\cos \left(k_{n}y\right)\left\{{\displaystyle \sum _{j=-\infty }^{\infty }\frac{\left|{\epsilon }_j^u\right|\exp \left[i\left({\Omega }_j^{u}t+{\varphi }_j\right)\right]}{\left(k_n^2+i{\Omega }_j^u\right)}}\right\}-{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}{k}_n\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left\{{\displaystyle \sum _{j=-\infty }^{\infty }\frac{\left|{\epsilon }_j^u\right|\exp \left(i{\varphi }_j\right)}{\left(k_n^2+i{\Omega }_j^u\right)}}\right\}}}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}{k}_n\cos \left(k_{n}y\right)\left\{{\displaystyle \sum _{p=-\infty }^{\infty }\frac{\left|{\epsilon }_p^l\right|\exp \left[i\left({\Omega }_p^{l}t+{\psi }_p\right)\right]}{\left(k_n^2+i{\Omega }_p^l\right)}}\right\}}-{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}{k}_n\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left\{{\displaystyle \sum _{p=-\infty }^{\infty }\frac{\left|{\epsilon }_p^l\right|\exp \left(i{\psi }_p\right)}{\left(k_n^2+i{\Omega }_p^l\right)}}\right\}}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\frac{\cos \left(k_{n}y\right)}{k_n}\left\{{\displaystyle \sum _{m=-\infty }^{\infty }\frac{4\left|P_m\right|\exp \left[i\left({\omega }_{m}t+{\phi }_m\right)\right]}{\left(k_n^2+i{\omega }_m\right)}}\right\}}-{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n\frac{\cos \left(k_{n}y\right)e^{-k_n^{2}t}}{k_n}\left\{{\displaystyle \sum _{m=-\infty }^{\infty }\frac{4\left|P_m\right|\exp \left(i{\phi }_m\right)}{\left(k_n^2+i{\omega }_m\right)}}\right\}}\end{array}$$

(36)

After tedious complex algebra reductions, one obtains:

$$\begin{array}{l}w\left(y,t\right)=-{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\cos \left({\varphi }_j-{\gamma }_{j,n}^u\right)}\right]\hspace{0.17em}}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[{\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\psi }_p-{\gamma }_{p,n}^l\right)}\right]}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[\frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\hspace{0.17em}\cos \left({\phi }_m-{\beta }_{m,n}\right)}\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)\left[{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\hspace{0.17em}\cos \left({\Omega }_j^{u}t+{\varphi }_j-{\gamma }_{j,n}^u\right)}\right]\hspace{0.17em}}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)\left[{\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\Omega }_p^{l}t+{\psi }_p-{\gamma }_{p,n}^l\right)}\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)\left[\frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\cos \left({\omega }_{m}t+{\phi }_m-{\beta }_{m,n}\right)}\right]}\end{array}$$

(37)

3.4. Global Solution for Symmetrically Driven Flows

A global solution of the symmetric IBVP is obtained by linearly superposing all three velocity components as given by Equation (19). After grouping all terms, the final dimensionless solution yields:

$$\begin{array}{l}u\left(y,t\right)=\left[P_{ss}\left(1-y^2\right)+\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{2}+\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)y\right]+{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\hspace{0.17em}\left[\frac{2{\epsilon }_i^0}{k_n}-\frac{4P_{ss}}{k_n^3}-\left(\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{k_n}\right)\right]}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\cos \left({\varphi }_j-{\gamma }_{j,n}^u\right)}\right]\hspace{0.17em}}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[{\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\psi }_p-{\gamma }_{p,n}^l\right)}\right]}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)e^{-k_n^{2}t}\left[\frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\hspace{0.17em}\cos \left({\phi }_m-{\beta }_{m,n}\right)}\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\hspace{0.17em}\cos \left(k_{n}y\right)\left[{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\hspace{0.17em}\cos \left({\Omega }_j^{u}t+{\varphi }_j-{\gamma }_{j,n}^u\right)}\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)\left[{\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\Omega }_p^{l}t+{\psi }_p-{\gamma }_{p,n}^l\right)}\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\cos \left(k_{n}y\right)\left[\frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\cos \left({\omega }_{m}t+{\phi }_m-{\beta }_{m,n}\right)}\right]}\end{array}$$

(38)

Phase and amplitudes (damping factors) of the oscillation response are:

$$\begin{array}{l}{\gamma }_{j,n}^u={\tan }^{-1}\left(\frac{{\Omega }_j^u}{k_n^2}\right)={\tan }^{-1}\left[\frac{{\left({\alpha }_j^u\right)}^2}{k_n^2}\right]={\tan }^{-1}{\delta }_{j,n}^u\hspace{1em}\hspace{1em}\hspace{1em}{C}_{j,n}^u=\frac{2}{k_n\sqrt{1+{\left({\delta }_{j,n}^u\right)}^2}}\\ {\gamma }_{p,n}^l={\tan }^{-1}\left(\frac{{\Omega }_p^l}{k_n^2}\right)={\tan }^{-1}\left[\frac{{\left({\alpha }_p^l\right)}^2}{k_n^2}\right]={\tan }^{-1}{\delta }_{p,n}^l\hspace{1em}\hspace{1em}\hspace{1em}{C}_{p,n}^l=\frac{2}{k_n\sqrt{1+{\left({\delta }_{p,n}^l\right)}^2}}\\ {\beta }_{m,n}={\tan }^{-1}\left(\frac{{\omega }_m}{k_n^2}\right)={\tan }^{-1}\left(\frac{{\alpha }_m^2}{k_n^2}\right)={\tan }^{-1}{\delta }_{m,n}^{pr}\hspace{1em}\hspace{1em}\hspace{1em}{D}_{m,n}=\frac{4}{k_n^3\sqrt{1+{\left({\delta }_{m,n}^{pr}\right)}^2}}\end{array}$$

(39)

The damping factors are positive-definite. Higher interface-motion harmonics will generate correspondingly higher Womersley (WO) numbers, ${\alpha }_n^u={\alpha }_0^u\sqrt{n}$, thus penetrating ever shallower into the fluid mass, but significantly affecting its phase lags. The solution consists of two traveling waves moving in the opposite directions, i.e., coming from each wall. Different regions of fluid are experiencing varying effects of resistance and inertia. Daidzic and Hossain [18] and Hossain and Daidzic [19] have demonstrated similar behavior for Newtonian and some shear-thinning non-Newtonian fluids.

An arbitrary periodic function $f\left(t\right)$ can be represented by converging the Fourier series. At least an L² or mean-square convergence is guaranteed for Lebesgue square-integrable functions:

$$\underset{N\to \infty }{\lim }\hspace{0.17em}{\displaystyle \underset{0}{\overset{T}{\int }}{\left|f\left(t\right)-S_N\left(t\right)\hspace{0.17em}\right|}^{2}dt}\hspace{0.17em}=0\hspace{1em}\hspace{1em}{S}_N\left(t\right)={\displaystyle \sum _{n=1}^N{f}_n\left(t\right)}\hspace{1em}\hspace{1em}\forall t\in \left[0,T\right]$$

Mean-square convergence does not imply pointwise or uniform convergence [39,40,41,42,43,51,52,53,54]. At the finite points of jump discontinuities of piecewise continuous functions, the Fourier series will converge to the average value of the left and the right limits with Gibbs phenomena and characteristic “ringing” appearing [39,40,41,42,43,51,52,53,54]. Accordingly, the convergence cannot possibly remain uniform. Once we have chosen oscillation waveforms, their Fourier coefficients can be directly substituted into the global solution (Equation (38)). The question of convergence is not the academic one. The nature of convergence strongly affects analytical solutions and computations. From the velocity distribution in the complex form, we can easily derive volume flow rate (VFR):

$$\begin{array}{l}Q\left(t\right)={\displaystyle \underset{-1}{\overset{+1}{\int }}u\left(y,t\right)dy}=\left[\left({\epsilon }_{ss}^u+{\epsilon }_{ss}^l\right)+\frac{4P_{ss}}{3}\right]+{\displaystyle \sum _{n=1}^{\infty }\frac{2\hspace{0.17em}{e}^{-k_n^{2}t}}{k_n}\hspace{0.17em}\left[\frac{2{\epsilon }_i^0}{k_n}-\frac{4P_{ss}}{k_n^3}-\left(\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{k_n}\right)\right]}\\ -{\displaystyle \sum _{n=1}^{\infty }\frac{2\hspace{0.17em}{e}^{-k_n^{2}t}}{k_n}\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\cos \left({\varphi }_j-{\gamma }_{j,n}^u\right)}+\\ {\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\psi }_p-{\gamma }_{p,n}^l\right)}+\\ \frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\hspace{0.17em}\cos \left({\phi }_m-{\beta }_{m,n}\right)}\end{array}\right\}}+{\displaystyle \sum _{n=1}^{\infty }\frac{2}{k_n}\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\hspace{0.17em}\cos \left({\Omega }_j^{u}t+{\varphi }_j-{\gamma }_{j,n}^u\right)}+\\ {\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\Omega }_p^{l}t+{\psi }_p-{\gamma }_{p,n}^l\right)}+\\ \frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\cos \left({\omega }_{m}t+{\phi }_m-{\beta }_{m,n}\right)}\end{array}\right\}}\end{array}$$

(40)

The velocity gradient of the axial velocity is:

$$\begin{array}{l}\frac{\partial u\left(y,t\right)}{\partial y}=\left[-2P_{ss}\cdot y+\left(\frac{{\epsilon }_{ss}^u-{\epsilon }_{ss}^l}{2}\right)\right]-{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}{k}_n{e}^{-k_n^{2}t}\sin \left(k_{n}y\right)\left[\frac{2{\epsilon }_i^0}{k_n}-\frac{4P_{ss}}{k_n^3}-\left(\frac{{\epsilon }_{ss}^u+{\epsilon }_{ss}^l}{k_n}\right)\right]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n{k}_n{e}^{-k_n^{2}t}\sin \left(k_{n}y\right)\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\cos \left({\varphi }_j-{\gamma }_{j,n}^u\right)}+\\ {\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\psi }_p-{\gamma }_{p,n}^l\right)}+\\ \frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\hspace{0.17em}\cos \left({\phi }_m-{\beta }_{m,n}\right)}\end{array}\right\}}\hspace{0.17em}\\ -{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n{k}_n\sin \left(k_{n}y\right)\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\left|{\epsilon }_j^u\right|C_{j,n}^u\hspace{0.17em}\cos \left({\Omega }_j^{u}t+{\varphi }_j-{\gamma }_{j,n}^u\right)}+\\ {\displaystyle \sum _{p=1}^{\infty }\left|{\epsilon }_p^l\right|C_{p,n}^l\cos \left({\Omega }_p^{l}t+{\psi }_p-{\gamma }_{p,n}^l\right)}+\\ \frac{4P_0}{k_n^3}+{\displaystyle \sum _{m=1}^{\infty }2\left|P_m\right|D_{m,n}\cos \left({\omega }_{m}t+{\phi }_m-{\beta }_{m,n}\right)}\end{array}\right\}}\end{array}$$

(41)

A word of caution is required. Term-by-term differentiation of a Fourier series is a delicate matter effectively multiplying a series by a factor (n), slowing its convergence, and is only allowed if the series and differentiated series converges uniformly [51,52,53,54,55,56]. Term-by-term integration is usually a less critical issue as a multiplication factor (1/n) is introduced accelerating the convergence [51,52,53,54,55,56]. The classical Fourier series can always be integrated term-by-term [51,52,53,54,55,56].

3.5. Numerical Method

To verify theoretical solutions, we employ numerical experiments using in-house-developed CFD programs. The finite-volume method (FVM) to discretize N-S equations [57,58,59,60] and to numerically calculate the flow field was employed. FVM is in many respects similar to the FD method, but spatial-temporal discretization is obtained while assuring conservation of flow properties in every small finite volume. To verify the accuracy, stability, consistency, and efficiency of numerical models, we test our CFD against several well-known analytic solutions [18,19]. Despite known shortcomings, the conditionally stable explicit FTCS discretization works very well for the 1D flow case. The problem with the Crank–Nicolson (C-N) or the fully implicit methods is that we must solve the system of linear algebraic equations at each time step. Fortunately, for 1D flow an efficient tri-diagonal matrix solver utilizing Thomas’ (TDMA) Gauss-elimination algorithm is used [18,19,57,58,59,60]. For 2D and 3D flows, we have to employ different methodologies. While it is often stated that the implicit schemes are unconditionally stable that may be useless feature when the accuracy suffers due to excessively large time steps.

The solution of the linear parabolic PDE given in Equation (5) is quite simple and can be most rapidly obtained by marching in time while complying with given IC and BCs. An added feature of the FVM is that it automatically satisfies conservation laws as each balance differential equation is integrated over a small control (finite) volume. It also accounts for all source (volumetric) and surface interactions. The numerical methodology implemented here was already sufficiently described earlier [18,19] and we only repeat here the most essential parts.

The “old” values of velocity at grid points can be denoted by $u_P^m\hspace{0.17em},\hspace{0.17em}{u}_E^m\hspace{0.17em},\mathrm{and},\hspace{0.17em}{u}_W^m$, and the updated values at the “new” time $t+\Delta t$, can be denoted by $u_P^{m+1}\hspace{0.17em},\hspace{0.17em}{u}_E^{m+1}\hspace{0.17em},\mathrm{and},\hspace{0.17em}{u}_W^{m+1}$. If we designate:

$$a_P^0\hspace{0.17em}=a_j^0\hspace{0.17em}=\hspace{0.17em}\rho \hspace{0.17em}\frac{\Delta y}{\Delta t}\hspace{1em}\hspace{1em}{a}_E\hspace{0.17em}=a_{j+1}\hspace{0.17em}=\hspace{0.17em}\frac{{\mu }_e}{{\left(\delta y\right)}_e}\hspace{1em}\hspace{1em}{a}_W\hspace{0.17em}=a_{j-1}\hspace{0.17em}=\hspace{0.17em}\frac{{\mu }_w}{{\left(\delta y\right)}_w}$$

(42)

The spatial-temporal discretized equation yields:

$$\begin{array}{l}{a}_j\hspace{0.17em}{u}_j^{m+1}=a_{j+1}\left\{f\cdot u_{j+1}^{m+1}+\left(1-f\right)\cdot u_{j+1}^m\right\}+a_{j-1}\left\{f\cdot u_{j-1}^{m+1}+\left(1-f\right)\cdot u_{j-1}^m\right\}\hspace{0.17em}\\ \hspace{1em}\hspace{1em}\hspace{1em}+\hspace{0.17em}\left\{a_j^0-a_{j+1}\left(1-f\right)-a_{j-1}\left(1-f\right)\right\}\cdot u_j^m\end{array}$$

(43)

where

$$a_j=a_j^0+f\cdot a_{j+1}+f\cdot a_{j-1}.$$

Here, f stands for the time-approximation weighting factor with the range [0, 1]. The value of f will determine whether the discretization scheme is explicit, fully implicit or C-N. The unconditionally stable popular implicit C-N (CTCS) method ($f\hspace{0.17em}=\hspace{0.17em}1/2$) is second-order accurate, both, in time and space:

$$\left(a_j^0+\frac{a_{j+1}+a_{j-1}}{2}\right)\cdot u_j^{m+1}\hspace{0.17em}=\hspace{0.17em}\frac{a_{j+1}}{2}\left(u_{j+1}^{m+1}+u_{j+1}^m\right)+\frac{a_{j-1}}{2}\left(u_{j-1}^{m+1}+u_{j-1}^m\right)+\hspace{0.17em}\left(a_j^0\hspace{0.17em}-\frac{a_{j+1}+\hspace{0.17em}{a}_{j-1}}{2}\right)\cdot u_j^m$$

(44)

A simple explicit (FTCS) discretization with the time-dependent PG as a source function yields:

$$u_j^{m+1}\hspace{0.17em}=\hspace{0.17em}\left(1-2\lambda \right)u_j^m+\lambda \left(u_{j+1}^m+u_{j-1}^m\right)+P^m\Delta \hspace{0.17em}t\hspace{1em}\hspace{1em}\hspace{1em}\lambda =\frac{\Delta \hspace{0.17em}t}{{\left(\Delta \hspace{0.17em}y\right)}^2}$$

(45)

Appropriate IC and BCs must be applied for the numerical computations. By choosing the proper stability parameter λ(=1/4) we can, in fact, obtain second-order accuracy in time.

3.6. Modeling of Pressure Gradients

The most comprehensive planar CPSW problem treated here is the one that includes both the periodic non-harmonic oscillations of two boundary wall plates in concurrence with the applied oscillatory periodic PGs. We assume that both plates move simultaneously and in phase and thus the problem is symmetric. The oscillatory PG is superposed to a steady-state PG when applicable. Clearly an infinite number of combinations of wall and PG dynamics are possible.

We chose a rectified-sine (absolute-sine) function to represent PG oscillatory dynamics. Another interesting choice of pressure oscillations would be in emulating heartbeat waveforms for pulsatile arterial hemodynamic applications. Digitally sampled waveforms can be utilized by employing harmonic analysis and calculating the Fourier coefficients of an arbitrary discrete excitation. A full-wave rectified-sine (absolute-sine) function delivering oscillating PGs on top of the steady-state PG for a horizontal slit is illustrated in Figure 2, and is described by:

$$q\left(t\right)=P_{ss}+f\left(t\right)=P_{ss}+P_{\max }\cdot \left|\sin \left(\frac{2\hspace{0.17em}\pi \hspace{0.17em}t}{T}\right)\right|\hspace{1em}\hspace{1em}{P}_{ss}\le \ge 0$$

(46)

This even real-valued function can be expressed in terms of Fourier cosine series. By using Euler formulas [39,40,41,42,43,51,52,53,54,55,56] for Fourier coefficients, one obtains:

$$f\left(t\right)=\frac{2\hspace{0.17em}{P}_{\max }}{\pi }-\frac{4\hspace{0.17em}{P}_{\max }}{\pi }{\displaystyle \sum _{n=2,4,6,\dots }^{\infty }\frac{\cos \left(\frac{2\hspace{0.17em}n\pi }{T}t\right)}{{\left(n^2-1\right)}^{}}}=\frac{2\hspace{0.17em}{P}_{\max }}{\pi }-\frac{4\hspace{0.17em}{P}_{\max }}{\pi }{\displaystyle \sum _{m=1}^{\infty }\frac{\cos \left(2\hspace{0.17em}m\omega \hspace{0.17em}t\right)}{\left(2\hspace{0.17em}m-1\right)\left(2\hspace{0.17em}m+1\right)}}$$

(47)

It is easy to demonstrate, using the Weierstrass convergence M-test [40,41,42], that this Fourier series converges uniformly and accordingly pointwise as well. Function in Equation (46) is square-integrable and its Fourier series thus converges in the mean (L² sense). This was already expected from the uniform convergence. Utilizing Parseval’s identity in the real form yields:

$$\frac{1}{T}{\displaystyle \underset{-T/2}{\overset{T/2}{\int }}{\left|f\left(t\right)\right|}^{2}dt}=\frac{A_0^2}{4}+\frac{1}{2}{\displaystyle \sum _{n=1}^{\infty }\left(A_n^2+B_n^2\right)}\hspace{1em}\hspace{1em}\Rightarrow \hspace{1em}\hspace{1em}{\displaystyle \sum _{n=1}^{\infty }\frac{1}{{\left(4\hspace{0.17em}{n}^2-1\right)}^2}}=\frac{{\pi }^2-8}{16}\approx \frac{1}{8}$$

The convergence rate is relatively fast (1/n²) since the function is continuous and the first derivative is piecewise continuous. For a real-valued even rectified-sine waveform, A_n = A_−n, and B_n = 0, one obtains:

$$C_0=\frac{2P_{\max }}{\pi }\hspace{1em}\hspace{1em}{C}_n=-\frac{2P_{\max }}{\pi {\left(4n^2-1\right)}^{}}\hspace{1em}\hspace{1em}{C}_n=C_{-n}\hspace{1em}\hspace{1em}\left|C_n\right|=\frac{2P_{\max }}{\pi {\left(4n^2-1\right)}^{}}\hspace{1em}\hspace{1em}{\phi }_n=\pi$$

In phasor form, the Fourier series yields the same result given in Equation (47) earlier:

$$f\left(t\right)=C_0+{\displaystyle \sum _{n=2,4,\dots }^{\infty }2\left|C_n\right|\cos \left(n\omega \hspace{0.17em}t+{\phi }_n\right)}=\frac{2P_{\max }}{\pi }-{\displaystyle \sum _{m=1}^{\infty }\frac{4P_{\max }\cos \left(2m\omega \hspace{0.17em}t\right)}{\pi {\left(4m^2-1\right)}^{}}}$$

The imaginary part of the complex Fourier coefficient is zero and the phase constant is π. Since the infinite series of squares of Fourier coefficients is bounded by the norm of the square-integrable function, the Riemann–Lebesgue lemma [44,45,46,47] assures us that:

$$\underset{\left|n\right|\to \infty }{\lim }{C}_n\to 0\hspace{1em}\mathrm{when}\hspace{1em}\underset{\left|n\right|\to \infty }{\lim }{A}_n\to 0\hspace{1em}\mathrm{and}\hspace{1em}\underset{\left|n\right|\to \infty }{\lim }{B}_n\to 0$$

Using Parseval’s (or Plancherel’s) identity in the complex form delivers the same result as with real form:

$${\displaystyle \sum _{n=-\infty }^{\infty }{\left|C_n\right|}^2{\Vert {\varphi }_n\Vert }^2}={\Vert f\Vert }^2\hspace{1em}\Rightarrow \hspace{1em}{C}_0^2+{\displaystyle \sum _{n=1}^{\infty }\left\{\hspace{0.17em}{\left|C_n\right|}^2+{\left|C_{-n}\right|}^2\right\}}=\frac{P_{\max }^2}{2}\hspace{1em}\Rightarrow \hspace{1em}{\displaystyle \sum _{n=1}^{\infty }\frac{1}{{\left(4\hspace{0.17em}{n}^2-1\right)}^2}}=\frac{{\pi }^2-8}{16}$$

(48)

The first 20 harmonic terms are sufficient to approximate rectified-sine function with high accuracy as the function is smooth everywhere except for the first derivatives being discontinuous at the half-period points.

3.7. Modeling of Wall Oscillations

The uniform wall movement such as in the classic planar CT flow is not easy to achieve in practice and especially not on microscales. Net flow rate induced by non-harmonic oscillatory motion of Newtonian and non-Newtonian fluids in a slit and in the absence of PGs was treated numerically [18,19]. However, in the numerical analysis presented in [18,19], the return stroke was assumed to be so fast that full slip occurred. In this work we are analytically treating, as well as numerically, the motion of the Newtonian fluid with no-slip enforced in the fast return stroke. A simple energy-conservation principle tells us that it makes no difference whether the excitation is harmonic or non-harmonic. All externally applied energies go into the fluid flow (isothermal assumption). However, the idealized flow with no-slip BC considered here does not consider change in physical properties with temperature and the strong non-linear behavior of non-Newtonian fluids. Apparent slip is often being reported in microscale flows [30,31].

The excitation of both wall interfaces is assumed to be in-phase (symmetric). In general, the period T of wall vibrations is not the same as the one for PG oscillations. The illustration of the square-like wall oscillation waveform used here is shown in Figure 3. The return stroke is very fast and of short duration. The area under forward and return strokes must be equal ($L_{fwd}=L_{ret}$) to ensure no net (zero-drift) wall motion:

$${\displaystyle \underset{0}{\overset{aT}{\int }}\left|u_{fwd}\left(t\right)\right|\hspace{0.17em}dt}={\displaystyle \underset{aT}{\overset{T}{\int }}\left|u_{ret}\left(t\right)\right|\hspace{0.17em}dt}\hspace{1em}\hspace{1em}\hspace{1em}{U}_{ret}=-U_{fwd}\left(\frac{a}{1-a}\right)\hspace{1em}\hspace{1em}\hspace{1em}\overline{U}=\frac{U_{fwd}-U_{ret}}{2}$$

(49)

The non-harmonic periodic asymmetric square-waveform function shown in Figure 3 can be expanded in a complex Fourier series with complex coefficients:

$$\begin{array}{l}{C}_0=0\hspace{1em}\hspace{1em}{C}_n=\frac{\overline{U}}{in\pi }\left(1-e^{-ia2n\pi }\right)\hspace{1em}\hspace{1em}{C}_{-n}=-\frac{\overline{U}}{in\pi }\left(1-e^{+ia2n\pi }\right)\\ \left|C_n\right|=\left|C_{-n}\right|=\frac{2\hspace{0.17em}\overline{U}}{n\pi }\sin \left(an\pi \right)\hspace{1em}\hspace{1em}{\varphi }_n=-an\pi \end{array}$$

(50)

Parseval’s identity delivers:

$$2T{\displaystyle \sum _{n=1}^{\infty }{\left|C_n\right|}^2}={\displaystyle \sum _{n=1}^{\infty }\frac{8\hspace{0.17em}{\overline{U}}^{2}T}{n^2{\pi }^2}{\sin }^2\left(an\pi \right)}=U_{fwd}^{2}T\frac{a}{1-a}\hspace{1em}\hspace{1em}\Rightarrow \hspace{1em}\hspace{1em}{\displaystyle \sum _{n=1}^{\infty }\frac{{\sin }^2\left(an\pi \right)}{n^2}}=\frac{{\pi }^{2}a}{2}\left(1-a\right)$$

(51)

which for a = 0.5 gives a well-known analytical result. Clearly this series will contain both cosine and sine harmonics.

By varying the time-fraction factor “a” which automatically defines the return stroke speed and by using the motion-conservation restriction given with Equation (49), we obtain a variety of square-like waveforms. Due to jump discontinuities in square-like excitations, the Gibbs and “ringing” phenomena arise, and the series cannot possibly converge uniformly and specifically not pointwise at jumps no matter how many Fourier coefficients we employ [51,52,53,54,55,56]. Accordingly, the convergence will be quite slower. The Fourier expansion of the above square function requires at least 100 higher-frequency harmonics to approximate the piecewise continuous function sufficiently well.

4. Results and Discussion

An infinite number of particular solutions can be obtained from the general velocity distribution given by Equation (38). The flow rate and the shear rate can be calculated exactly for the known wall and PG excitations using Equations (40) and (41).

4.1. Unsteady Dynamics and Model Verification

Unsteady laminar flow analytical and numerical computations were coded and compiled in Fortran 95 (standard ISO/IEC 1539-1:1997) high-performance optimizing compiler (Lahey Computer Systems^TM) for rapid calculations. Additionally, slower interactive execution utilizing MATLAB^TM 7.14 (R2012a, The MathWorks^TM) interpreter was used for computations, testing and plotting of results.

The numerical (FVM-CFD) and analytical results for sinusoidal oscillatory shear-driven Couette–Stokes (CS) flow with an arbitrary phase ψ and unit amplitude in the absence of static and dynamic PGs and for $\alpha =3.126$ ($\Omega =10$) and, $\alpha =10$ ($\Omega =100$), are shown in Figure 4. In addition to testing the accuracy of the analytical and the numerical model, many useful conclusions can be made regarding the transient behavior and the effect of the WO number on the penetration depth (zone/radius of influence). The analytical solution and the numerical computations agree within 8+ significant digits throughout the entire spatial-temporal domain.

Lower WO number wall oscillations penetrate deeper into the fluid. For WO numbers less than one, there is essentially no transient behavior, and the flow is basically quasi-steady at all times and faithfully responding to moving boundary with zero phase lag. Transient behavior is best observed in the velocity profiles at T/16 and 17T/16, for both, low and high WO numbers. For higher WO numbers, the transient effects are also apparent in the 2nd and 3rd period. For WO = 10 flow, the lower ¾ of fluid domain is basically unaffected by the upper-wall oscillations. In the case of lower WO = 3.26, only the lower ¼ of fluid depth is unaffected by the upper boundary. The analytical result utilizing the shifting function for dimensionless velocity of harmonically oscillating (upper plate) flow in the absence of PGs yields:

$$\begin{array}{l}u\left(y,t;{\alpha }_{\Omega },\psi \right)=\frac{y+1}{2}\sin \left({\alpha }_{\Omega }^{2}t+\psi \right)+2{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n\sin \left[k_n\left(\frac{y+1}{2}\right)\right]\frac{1}{k_n}\frac{{\delta }_n}{\sqrt{1+{\delta }_n^2}}\cos \left({\alpha }_{\Omega }^{2}t+\psi -{\phi }_n\right)}\\ +2{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n\sin \left[k_n\left(\frac{y+1}{2}\right)\right]\frac{1}{k_n}{e}^{-k_n^{2}t}\left[\sin \psi -\frac{{\delta }_n}{\sqrt{1+{\delta }_n^2}}\cos \left(\psi -{\phi }_n\right)\right]}\\ {\delta }_n=\frac{\Omega }{k_n^2}=\frac{{\alpha }_{\Omega }^2}{k_n^2}\hspace{1em}\hspace{1em}{\phi }_n={\tan }^{-1}{\delta }_n\hspace{1em}\hspace{1em}{k}_n=n\pi \hspace{1em}\hspace{1em}-1\le y\le +1\end{array}$$

(52)

In subsequent computations, analytical and computational efforts were only compared in VFR results. This is to avoid busy graphs. The agreement between the theory and the numerical experiment is always perfect. We have successfully performed transient analysis of all basic flows, such as CT, PO, WO, and ST flows which due to the lack of space cannot be shown here. However, more complex combined flows will reveal all the unsteady dynamics.

The CP flow transient is shown in Figure 5. This flow incorporates shear and pressure driving. We already know from the analytical CP flow solution that it will result in stalled or zero steady-state VFR. Interestingly, the shear effect overcomes pressure-driven forces initially. One of the reasons is the smaller inertia of a liquid layer in the wall vicinity. The shear rate approaches the steady-state value of $\dot{\gamma }\left(-1\right)=-3,\hspace{0.17em}\dot{\gamma }\left(1\right)=3$ rapidly. Since the flow is asymmetric, we cannot use cosine eigenfunctions alone. Using the EEM, which converges everywhere except at the moving boundary (non-compatible BCs), we obtain:

$$\begin{array}{l}u\left(y,t;{\alpha }_{\Omega },\psi \right)=-2{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n\sin \left[k_n\left(\frac{y+1}{2}\right)\right]\frac{\sin \left({\alpha }_{\Omega }^{2}t+\psi -{\phi }_n\right)}{k_n\sqrt{1+{\delta }_n^2}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+2{\displaystyle \sum _{n=1}^{\infty }{\left(-1\right)}^n\sin \left[k_n\left(\frac{y+1}{2}\right)\right]\frac{\sin \left(\psi -{\phi }_n\right)}{k_n\sqrt{1+{\delta }_n^2}}{e}^{-k_n^{2}t}}\end{array}$$

(53)

Surprisingly, as it may seem, both series give the same result except at the boundary and the EEM does not converge uniformly. However, it is more general. The unsteady pressure-driven PW flow transient is shown in Figure 6. The pumping action comes from the PG DC-component (0.05) resulting in pulsatile and positive net flow. Flow reversals are similar to pulsatile flow in arterial vessels caused by reciprocating heart pumping. Oscillating shear stresses are essential information in hemodynamic analysis. The classical solution of this problem in terms of elementary functions is given in [4,5,13,14,15], but the solution derived here using Fourier series expansion is perhaps more elegant and more general.

Next, test results treat full transient and QSS solution involving symmetric (ST2) harmonic wall excitations with uniform (CT) components as well as harmonic (WO) PG superimposed onto the steady-state (PO) PG. The analytic solution originates directly from the general solution for the simple wall sinusoidal oscillations. The CPW flow transient is shown in Figure 7. The agreement with CFD is complete. The Womersley component communicates oscillatory fluid behavior.

Oscillatory PSW flow is shown in Figure 8. The interaction of the wall and the PG oscillations results in periodic yet non-harmonic flow. The last verification and testing result is reserved for the most complicated harmonic flow so far, i.e., CPSW flow transient and is shown in Figure 9. Again, no net flow exists as we just added simple harmonic (sinusoidal with phase) ST and WO components to steady CP flow. However, the transient and the QSS flow do not look trivial.

Numerical and analytical computations run under different clocks. While we can choose any future time with the analytical solution, the stability- and accuracy-limited numerical solution has to drudgingly march toward it. Since numerical and analytical computations work under different “clocks”, it was sometimes difficult to synchronize discrete times resulting in small discrepancies. However, we already demonstrated excellent agreement between the numerical experiments and theoretical computations involving generalized Fourier series. Many tens of hours were spent investigating various transient and quasi-steady-state oscillatory flows caused by diverse wall and PG dynamics. An article ten-fold in size would be required to show all relevant tests and comparisons between the analytical and numerical results. It suffices to say that all tests and computations were successful and revealed rich and interesting flow patterns even for such linear flows of Newtonian fluid.

4.2. Results for Unsteady Flows with Nonharmonic Wall Oscillation

In the final set of analytical-only results, we present planar non-harmonic shear driving of Newtonian fluids in the presence of steady and dynamic PGs. Rectified-sine function was used to model PGs (Equation (47)), while the non-harmonic wall oscillations were modeled using square-waveform given by Equation (50). Various alternative harmonic and non-harmonic waveforms could be used, which is left only to our imagination. The Fourier series approach works very well, sometimes even with the most bizarre of functions. Numerical solutions always fully agreed with the analytical results.

A result for symmetric square-waveform shear-driven flow in the absence of any PGs and period of 2π is shown in Figure 10. The fluid responds to wall-dynamics switching logic between forward and return stroke with inertia and resistance. Integral volume, obtained by analytic integration of Equation (40), shows no net fluid flow other than for small initial transient. Non-harmonic shear-driven planar flow with no PGs is shown in Figure 11. No net flow is obtained with non-harmonic driving except for small initial transient. Finally, the non-harmonic shear-driven flow is being exposed to adverse steady-state PG and proverse rectified-sine dynamic PG with finite Womersley number (WO = 3) is shown in Figure 12. The DC-component of proverse pulsatile PG is able to overcome the adverse negative PO pressure component resulting in positive and somewhat interesting oscillating flow patterns. All numerical experiments were in perfect agreement with the analytical model.

Even though Stokes’ penetration layer decreased significantly for the high WO flows, non-harmonic excitation did not result in net flow rates. However, this is chiefly because of the assumed flow idealizations. No energy equation was solved, and the fluid is linear. Since all the external energies applied in the forward and the return strokes go completely into fluid motion, we could not have expected differently. However, high shear rates in the rapid wall return strokes will result in quadratically higher local dissipation and heating of the liquid, which in turn will lower its viscosity, provided that the heat is being conducted away before the forward strokes. Non-Newtonian fluids (such as blood) are especially well suited for non-harmonic shear driving as their dynamic viscosity changes significantly with the shear rate (increasing or decreasing). Special surface roughness (micro-geometry) and/or coatings may be used to increase shear stress in the forward stroke while decreasing it in the reverse stroke. Last, but not least, rapid return stroke may result in partial or full (molecular and/or apparent) wall slip. Work is in progress to solve similar problem in rigid cylindrical geometry by utilizing Newtonian and various non-Newtonian fluids with or without slip which will be reported in subsequent publications.

Works by Rajagopal et al. [61] and Fetecau and Narahari [62] recently came to our attention. They address IBVP problems somewhat related and similar to ours. Rajagopal et al. [61] investigate unsteady flow with pressure-dependent viscosity and the effect of gravity. The hydromagnetic flow of incompressible viscous fluids between horizontal infinite parallel plates embedded in a porous medium has been analytically studied by Fetecau and Narahari [62]. In particular, the Couette and Stokes problem was considered for known magnetic parameters of fluid and porous medium.

Our future studies will include cylindrical and other geometries with and without boundary slip and various other generalized Newtonian fluid models. Some of these problems will be nonlinear for which some special solution techniques will be applied. Numerically, such problems can be solved with ease, but do not provide general solutions.

5. Conclusions

A unified theory of diffusion-controlled laminar slit flow where the driving force was provided by both variable pressure gradients and wall motion was developed. The eigenfunction expansion method was used to solve the global non-homogeneous linear initial and boundary value problem analytically. Simultaneously, a finite-volume method was used for numerical experiment. The classical results for the steady-state, transient and quasi-steady-state flows were recovered both analytically and numerically with high accuracy and fidelity. Both adverse and proverse pressure gradients were used for fluid pressure driving. Harmonically and non-harmonically oscillating infinite rigid impenetrable flat plates were used for fluid shear driving. We used rectified-sine excitations for pressure gradients in combination with the square-waveform non-harmonic zero-drift wall motion. The effect of wall- and pressure-gradient Womersley numbers on flow dynamics has been studied. No net flow rates could be achieved by non-harmonic oscillations for Newtonian fluid, but that is mostly because several physical effects were neglected, and the energy equation was not coupled in. Ongoing analytical, experimental, and numerical studies will attempt to answer some important questions regarding various microscale flow geometries, partial or full apparent boundary slip condition, change in physical parameters due to local viscous dissipation and heat transfer, shear pumping of non-Newtonian fluids, flow stability, interactions and zone-of-influence from multiple moving surfaces and others.