Analytical and Computer Modelling of Chaotic Advection of Stokes Flow Inside a Rectangular Cavity under Constant Velocities

1. Introduction

The mixing problem of laminar fluid flow is important in microtechnology when developing mixing devices in chemical engineering and microbiology. The interest of researchers in this problem arose not only because of the prevalence of the mixing phenomenon in nature, but also due to the use of various mixing technologies in industry and technology (aerodynamics, oceanology, biology, chemistry). Monitoring the conditions of unstable chemical reactions requires the creation of special microchemical processes that maintain the conditions of chemical reactions in a narrow range of physical parameters. They are capable of mixing high molecular weight compounds in a liquid, the molecules of which cannot withstand high stress.

Analysis of various motions allows to conclude that at the initial stage of flow with different geometric scales, convective mixing mechanisms prevail, and the problem of mixing is reduced to analyzing the processes of deformation of the volume of liquid in the selected velocity field. This problem is related to the Lagrangian description of motion in hydrodynamics and is reduced to the analysis of the system of motion of Lagrangian particles in the Eulerian velocity field. Such a problem in the literature has been called the problem of fluid advection.

Liquid advection is associated with the study of the time deformation of selected volumes or regions (in a two-dimensional case) of a liquid that consists of a large number of particles. Despite the fact that such systems have infinitely many degrees of freedom, special attention is paid to the consolidation of individual fluid parts in hydrodynamic flows. Reducing the number of degrees of freedom allows not only to simplify the analysis of the phenomenon, but also to reveal the main properties and patterns in intensive modes of advection.

The study of chaotic mixing regimes arising in some areas of the flow is reduced to the sequential decoupling of two problems. The first task is to determine the velocity field inside a rectangular cavity for a given fluidity at the boundaries. It is necessary to obtain exact solutions, because the second task is associated with finding the trajectory of movement of individual particles of the fluid in the selected area and the deformation of the boundaries of the fluid associated with a dynamically unstable fluid of movement inside the rectangle. To solve the first problem, let’s use the superposition method, which allows to find the asymptotics of unknown coefficients in a number of general solutions, and thereby allows to obtain the desired accuracy.

2. Literature Analisys

In the Lagrangian representation, the problem of advection of a passive marker particle by a prescribed flow defines a dynamical system. For two-dimensional incompressible flow this system is Hamiltonian and has just one degree of freedom. For unsteady flow the system is non-autonomous and one must in general expect to observe chaotic particle motion. These ideas are developed in [1] and subsequently corroborated through the study of a model which provides an idealization of a stirred tank. In the model the fluid is assumed incompressible and inviscid and its motion wholly two-dimensional.

The heat transfer rate from a solid boundary to a highly viscous fluid can be enhanced significantly by a phenomenon which is called chaotic advection or Lagrangian turbulence. Although the flow is laminar and dominated by viscous forces, some fluid particle trajectories are chaotic due either to a suitable boundary displacement protocol or to a change in geometry. As in turbulent flow, the heat transfer rate enhancement between the boundary and the fluid is intimately linked to the mixing of fluid in the system. Chaotic advection in real Stokes flows, i.e., flows governed by viscous forces and that can be constructed experimentally, is reviewed in paper [2].

The Stokes flow inside a two-dimensional rectangular cavity is analyzed in [3] for a highly viscous, incompressible fluid flow, driven by a single rotlet. Specifically, a rigorous solution of the governing two-dimensional biharmonic equation for the stream function is constructed analytically by means of the superposition principle. With this solution, multicellular flow patterns can be described for narrow cavities, in which the number of flow cells is directly related to the value of the aspect ratio.

Two-dimensional chaotic mixing of similar Newtonian fluids in the presence of an advected dissimilar minor phase fluid body with specified size, interfacial tension, and viscosity ratio was numerically investigated in [4]. Interfacial tension was sufficiently high to allow only small deformations in the dissimilar minor phase body. Mixing was confined to a rectangular cavity with periodically driven upper and lower surfaces. Regions of regular motion (i.e., islands) of comparable size to the minor phase body were eventually destroyed or replaced by the minor phase body. Islands persisted for longer times when the initial separation distance between the minor phase body and island was large or when the viscosity ratio was small. When interfacial tension was small enough to deform the minor phase body more readily, islands showed little indication of instability. Results suggest opportunities for improving mixing uniformity in practical processes and disclose how interactions between dissimilar fluids affect mixing.

In [5,6] considered advection of neutrally buoyant discs in two-dimensional chaotic Stokes flow. The goal of the study is to explore a possibility to enhance laminar mixing in batch-flow mixers. Addition of freely moving bodies to periodically driven chaotic flow renders the flowfield nonperiodic [4], i.e., the Lagrangian chaos of the bodies motion induces Eulerian chaos of the flow that makes mixing more intensive. Simulations were performed using a new variant of the immersed boundaries method that allows the direct numerical simulation of fluid–solid flows on a regular rectangular grid without explicit calculation of the forces that the particles exert on the fluid.

The simultaneous effects of flow pulsation and geometrical perturbation on laminar mixing in curved ducts have been numerically studied in [7] by three different metrics: analysis of the secondary flow patterns, Lyapunov exponents and vorticity vector analysis.

In [8] present the paradigm of "designing for chaos" as a general framework for enhancing mixing in microfluidic applications. Designing for chaos is based on a fundamental understanding of the kinematics underlying the mixing process. Computational and experimental analyses demonstrate the effectiveness of the resulting design in generating chaos in the flow.

In paper [9] modified three-dimensional Navier-Stokes system in unbounded domain, which holds the Poincare inequality, is considered. Unique global solvability is proven, for the corresponding semigroup the existence of the global attractor in strong topology of phase space was ascertained, convergence of the obtained attractors to the set of bounded trajectories of the Navier-Stokes 3D-system is shown.

The present a numerical scheme for approximating the incompressible Navier-Stokes equations based on an auxiliary variable associated with the total system energy in [10]. By introducing a dynamic equation for the auxiliary variable and reformulating the Navier-Stokes equations into an equivalent system, the scheme satisfies a discrete energy stability property in terms of a modified energy and it allows for an efficient solution algorithm and implementation.

The purpose of the work is to perform the analysis of the process of advection of the allocated volume of liquid in a rectangular cavity under the action of tangential velocities applied to the upper and lower walls, using the method of piecewise spline interpolation [11]. The use of cubic spline interpolation on smooth contour segments allows to determine the coordinates of additional points on the boundary of the selected volume with high accuracy. The motion of the walls is set by the law, which was considered in the experimental work of Ottino [12].

3. Analytical Determination of the Velocity Field of a Flow of a Viscous Fluid in a Rectangular Cavity

3.1. The Equation of Motion of a Viscous Fluid in the Stokes Approximation

For consideration, let’s take a two-dimensional flow of a viscous incompressible homogeneous fluid which density is equal to $\rho$. Inside the cavity of dimensions $2a\times 2b$, where the upper and lower boundaries are movable and the side walls are stationary, the flow moves with constant kinematic viscosity $\nu$. The Cartesian coordinate system introduced is related to the geometric center of the rectangle (Figure 1). The velocity of the top wall is denoted by $U_{top}\left(x\right)$, and the bottom one is $U_{bot}\left(x\right)$. The problem is to determine the distribution of the velocity field $\overrightarrow{U}(x,y)$ inside the cavity under consideration.

The flow regime is determined by the Reynolds number Re

$$Re={\displaystyle \frac{UL}{\nu }},$$

(1)

where U is a characteristic flow rate $U=max\{U_{top}\left(x\right),U_{bot}\left(x\right)\}$, L - geometric scale of the flow, $L=2a$. For this case, the inertial and viscous terms are proportional to $u^2/L$ and $u\nu /L^2$, respectively. The Reynolds number [13] shows the ratio of inertial forces to viscosity forces.

The velocity index for the flow under consideration must be chosen in such a way that the Reynolds number $Re\ll 1$. This approximation in the available sources [14] is called the Stokes approximation.

The motion of a viscous incompressible fluid, that has constant physical properties, can be described by the Navier - Stokes equation and the continuity equation

$${\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\nabla \overrightarrow{u}=\overrightarrow{F}-{\displaystyle \frac{1}{\rho }}\nabla p+\nu \Delta \overrightarrow{u},$$

(2)

$$\overrightarrow{\nabla }·\overrightarrow{u}=0,$$

(3)

where $\overrightarrow{u}$ - velocity vector, $\overrightarrow{F}$ - the main vector of volume forces, p - pressure.

The left part of equation (2) can be written as

$${\displaystyle \frac{d\overrightarrow{u}}{dt}}={\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\nabla \overrightarrow{u}.$$

(4)

This particle is called the inertial term, it represents the acceleration of a fluid particle or substantial acceleration. This term corresponds to the force of inertia and shows the rate of formation of the momentum of a fluid element on which external forces act. On the right side, the first term is responsible for the volume forces. The second term is responsible for the pressure, the third for the viscosity force.

For small Reynolds values, flows are characteristic that are formed in small-scale geometric regions, with small scales of velocities and flows with high viscosity. When the Re numbers are small, the flows begin to develop under the action of viscosity forces, prevail, then the inertial properties of the fluid can be neglected. The formation of such flows occurs under the action of an applied pressure gradient, or they are affected by the motion of the boundaries of the Stokes approximation to an instantaneous change in the velocity field throughout the entire flow region.

The Navier-Stokes equations in the case when the inertial components are neglected for an incompressible fluid, when there are no external volume forces and, taking into account all conservative mass pressure forces, is simplified to the form:

$$-{\displaystyle \frac{1}{\rho }}\overrightarrow{\nabla }p+\nu \Delta \overrightarrow{u}=0,$$

(5)

together with the continuity equation (3), it is called the Stokes equation.

This equation is linear with respect to the flow velocity field and pressure distribution, and the time component is not included. So the Stokes equations describe a stationary flow of a viscous incompressible fluid, with constant boundary conditions and a stationary pressure gradient.

For the two-dimensional case, the system of Stokes equations can be represented in the form of three differential equations

$$\left\{\begin{array}{c}-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+\nu ({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}})=0,\hfill \\ -{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial y}}+\nu ({\displaystyle \frac{{\partial }^2\nu }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\nu }{\partial y^2}})=0,\hfill \\ {\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial \nu }{\partial y}}=0.\hfill \end{array}\right.$$

(6)

Let’s differentiate the first and second equations of system (6) with respect to x and y, making the addition operation, get rid of the pressure, replace the first two equations of system (6) with one third-order differential equation

$${\displaystyle \frac{{\partial }^{3}u}{\partial y^3}}+{\displaystyle \frac{{\partial }^{3}u}{\partial x^2\partial y}}-{\displaystyle \frac{{\partial }^3\nu }{\partial x^3}}-{\displaystyle \frac{{\partial }^3\nu }{\partial x\partial y^2}}=0.$$

(7)

Let’s introduce a vector stream function $\overrightarrow{\psi }$, by its expression $\overrightarrow{u}=\mathrm{rot}\overrightarrow{\psi }$, that connects the distribution of the velocity field of an incompressible flow [14]. In the two-dimensional case, the vector function $\overrightarrow{\psi }$ has only one component ${\psi }_z=\psi (x,y)$, which connects the components of the velocity field in the following form

$$u={\displaystyle \frac{\partial \psi }{\partial y}},\nu =-{\displaystyle \frac{\partial \psi }{\partial x}}.$$

(8)

Substituting (8) into (7) yields a biharmonic equation with respect to the stream function

$${\displaystyle \frac{{\partial }^4\psi }{\partial x^4}}+2{\displaystyle \frac{{\partial }^4\psi }{\partial x^2\partial y^2}}+{\displaystyle \frac{{\partial }^4\psi }{\partial y^4}}=0.$$

(9)

Taking into account (8), we conclude that the latter equation describes the distribution of the velocity field in the flow that is considered.

Let’s supplement equation (9) with boundary conditions on the lower and upper surfaces. For a viscous fluid, the no-slip conditions [14] on the moving boundary

$$\overrightarrow{u}{|}_S=U_S$$

(10)

can be written in the form of derivatives along the normal and tangent, respectively

$${\left.{\displaystyle \frac{\partial \psi }{\partial n}}\right|}_S=U_S,{\left.{\displaystyle \frac{\partial \psi }{\partial \tau }}\right|}_S{=0\mathrm{or}\psi |}_S=const.$$

(11)

For further calculations, let’s use the biharmonic equation (9) and boundary conditions (11) to construct an analytical solution of the hydrodynamic problem of the flow of a viscous fluid in a rectangular cavity.

3.2. Stokes Flow in a Rectangular Cavity

The two-dimensional Stokes flow of an incompressible viscous fluid can be described as a biharmonic problem. If the movement is too slow, then the inertial forces with a quadratic velocity can be neglected in comparison with the viscous forces, then the stream function can be satisfied with the biharmonic equation

$$\Delta \Delta \psi =0.$$

(12)

Euler components of the velocity vector in rectangular coordinates are determined by (8).

The flow in a rectangular cavity $\left|x\right|<a$, $\left|y\right|<b$ arises with the participation of tangential velocities $U_{+}$ and $U_{-}$ on the upper $(y=b)$ and lower $(y=-b)$ walls, while the side walls remain stationary $x=\pm a$, and therefore, the boundary conditions for equation (12) have the form:

$$\begin{array}{cc}\hfill & \psi =0,{\displaystyle \frac{\partial \psi }{\partial x}}=0,\mathrm{if}x=\pm a,\left|y\right|\le b\hfill \\ \hfill & \psi =0,{\displaystyle \frac{\partial \psi }{\partial y}}=U_{\pm },\mathrm{if}y=\pm b,\left|x\right|\le a.\hfill \end{array}$$

(13)

The velocities $U_{+}$ and $U_{-}$ can both depend on coordinate x and on time t simultaneously. Periodic functions with a definite period T and its dependence on time are of particular interest. So the stream function $\psi$ is a solution to the quasi-stationary limiting problem will also be periodic with a period T.

Using the linear boundary value problem (12), (13), the representation for the stream function in a rectangular cavity can be presented as a sum of two functions ${\psi }_1$ and ${\psi }_2$ describing a symmetric and antisymmetric velocity field. The tangential velocities $U_{+}^{\left(1\right)}=-U_{-}^{\left(1\right)}=U_1$ and $U_{+}^{\left(2\right)}=U_{-}^{\left(2\right)}=U_2$ applied on the upper and lower walls of the cavity excite these fields. In this case, an arbitrary load velocity with even on X functions $U_{+}\left(x\right)$ and $U_{-}\left(x\right)$ can be represented in the form

$$U_1\left(x\right)={\displaystyle \frac{U_{+}\left(x\right)-U_{-}\left(x\right)}{2}},U_2\left(x\right)={\displaystyle \frac{U_{+}\left(x\right)+U_{-}\left(x\right)}{2}}.$$

(14)

3.2.1. Construction of a Solution to a Symmetric Problem with Constant Velocities

Using the superposition method, a solution for a symmetric problem with constant velocities can be constructed as a sum of two Fourier series in the complete orthogonal system of functions [15,16,17]

$$\begin{array}{cc}\hfill {\psi }_{EE}=& b\sum _{m=1}^{\infty }{(-1)}^m{\displaystyle \frac{X_m^{\left(1\right)}}{{\alpha }_m}}(btanh{\alpha }_{m}b{\displaystyle \frac{cosh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-y{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}})cos{\alpha }_{m}x-\hfill \\ \hfill & -a\sum _{l=1}^{\infty }{(-1)}^l{\displaystyle \frac{Y_l^{\left(1\right)}}{{\beta }_l}}(atanh{\beta }_{l}a{\displaystyle \frac{cosh{\beta }_{l}x}{cosh{\beta }_{l}a}}-x{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}})cos{\beta }_{l}y,\hfill \end{array}$$

(15)

where ${\alpha }_m={\displaystyle \frac{(2m-1)\pi }{2a}}$, ${\beta }_l={\displaystyle \frac{(2l-1)\pi }{2b}}$.

The first part of the equation in (15) is the solution for the strip $\left|y\right|\le b$, $-\infty <x<\infty$, while the second part is the solution for the strip $\left|x\right|\le a$, $-\infty <y<\infty$. For the first and second series, there is the required number of degrees of functional freedom in order to satisfy the boundary conditions on two opposite boundaries of the rectangle and to satisfy the requirement of completeness of the system of functions on different sides.

The unknown coefficients of the Fourier series $X_m^{\left(1\right)}$, $Y_l^{\left(1\right)}$ can be calculated by solving the system

$$\left\{\begin{array}{c}{X}_m^{\left(1\right)}b{\Delta }_1\left({\alpha }_{m}b\right)-\sum _{l=1}^{\infty }{Y}_l^{\left(1\right)}{\displaystyle \frac{4{\alpha }_m^2{\beta }_l}{{({\beta }_l^2+{\alpha }_m^2)}^2}}={(-1)}^{m-1}{\alpha }_m{U}_m^{\left(1\right)},1\le m\le \infty ,\hfill \\ Y_l^{\left(1\right)}a{\Delta }_1\left({\beta }_{l}a\right)-\sum _{m=1}^{\infty }{X}_m^{\left(1\right)}{\displaystyle \frac{4{\beta }_l^2{\alpha }_m}{{({\alpha }_m^2+{\beta }_l^2)}^2}}=0,1\le l\le \infty ,\hfill \end{array}\right.$$

(16)

where

$${\Delta }_1\left(\xi \right)=tanh\xi +{\displaystyle \frac{\xi }{{cosh}^2\xi }}$$

(17)

and $U_m^{\left(1\right)}$ — Fourier coefficients of the even function $U_1\left(x\right)$.

Using the boundary conditions: the first equation is satisfied by the boundary condition ${\displaystyle \frac{\partial {\psi }_{EE}}{\partial y}}=\pm U_1\left(x\right)$ at $y=\pm b$, multiplied by $cos{\alpha }_{m}x$ and integrated over x from $-a$ to a, and the second – by fulfilling the boundary condition ${\displaystyle \frac{\partial {\psi }_{EE}}{\partial x}}=0$ at $x=\pm a$ for the desired flow function ${\psi }_{EE}(x,y)$, multiplied by $cos{\beta }_{l}y$ and integrated over y from $-b$ to b.

At infinity, the asymptotics of the coefficients $X_m^{\left(1\right)}$, $Y_l^{\left(1\right)}$ takes the form

$$\begin{array}{cc}\hfill & X_m^{\left(1\right)}={\displaystyle \frac{2{\pi }^2{U}_1\left(\alpha \right)}{ab({\pi }^2-4)}}-{\displaystyle \frac{\pi U_1^{\prime }\left(a\right)}{2ab{\alpha }_m}}+\Re \left(E_{\lambda }^{\left(1\right)}{\alpha }_m^{-{\lambda }_1}\right)+o\left(m^{-3}\right),m\to \infty ,\hfill \\ \hfill & Y_l^{\left(1\right)}={\displaystyle \frac{4\pi U_1\left(a\right)}{ab({\pi }^2-4)}}-{\displaystyle \frac{\pi U_1^{\prime }\left(a\right)}{2ab{\beta }_k}}-\Re \left(E_{\lambda }^{\left(1\right)}{\beta }_l^{-{\lambda }_1}\right)+o\left(l^{-3}\right),l\to \infty \hfill \end{array}$$

(18)

In the case when the function $U_1\left(x\right)$ is continuous together with the first and second derivatives on the interval $\left|x\right|\le a$ we have

$$U_1\left(x\right)-U_1\left(a\right)=\sum _{m=1}^{\infty }{(-1)}^{m-1}{\overline{U}}_{m}cos{\alpha }_{m}x,{\overline{U}}_m=O\left({\alpha }_m^{-3}\right),m\to \infty ,$$

(19)

$$U_m-{\displaystyle \frac{2}{a{\alpha }_m}}{U}_1\left(a\right)={\overline{U}}_m.$$

(20)

Let’s obtain the ratio in this form:

$$U_1^{\prime }\left(a\right)=\sum _{m=1}^{\infty }{\alpha }_m{\overline{U}}_m.$$

(21)

When $U_1\left(x\right)\ne 0$, then with the increase of m the value of $|{\alpha }_m{U}_m|$ does not fall. This condition affects the stability of the numerical solution of an infinite system.

Let’s consider the case when $U_1\left(x\right)=U=const$ on the interval $\left|x\right|\le a$. According to the formula

$$U_m^{\left(1\right)}={\displaystyle \frac{1}{a}}{\int }_{-a}^a{U}_1\left(x\right)cos{\alpha }_{m}x,$$

(22)

the Fourier coefficients for this case take the form $U_m^{\left(1\right)}={\displaystyle \frac{2U}{a{\alpha }_m}}$. Therefore, the peculiarity of the system (16) will be that the free members will grow. Therefore, it is impossible to use existing data in the study of the asymptotic nature of the behavior of unknowns in systems of this type, because free members are not satisfied with the necessary decreasing condition $O\left({\displaystyle \frac{1}{{\alpha }_m^2}}\right)$ [18].

In order to be able to use the method of analysis of regular systems, let’s apply new unknowns $x_m^{\left(1\right)}$ and $y_l^{\left(1\right)}$ [21] according to formulas

$$X_m^{\left(1\right)}=X_1+x_m^{\left(1\right)},Y_l^{\left(1\right)}=Y_1+y_l^{\left(1\right)},$$

(23)

where $X_1$ and $Y_1$ are some constants. Substituting (23) into (16), let’s obtain one more system for determining the unknown coefficients.

$$\begin{array}{cc}\hfill & x_m^{\left(1\right)}b{\Delta }_1\left({\alpha }_{m}b\right)-\sum _{l=1}^{\infty }{y}_l^{\left(1\right)}{\displaystyle \frac{4{\beta }_l{\alpha }_m^2}{{({\alpha }_m^2+{\beta }_l^2)}^2}}=F_m^{\left(1\right)},m=1,2,\cdots \hfill \\ \hfill & y_l^{\left(1\right)}a{\Delta }_1\left({\beta }_{l}a\right)-\sum _{m=1}^{\infty }{x}_m^{\left(1\right)}{\displaystyle \frac{4{\beta }_l^2{\alpha }_m}{{({\alpha }_m^2+{\beta }_l^2)}^2}}=G_l^{\left(1\right)},l=1,2,\cdots \hfill \end{array}$$

(24)

where

$$F_m^{\left(1\right)}=4Y_1{\alpha }_m^2\sum _{l=1}^{\infty }{\displaystyle \frac{{\beta }_l}{{({\alpha }_m^2+{\beta }_l^2)}^2}}-X_{1}b{\Delta }_1\left({\alpha }_{m}b\right)+{\displaystyle \frac{2}{a}}U,$$

$$G_l^{\left(1\right)}=4X_1{\beta }_l^2\sum _{m=1}^{\infty }{\displaystyle \frac{{\alpha }_m}{{({\alpha }_m^2+{\beta }_l^2)}^2}}-Y_{1}a{\Delta }_1\left({\beta }_{l}a\right).$$

Let’s impose a decreasing condition $O\left({\displaystyle \frac{1}{{\alpha }_m^2}}\right)$ on $F_m^{\left(1\right)}$ and $G_l^{\left(1\right)}$. Set the coefficients

$$\begin{array}{cc}\hfill & F_m^{\left(1\right)}=\left({\displaystyle \frac{2b}{\pi }}{Y}_1-b{X}_1+{\displaystyle \frac{2}{a}}U\right)+{\displaystyle \frac{\pi }{6b{\alpha }_m^2}}{Y}_1+O\left(m^{-4}\right),m\to \infty ,\hfill \\ \hfill & G_l^{\left(1\right)}=\left({\displaystyle \frac{2a}{\pi }}{X}_1-a{Y}_1\right)+{\displaystyle \frac{\pi }{6a{\beta }_l^2}}{X}_1+O\left(l^{-4}\right),l\to \infty .\hfill \end{array}$$

(25)

In the equations, the condition of decreasing free terms must be satisfied. As a result, let’s obtain the values of the constants $X_1$ and $Y_1$ as a result of equating the main terms to zero

$$X_1={\displaystyle \frac{2{\pi }^2}{ab({\pi }^2-4)}}U,Y_1={\displaystyle \frac{4\pi }{ab({\pi }^2-4)}}U.$$

(26)

Therefore, the infinite system of linear algebraic equations (24) with respect to unknown coefficients $x_m^{\left(1\right)}$ and $y_l^{\left(1\right)}$, in general, is satisfied by the regularity conditions [18] and for such a system the asymptotic law of solutions is true

$$\underset{m\to \infty }{lim}{x}_m^{\left(1\right)}=\underset{l\to \infty }{lim}{y}_l^{\left(1\right)}=A.$$

(27)

It follows from this condition that there exists a main solution which is bounded and it can be found using the reduction method. The number of unknowns $x_m^{\left(1\right)}$ and $y_l^{\left(1\right)}$ is finite (M and L), if to add all the terms following the first M and L, then their sum will be equal to zero. Therefore, in order to find the unknown coefficients $x_m^{\left(1\right)}$ and $y_l^{\left(1\right)}$ it is necessary to solve a finite system of $M+L$ linear equations.

Next, let’s calculate the unknown stream function ${\psi }_{EE}(x,y)$. Substituting (23) in (15) and obtain the relation

$${\psi }_{EE}(x,y)={\psi }_{EE}^{(a,b)}(x,y)+{\psi }_{EE}^{(x,y)}(x,y),$$

(28)

where

$${\psi }_{EE}^{(a,b)}(x,y)=X_{l}b{S}_a^{\left(1\right)}(x,y)-Y_{1}a{S}_b^{\left(1\right)}(x,y),$$

$${\psi }_{EE}^{(x,y)}(x,y)=b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{x}_m^{\left(1\right)}{p}_m^{\left(1\right)}\left(y\right)cos{\alpha }_{m}x-\alpha \sum _{l=1}^{\infty }{\displaystyle \frac{{(-1)}^l}{{\beta }_l}}{x}_l^{\left(1\right)}{q}_l^{\left(1\right)}\left(x\right)cos{\beta }_{l}y,$$

$$\begin{array}{cc}\hfill & S_a^{\left(1\right)}(x,y)=\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{x}_m^{\left(1\right)}{p}_m^{\left(1\right)}\left(y\right)cos{\alpha }_{m}x,\hfill \\ \hfill & S_b^{\left(1\right)}(x,y)=\sum _{l=1}^{\infty }{\displaystyle \frac{{(-1)}^l}{{\beta }_l}}{x}_l^{\left(1\right)}{q}_l^{\left(1\right)}\left(x\right)cos{\beta }_{l}y.\hfill \end{array}$$

In a rectangular cavity, the terms of the equation $cosh{\alpha }_{m}y/cosh{\alpha }_{m}b$ and $sinh{\alpha }_{m}y/cosh{\alpha }_{m}b$ are quantities of the order $e^{{\alpha }_m\left(\right|y|-1)}$ of that, the series in (24) inside the cavity and on the boundary of the region converge uniformly [15,16,18]. There is no need to improve convergence for these series. But on the border, there are diverging rows $S_a^{\left(1\right)}(x,y)$ and $S_b^{\left(1\right)}(x,y)$.

To transform the expressions $p_m^{\left(1\right)}\left(y\right)$ and $q_l^{\left(1\right)}\left(x\right)$, let’s use the relation

$${\displaystyle \frac{1}{cosh\xi }}=2\sum _{r=0}^{\infty }{(-1)}^r{e}^{-(2r+1)\xi },{\displaystyle \frac{sinh\xi }{{cosh}^2\xi }}=2\sum _{r=0}^{\infty }{(-1)}^r(2r+1)e^{-(2r+1)\xi }$$

(29)

Let’s perform certain transformations to obtain expressions for $S_a^{\left(1\right)}(x,y)$ and $S_b^{\left(1\right)}(x,y)$, namely

$$\begin{array}{cc}\hfill & S_a^{\left(1\right)}(x,y)=\sum _{r=0}^{\infty }{(-1)}^r\left[(2rb+b-y)S_1(x,2rb+b-y)+(2rb+b+y)S_1(x,2rb+b+y)\right],\hfill \\ \hfill & S_a^{\left(1\right)}(x,y)=\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_2(y,2ra+a-x)+(2ra+a+x)S_2(y,2ra+a+x)\right],\hfill \end{array}$$

(30)

where

$$\begin{array}{ccc}\hfill & S_1(\xi ,\eta )=\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{e}^{-{\alpha }_m\eta }cos{\alpha }_m\xi \hfill & \hfill =-{\displaystyle \frac{a}{\pi }}arctan\left[{\displaystyle \frac{cos(\pi \xi /2a)}{sinh(\pi \eta /2a)}}\right],\\ \hfill & S_2(\xi ,\eta )=\sum _{l=1}^{\infty }{\displaystyle \frac{{(-1)}^l}{{\beta }_l}}{e}^{-{\beta }_l\eta }cos{\beta }_l\xi \hfill & \hfill =-{\displaystyle \frac{b}{\pi }}arctan\left[{\displaystyle \frac{cos(\pi \xi /2b)}{sinh(\pi \eta /2b)}}\right].\end{array}$$

Finding the terms slowly coincide, and adding them, let’s obtain the functions $S_1(\xi ,\eta )$ and $S_2(\xi ,\eta )$. Now, in expression (24), the series uniformly coincide both at all points of the region and on the boundary. For the function $S_1(\xi ,\eta )$ and $S_2(\xi ,\eta )$ at the boundary of the region, let’s encountered a singularity, so the calculations were carried out on the basis of the passage to the limit.

Let’s analyze the convergence of the series for the velocity field. Velocity components can be calculated using the formulas given above, as the corresponding derivatives with respect to and from the stream function

$$\begin{array}{cc}\hfill u_{EE}(x,y)=& b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{X}_m^{\left(1\right)}\left(btanh{\alpha }_{m}b{\displaystyle \frac{{\alpha }_{m}sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-y{\displaystyle \frac{{\alpha }_{m}cosh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}\right)cos{\alpha }_{m}x+\hfill \\ \hfill & +a\sum _{l=1}^{\infty }{(-1)}^l{Y}_l^{\left(1\right)}\left(atanh{\beta }_{l}a{\displaystyle \frac{cosh{\beta }_{l}x}{cosh{\beta }_{l}a}}-x{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}\right)sin{\beta }_{l}y=\hfill \\ \hfill =& b\sum _{m=1}^{\infty }{(-1)}^m{X}_m^{\left(1\right)}\left(btanh{\alpha }_{m}b{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-y{\displaystyle \frac{cosh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}\right)cos{\alpha }_{m}x-\hfill \\ \hfill & -b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{X}_m^{\left(1\right)}{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}cos{\alpha }_{m}x-a\sum _{l=1}^{\infty }{(-1)}^l{Y}_l^{\left(1\right)}\left(atanh{\beta }_{l}a{\displaystyle \frac{cosh{\beta }_{l}x}{cosh{\beta }_{l}a}}-x{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}\right)\hfill \\ \hfill sin{\beta }_{l}y,\\ \hfill v_{EE}(x,y)=& b\sum _{m=1}^{\infty }{(-1)}^m{X}_m^{\left(1\right)}\left(btanh{\alpha }_{m}b{\displaystyle \frac{cosh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-y{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}\right)sin{\alpha }_{m}x+\hfill \\ \hfill & +a\sum _{l=1}^{\infty }{\displaystyle \frac{{(-1)}^l}{{\beta }_l}}{Y}_l^{\left(1\right)}\left(atanh{\beta }_{l}a{\displaystyle \frac{{\beta }_{l}sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}-x{\displaystyle \frac{{\beta }_{l}cosh{\beta }_{l}x}{cosh{\beta }_{l}a}}-{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}\right)cos{\beta }_{l}y=\hfill \\ \hfill =& b\sum _{m=1}^{\infty }{(-1)}^m{X}_m^{\left(1\right)}\left(btanh{\alpha }_{m}b{\displaystyle \frac{cosh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}-y{\displaystyle \frac{sinh{\alpha }_{m}y}{cosh{\alpha }_{m}b}}\right)sin{\alpha }_{m}x-\hfill \\ \hfill & -a\sum _{l=1}^{\infty }{\displaystyle \frac{{(-1)}^l}{{\beta }_l}}{Y}_l^{\left(1\right)}{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}cos{\beta }_{l}y+a\sum _{l=1}^{\infty }{(-1)}^l{Y}_l^{\left(1\right)}\left(atanh{\beta }_{l}a{\displaystyle \frac{sinh{\beta }_{l}x}{cosh{\beta }_{l}a}}-x{\displaystyle \frac{cosh{\beta }_{l}x}{cosh{\beta }_{l}a}}\right)\hfill \\ \hfill cos{\beta }_{l}y.\end{array}$$

(31)

Based on formula (28), there are expressions for the velocity component

$$\begin{array}{cc}\hfill & u_{EE}(x,y)=u_{EE}^{(a,b)}(x,y)+u_{EE}^{(x,y)}(x,y),\hfill \\ \hfill & v_{EE}(x,y)=v_{EE}^{(a,b)}(x,y)+v_{EE}^{(x,y)}(x,y),\hfill \end{array}$$

(32)

where $u_{EE}^{(x,y)}(x,y)$ and $v_{EE}^{(x,y)}(x,y)$ is possible to express from formula (31). The coefficients take the form $x_m^{\left(1\right)}$ and $y_l^{\left(1\right)},u_{EE}^{(a,b)}(x,y)$ and $v_{EE}^{(a,b)}(x,y)$ and have the sums of terms with constant coefficients $X_1$ and $Y_1$. Expressions for $u_{EE}^{(a,b)}(x,y)$ and $v_{EE}^{(a,b)}(x,y)$ are as follows

$$\begin{array}{cc}\hfill u_{EE}^{(a,b)}(x,y)=& X_{1}b\sum _{r=0}^{\infty }{(-1)}^r\left[(2rb+b-y)S_3(x,2rb+b-y)-(2rb+b+y)S_3(x,2rb+b+y)\right]-\hfill \\ \hfill & -X_{1}b\sum _{r=0}^{\infty }{(-1)}^r\left[S_1(x,2rb+b-y)-S_1(x,2rb+b+y)\right]+\hfill \\ \hfill & +Y_{1}a\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_4(y,2ra+a-x)+(2ra+a+x)S_4(y,2ra+a+x)\right],\hfill \\ \hfill v_{EE}^{(a,b)}(x,y)=& X_{1}b\sum _{r=0}^{\infty }{(-1)}^r\left[(2rb+b-y)S_5(x,2rb+b-y)+(2rb+b+y)S_5(x,2rb+b+y)\right]-\hfill \\ \hfill & -Y_{1}a\sum _{r=0}^{\infty }{(-1)}^r\left[S_2(y,2ra+a-x)+S_2(y,2ra+a+x)\right]+\hfill \\ \hfill & +Y_{1}a\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_6(y,2ra+a-x)-(2ra+a+x)S_6(y,2ra+a+x)\right],\hfill \end{array}$$

(33)

where

$$S_3(\xi ,\eta )=\sum _{m=1}^{\infty }{(-1)}^m{e}^{-{\alpha }_m\eta }cos{\alpha }_m\xi =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{cos(\pi \xi /2a)cosh(\pi \eta /2a)}{{cos}^2(\pi \xi /2a)+{sinh}^2(\pi \eta /2a)}},$$

$$S_4(\xi ,\eta )=\sum _{l=1}^{\infty }{(-1)}^l{e}^{-{\beta }_l\eta }sin{\beta }_l\xi =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{sin(\pi \xi /2b)sinh(\pi \eta /2b)}{{cos}^2(\pi \xi /2b)+{sinh}^2(\pi \eta /2b)}},$$

$$S_5(\xi ,\eta )=\sum _{m=1}^{\infty }{(-1)}^m{e}^{-{\alpha }_m\eta }sin{\alpha }_m\xi =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{sin(\pi \xi /2a)sinh(\pi \eta /2a)}{{cos}^2(\pi \xi /2a)+{sinh}^2(\pi \eta /2a)}},$$

$$S_6(\xi ,\eta )=\sum _{l=1}^{\infty }{(-1)}^l{e}^{-{\beta }_l\eta }cos{\beta }_l\xi =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{cos(\pi \xi /2b)cosh(\pi \eta /2b)}{{cos}^2(\pi \xi /2b)+{sinh}^2(\pi \eta /2b)}},$$

components $S_1(\xi ,\eta )$ and $S_2(\xi ,\eta )$ are expressed by formulas (30), and $S_1(\xi ,0)=-a/2,S_1(\pm a,\eta )=0,S_2(\xi ,0)=-b/2,S_2(\pm b,\eta )=0$.

The series that are included in the solution coincide quickly, so it is possible to select a part to sum to the final values. The accuracy of the fulfillment of the boundary conditions is associated with an infinite number of members of the series, when using the reduced series, there is an error in solving the boundary value problem.

3.2.2. Construction of a Solution to an Antisymmetric Problem with Constant Velocities

Similarly to how let’s found the general solution ${\psi }_{EE}(x,y)$ using the superposition method, let’s find the solution to the boundary value problem [15,16]

$$\begin{array}{cc}\hfill {\psi }_{EO}=& b\sum _{m=1}^{\infty }{(-1)}^m{\displaystyle \frac{X_m^{\left(2\right)}}{{\alpha }_m}}\left(bcoth{\alpha }_{m}b{\displaystyle \frac{sinh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}-y{\displaystyle \frac{cosh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}\right)cos{\alpha }_{m}x-\hfill \\ \hfill & -a\sum _{k=1}^{\infty }{(-1)}^k{\displaystyle \frac{Y_k^{\left(2\right)}}{{\gamma }_k}}\left(a\tan \mathrm{h}{\gamma }_{k}a{\displaystyle \frac{cosh{\gamma }_{k}x}{cosh{\gamma }_{k}a}}-x{\displaystyle \frac{sinh{\gamma }_{k}x}{sinh{\gamma }_{k}a}}\right)sin{\gamma }_{k}y\hfill \end{array}$$

(34)

where ${\gamma }_k={\displaystyle \frac{k\pi }{b}}$.

When the boundary conditions are satisfied, let’s obtain a system of linear algebraic equations for the unknown Fourier coefficients $X_m^{\left(2\right)}$ and $Y_k^{\left(2\right)}$ and write it in the form

$$\begin{array}{cc}\hfill & X_m^{\left(2\right)}b{\Delta }_2\left({\alpha }_{m}b\right)-\sum _{k=1}^{\infty }{Y}_k^{\left(2\right)}{\displaystyle \frac{4{\delta }_k{\alpha }_m^2}{{({\alpha }_m^2+{\delta }_k^2)}^2}}={\alpha }_m{U}_m^{\left(2\right)},m=1,2,\cdots ,\hfill \\ \hfill & Y_k^{\left(2\right)}a{\Delta }_1\left({\delta }_{k}a\right)-\sum _{m=1}^{\infty }{X}_m^{\left(2\right)}{\displaystyle \frac{4{\delta }_k^2{\alpha }_m}{{({\alpha }_m^2+{\delta }_k^2)}^2}}=0,l=1,2,\cdots ,\hfill \end{array}$$

(35)

$${\Delta }_2\left(\xi \right)=coth\xi -{\displaystyle \frac{\xi }{{sinh}^2\xi }}$$

(36)

$U_m^{\left(2\right)}$ — Fourier coefficients of an even function $U_2\left(x\right)$.

Let’s suppose that $U_2\left(x\right)=U=const$ in between $\left|x\right|\le a$. The Fourier coefficients in this case will have the form

$$U_m^{\left(2\right)}={\displaystyle \frac{{(-1)}^{m-1}}{a}}{\int }_{-a}^a{U}_2\left(x\right)cos{\alpha }_{m}xdx={\displaystyle \frac{2U}{a{\alpha }_m}}.$$

(37)

To neglect the influence of a non-decreasing free term, let’s introduce new variables $x_m^{\left(2\right)}$ and $y_k^{\left(2\right)}$

$$X_m^{\left(2\right)}=X_2+x_m^{\left(2\right)},Y_k^{\left(2\right)}=Y_2+y_k^{\left(2\right)},$$

(38)

At infinity, the asymptotic of the coefficients $X_m^{\left(2\right)}$ and $Y_k^{\left(2\right)}$ takes the form

$$\begin{array}{cc}\hfill & x_m^{\left(2\right)}b{\Delta }_2\left({\alpha }_{m}b\right)-\sum _{k=1}^{\infty }{y}_k^{\left(2\right)}{\displaystyle \frac{4{\delta }_k{\alpha }_m^2}{{({\alpha }_m^2+{\delta }_k^2)}^2}}=F_m^{\left(2\right)},m=1,2,\cdots \hfill \\ \hfill & y_k^{\left(2\right)}a{\Delta }_1\left({\delta }_{k}a\right)-\sum _{m=1}^{\infty }{x}_m^{\left(2\right)}{\displaystyle \frac{4{\delta }_k^2{\alpha }_m}{{({\alpha }_m^2+{\delta }_k^2)}^2}}=G_k^{\left(2\right)},l=1,2,\cdots \hfill \end{array}$$

(39)

In the right part there are

$$\begin{array}{cc}\hfill & F_m^{\left(2\right)}=4Y_2{\alpha }_m^2\sum _{k=1}^{\infty }{\displaystyle \frac{{\delta }_k}{{({\alpha }_m^2+{\delta }_k^2)}^2}}-X_{2}b{\Delta }_2\left({\alpha }_{m}b\right)+{\displaystyle \frac{2}{a}}U,\hfill \\ \hfill & G_k^{\left(2\right)}=4X_2{\delta }_k^2\sum _{m=1}^{\infty }{\displaystyle \frac{{\alpha }_m}{{({\alpha }_m^2+{\delta }_k^2)}^2}}-Y_{2}a{\Delta }_1\left({\delta }_{k}a\right).\hfill \end{array}$$

(40)

For $F_m^{\left(2\right)}$ and $G_k^{\left(2\right)}$, the decrease requirement $O\left({\displaystyle \frac{1}{{\alpha }_m^2}}\right)$ is performed. For system (39), the law of asymptotic expressions can be applied, since it satisfies the regularity conditions [18].

$$\underset{m\to \infty }{lim}{x}_m^{\left(2\right)}=\underset{k\to \infty }{lim}{y}_k^{\left(2\right)}=B$$

(41)

The solution is sought by a simple reduction method. The number of unknowns $x_m^{\left(2\right)}$ and $y_k^{\left(2\right)}$ is finite (M and K), if to add all the terms following the first and then their sum is equal to zero. To find the unknown coefficients $x_m^{\left(2\right)}$ and $y_k^{\left(2\right)}$ it is necessary to solve a finite system of linear equations.

Let’s also take $x_m^{\left(2\right)}$ and $y_k^{\left(2\right)}$ for the known values. Let’s calculate the stream function

$${\psi }_{EO}(x,y)={\psi }_{EO}^{(a,b)}(x,y)+{\psi }_{EO}^{(x,y)}(x,y),$$

(42)

where

$$\begin{array}{cc}\hfill & {\psi }_{EO}^{(x,y)}(x,y)=b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^n}{{\alpha }_m}}{x}_m^{\left(2\right)}{p}_m^{\left(2\right)}\left(y\right)cos{\alpha }_{m}x-a\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\delta }_k}}{y}_k^{\left(2\right)}{q}_k^{\left(2\right)}\left(x\right)sin{\delta }_{k}y,\hfill \\ \hfill & {\psi }_{EO}^{(a,b)}(x,y)=X_{2}b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{p}_m^{\left(2\right)}\left(y\right)cos{\alpha }_{m}x-Y_{2}a\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\delta }_k}}{q}_k^{\left(2\right)}\left(x\right)sin{\delta }_{k}y.\hfill \end{array}$$

Let’s use the ratios

$${\displaystyle \frac{1}{sinh\xi }}=2\sum _{r=0}^{\infty }{e}^{-(2r+1)\xi },{\displaystyle \frac{cosh\xi }{{sinh}^2\xi }}=2\sum _{r=0}^{\infty }(2r+1)e^{-(2r+1)\xi },$$

(43)

to transform $p_m^{\left(2\right)}\left(y\right)$ and $q_k^{\left(2\right)}\left(x\right)$. Performing mathematical transformations, let’s find an expression for ${\psi }_{EO}^{(a,b)}(x,y)$

$$\begin{array}{cc}\hfill {\psi }_{EO}^{(a,b)}(x,y)=& X_{2}b\sum _{r=0}^{\infty }\left[(2rb+b-y)S_1(x,2rb+b-y)-(2rb+b+y)S_1(x,2rb+b+y)\right]-\hfill \\ \hfill & -Y_{2}a\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_{14}(y,2ra+a-x)+(2ra+a+x)S_{14}(y,2ra+a+x)\right],\hfill \end{array}$$

(44)

where

$$S_{14}(\xi ,\eta )=\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\delta }_k}}{e}^{-{\delta }_k\eta }sin{\delta }_k\xi =-{\displaystyle \frac{b}{\pi }}arctan\left[{\displaystyle \frac{sin(\pi \xi /b)}{exp(\pi \eta /b)+cos(\pi \xi /b)}}\right].$$

The velocity components, as derivatives with respect to x and y of the stream function, can be represented as

$$\begin{array}{cc}\hfill & u_{EO}(x,y)=u_{EO}^{(a,b)}(x,y)+u_{EO}^{(x,y)}(x,y),\hfill \\ \hfill & v_{EO}(x,y)=v_{EO}^{(a,b)}(x,y)+v_{EO}^{(x,y)}(x,y),\hfill \end{array}$$

(45)

where $u_{EO}^{(x,y)}(x,y)$ and $v_{EO}^{(x,y)}(x,y)$ are expressed from the following formulas with coefficients $x_m^{\left(2\right)}$ and $y_k^{\left(2\right)}$

$$\begin{array}{cc}\hfill u_{EO}^{(x,y)}(x,y)=& b\sum _{m=1}^{\infty }{(-1)}^m{x}_m^{\left(2\right)}\left(bcoth{\alpha }_{m}b{\displaystyle \frac{cosh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}-y{\displaystyle \frac{sinh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}\right)cos{\alpha }_{m}x-\hfill \\ \hfill & -b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{x}_m^{\left(2\right)}{\displaystyle \frac{cosh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}cos{\alpha }_{m}x-a\sum _{k=1}^{\infty }{(-1)}^k{y}_k^{\left(2\right)}{q}_k^{\left(2\right)}\left(x\right)cos{\delta }_{k}y,\hfill \\ \hfill v_{EO}^{(x,y)}(x,y)=& a\sum _{k=1}^{\infty }{(-1)}^k{y}_k^{\left(2\right)}\left(atanh{\delta }_{k}a{\displaystyle \frac{sinh{\delta }_{k}x}{cosh{\delta }_{k}a}}-x{\displaystyle \frac{cosh{\delta }_{k}x}{cosh{\delta }_{k}a}}\right)sin{\delta }_{k}y-\hfill \\ \hfill & -a\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\delta }_k}}{y}_k^{\left(2\right)}{\displaystyle \frac{sinh{\delta }_{k}x}{cosh{\delta }_{k}a}}sin{\delta }_{k}y+b\sum _{m=1}^{\infty }{(-1)}^m{x}_m^{\left(2\right)}{p}_m^{\left(2\right)}\left(y\right)sin{\alpha }_{m}x,\hfill \end{array}$$

(46)

and $u_{EO}^{(a,b)(x,y)}$ and $v_{EO}^{(a,b)}(x,y)$ have the sums of all components with constant coefficients $X_2$ and $Y_2$, they can be calculated as derivatives with respect to x and y of the component of the stream function ${\psi }_{EO}^{(a,b)}(x,y)$,

$$\begin{array}{cc}\hfill u_{EO}^{(a,b)}(x,y)=& X_{2}b\sum _{m=1}^{\infty }{(-1)}^m\left(bcoth{\alpha }_{m}b{\displaystyle \frac{cosh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}-y{\displaystyle \frac{sinh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}\right)cos{\alpha }_{m}x-\hfill \\ \hfill & -X_{2}b\sum _{m=1}^{\infty }{\displaystyle \frac{{(-1)}^m}{{\alpha }_m}}{\displaystyle \frac{cosh{\alpha }_{m}y}{sinh{\alpha }_{m}b}}cos{\alpha }_{m}x-Y_{2}a\sum _{k=1}^{\infty }{(-1)}^k{q}_k^{\left(2\right)}\left(x\right)cos{\delta }_{k}y=\hfill \\ \hfill =& X_{2}b\sum _{r=0}^{\infty }\left[(2rb+b-y)S_3(x,2rb+b-y)+(2rb+b+y)S_3(x,2rb+b+y)\right]-\hfill \\ \hfill & -X_{2}b\sum _{r=0}^{\infty }\left[S_1(x,2rb+b-y)+S_1(x,2rb+b+y)\right]-\hfill \\ \hfill & -Y_{2}a\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_{15}(y,2ra+a-x)+(2ra+a+x)S_{15}(y,2ra+a+x)\right],\hfill \end{array}$$

(47)

where

$$S_{15}(\xi ,\eta )=\sum _{k=1}^{\infty }{(-1)}^k{e}^{-{\delta }_k\eta }cos{\delta }_{k\xi }={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{sinh(\pi \eta /b)}{cosh(\pi \eta /b)+cos(\pi \xi /b)}}-1\right]$$

$$\begin{array}{cc}\hfill v_{EO}^{(a,b)}(x,y)=& Y_{2}a\sum _{k=1}^{\infty }{(-1)}^k\left(a\tan \mathrm{h}{\delta }_{k}a{\displaystyle \frac{sinh{\delta }_{k}x}{cosh{\delta }_{k}a}}-x{\displaystyle \frac{cosh{\delta }_{k}x}{cosh{\delta }_{k}a}}\right)sin{\delta }_{k}y-\hfill \\ \hfill & -Y_{2}a\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\delta }_k}}{\displaystyle \frac{sinh{\delta }_{k}x}{cosh{\delta }_{k}a}}sin{\delta }_{k}y+X_{2}b\sum _{m=1}^{\infty }{(-1)}^m{p}_m^{\left(2\right)}\left(y\right)sin{\alpha }_{m}x=\hfill \\ \hfill =& X_{2}b\sum _{r=0}^{\infty }\left[(2rb+b-y)S_5(x,2rb+b-y)-(2rb+b+y)S_5(x,2rb+b+y)\right]-\hfill \\ \hfill & -Y_{2}a\sum _{r=0}^{\infty }\left[S_{14}(y,2ra+a-x)-S_{14}(y,2ra+a+x)\right]+\hfill \\ \hfill & +Y_{2}a\sum _{r=0}^{\infty }{(-1)}^r\left[(2ra+a-x)S_{16}(y,2ra+a-x)-(2ra+a+x)S_{16}(y,2ra+a+x)\right]\hfill \end{array}$$

(48)

where

$$S_{16}(\xi ,\eta )=\sum _{k=1}^{\infty }{(-1)}^k{e}^{-{\delta }_k\eta }sin{\delta }_k\xi =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{sin(\pi \xi /b)}{cosh(\pi \eta /b)+cos(\pi \xi /b)}}.$$

With a constant velocity distribution on the boundary, the stream function ${\psi }_{EO}(x,y)$ can be calculated by combining formulas (42) and (44). Velocity components $u_{EO}(x,y)$ and $v_{EO}(x,y)$ are calculated by combining (45) - (48).

4. Advection of Fluid

4.1. Advection Equations of a Moving Fluid Particle

Any small, non-inert fluid particle moves with a velocity equal to the flow velocity at the point at which it is located. This means that the equality

$$V_{particle}=V_{flow}$$

(49)

and this will be the formal expression of advection. This equation is defining in Lagrangian hydromechanics [14]: the kinematics of the fluid itself is such that each particle of the fluid undergoes passive advection.

The velocity of a particle in the two-dimensional case $V_{particle}$ is, of course, given by the rate of change in its position:

$${\overrightarrow{V}}_{particle}=\left({\displaystyle \frac{dx}{dt}},{\displaystyle \frac{dy}{dt}}\right),$$

(50)

where $(x,y)$ — coordinates of the radius vector of the particle in the Cartesian coordinate system.

The fluid velocity is specified according to other considerations, including solving some system of partial differential equations, for example, the Euler equations, the Navier-Stokes equations, or the Stokes equations. In other words, the analysis of the motion of individual fluid particles provides that a hydrodynamic problem has been solved, which allows for an arbitrary point at an arbitrary time instant to determine the magnitude of the velocity vector of the flow under consideration

$$V_{fluid}=(u(x,y,t),v(x,y,t)).$$

(51)

Then condition (49) leads to a system of ordinary differential equations, which are called advection equations [19,20]

$${\displaystyle \frac{dx}{dt}}=u(x,y,t),{\displaystyle \frac{dy}{dt}}=v(x,y,t),$$

(52)

In other words, equation (52) describes the motion of a Lagrangian liquid particle in an Eulerian velocity field [14].

From the point of view of the theory of dynamical systems, two ordinary differential equations (52) are more than enough to obtain a non-integral or chaotic dynamical system. The right-hand sides of the equations don’t even need to be very complicated.

In two-dimensional space, to obtain a chaotic motion of a particle, it is necessary that the flow depends on time. Stationary 2D advection is integrated.

From the theoretical point of view, the greatest interest from the point of view of the randomness of the motion of individual fluid particles in hydrodynamic systems are two-dimensional flows, since in this case the velocity is the derivative of the stream function $\psi$ and is expressed through formulas.

If equation (8) is combined with advection equations (52), then they turn into Hamiltonian canonical equations for a system with one degree of freedom

$${\displaystyle \frac{dx}{dt}}={\displaystyle \frac{\partial \psi }{\partial y}},{\displaystyle \frac{dy}{dt}}=-{\displaystyle \frac{\partial \psi }{\partial x}},$$

(53)

with initial conditions $x=x_0,y=y_0$ at $t=0$.

Kinematic equations (53) can be interpreted as a Hamiltonian dynamical system with one degree of freedom, where both play the role of canonical variables and are conjugate coordinates. Any of these coordinates can be taken as a generalized coordinate. Then the second Cartesian coordinate will be the conjugate generalized impulse. The stream function $\psi$ plays the role of the Hamiltonian. The phase space in this problem is configuration space.

The system of differential equations (53) describes the motion of passive non-inertial particles in the two-dimensional case in terms of the stream function, regardless of whether the fluid is viscous or not. And here there is neither contradiction nor paradox, because the Hamiltonian nature of kinematics arises precisely from incompressibility. This does not depend on whether the motion is dynamically dissipative or not. To determine the position of the particles, it is necessary to integrate the equations of motion (53) of each particle within a certain time interval.

4.2. Numerical Modeling of Closed Loop Advection Process

After finding the velocity field in a rectangular cavity, one can begin to study the problem of mixing a viscous fluid with slow motions of the Stokes flow. Below there are examples of calculations obtained on the basis of understanding the analytical velocity field. The periodic motion of the rectangular cavity walls is given by

$$\begin{array}{cc}\hfill U_{top}=& \left\{\begin{array}{cc}{U}_0,\hfill & nT\le t<\left(n+{\displaystyle \frac{1}{2}}\right)T\hfill \\ 0,\hfill & \left(n+{\displaystyle \frac{1}{2}}\right)T\le t<(n+1)T\hfill \end{array}\right.\hfill \\ \hfill U_{bot}=& \left\{\begin{array}{cc}0,\hfill & nT\le t<\left(n+{\displaystyle \frac{1}{2}}\right)T\hfill \\ -U_0,\hfill & \left(n+{\displaystyle \frac{1}{2}}\right)T\le t<(n+1)T\hfill \end{array}\right.\hfill \end{array}$$

(54)

Let the upper wall move from left to right with velocity $U_0$ during the first half of the period, the lower wall is motionless, and the second half of the period, the lower wall moves from right to left with velocity $-U_0$, and the upper wall is motionless.

In the future, it is convenient to normalize the problem to the height of the cavity $2b$ and to the period T of motion of the boundaries of the rectangular cavity. Let’s introduce a dimensionless parameter $D={\displaystyle \frac{U_{0}T}{2a}}$ characterizing the movement of each of the movable walls, where $U_0$ — modulus of the maximum wall velocity.

As an example, let’s consider the mixing process in a rectangular cavity with parameters $U_0$ and D corresponding to the data of the Ottino experiment [12].

Let’s consider the process of advection of a selected volume of liquid in the form of a circular region of radius $r=0.1$ centered at the origin, the dimensions of which coincide with the dimensions of the spot from the Ottino experiment. The initial position of the circular spot is shown in Figure 2(a).

In order to numerically simulate the fluid mixing process depicted in Figure 2(f) at $D=6.24$, the following procedure was chosen. Let’s fix the value $U_0=10.24$ and increase T from 0 to $8T$ (corresponds to the value $D=6.24$) and monitor the advection process of the selected spot. This approach minimizes the accumulation of error, leads to an erroneous result in the case of a sufficiently large value $U_0$ and consideration of the advection process during only one period T.

It should be noted that the process of advection of the contour is based on piecewise spline cubic interpolation method, developed by Prof. A. Gurzhiy. [11].

The images obtained illustrate that numerical calculations in the process of solving the problem quite accurately illustrate the mixing process. The results are in agreement with the experiment carried out by Ottino [12]. As in the experiment, we obtained the Period-1 island, where the mixing process is not take place (Ottino’s article [12] demonstrates this phenomenon).

5. Conclusions

Using the superposition method, an analytical distribution of the velocity field of the flow of a viscous fluid moving in a rectangular cavity under the action of constant tangential velocities applied to the upper and lower walls is obtained.

The process of mixing the selected volume of fluid with periodic movement of the cavity walls was numerically modeled using the contour advection technique developed by Prof. A. Gurzhiy, and the results obtained are in good agreement with the experimental data.