Matrix logarithmic wave equation and multi-channel systems in fluid mechanics

TL;DR

This paper develops a mathematical framework linking nonlinear wave equations with matrix logarithmic nonlinearity to multi-channel fluid flow systems, providing analytical solutions and insights into complex fluid behaviors.

Contribution

It introduces a novel mapping between matrix logarithmic wave equations and multi-channel fluid flow equations, extending previous models to include matrix and multi-component systems.

Findings

Derived Gaussian-type matrix solutions for specific cases

Mapped nonlinear wave equations to multi-channel fluid flow models

Provided analytical insights into complex fluid behaviors

Abstract

We formulate the mapping between a large class of nonlinear wave equations and flow equations for barotropic fluid with internal surface tension and capillary effects. Motivated by statistical mechanics and multi-channel physics arguments, we focus on wave equations with logarithmic nonlinearity, and further generalize them to matrix equations. We map the resulting equation to flow equations of multi-channel or multi-component Korteweg-type materials. For some special cases, we analytically derive Gaussian-type matrix solutions and study them in the context of fluid mechanics.

Equations

$$\Psi =\sqrt{\rho }\exp (iS),$$

$${\int }_V|\Psi {|}^{2}dV={\int }_V\rho dV=M>0,$$

$$i{\partial }_t\Psi =\left[-\frac{\tilde{D}}{2}{\nabla }^2-b\ln \left(|\Psi {|}^2/{\rho }_0\right)\right]\Psi ,$$

$$b\sim T\sim T_{\Psi },$$

$$i{\partial }_t\Psi =\left[-\frac{\tilde{D}}{2}{\nabla }^2-F(|\Psi {|}^2)\right]\Psi ,$$

$$\displaystyle\begin{array}[]{l}\partial_{t}\rho+{\tilde{D}}\boldsymbol{\nabla}\cdot(\rho\boldsymbol{\nabla}S)=0,\\ \partial_{t}S-\frac{{\tilde{D}}}{2}\left[\boldsymbol{\nabla}\cdot\left(\frac{\boldsymbol{\nabla}\sqrt{\rho}}{\sqrt{\rho}}\right)+\frac{(\boldsymbol{\nabla}\sqrt{\rho})^{2}}{\rho}-(\boldsymbol{\nabla}S)^{2}\right]=F(\rho),\end{array}$$

$$\text{u}=\tilde{D}\nabla S,$$

$$\displaystyle\begin{array}[]{l}\partial_{t}\rho+\boldsymbol{\nabla}\cdot(\rho\textbf{u})=0,\\ \partial_{t}\textbf{u}+\textbf{u}\cdot\boldsymbol{\nabla}\textbf{u}-\frac{1}{\rho}\boldsymbol{\nabla}\cdot T=0,\end{array}$$

$$T=-\frac{{\tilde{D}}^2}{4\rho }\nabla \rho \otimes \nabla \rho -\tilde{p}I,$$

$$\tilde{p}=p(\rho )-\frac{1}{4}{\tilde{D}}^2{\nabla }^2\rho ,$$

$$p(\rho )=-\tilde{D}\int \rho F^{\prime }(\rho )d\rho =\tilde{D}\left[\int F(\rho )d\rho -\rho F(\rho )\right],$$

$$c_s^2=\frac{\partial \tilde{p}}{\partial \rho }=-\tilde{D}\rho F^{\prime }(\rho )-\frac{{\tilde{D}}^2}{4}\frac{\partial }{\partial \rho }{\nabla }^2\rho ,$$

$$\tilde{p}\approx -\tilde{D}\int \rho F^{\prime }(\rho )d\rho ,c_s\approx \sqrt{-\tilde{D}\rho F^{\prime }(\rho )}.$$

$$F(\rho )=b\ln \rho ,$$

$$p(\rho )=-\tilde{D}b\rho ,$$

$$\tilde{p}\approx -\tilde{D}b\rho \propto (T-T_0)\rho ,c_s\approx \sqrt{-\tilde{D}b}\propto \sqrt{T-T_0},$$

$$\Psi =\hat{\mathcal{M}}{\Psi }_0,$$

$$i{\partial }_t\mathcal{M}+\frac{\tilde{D}}{2}{\nabla }^2\mathcal{M}+b\ln \left({\rho }_0^{-1}\mathcal{M}{\mathcal{M}}^{†}\right)\mathcal{M}=0,$$

$$\mathrm{Tr}({\mathcal{M}}^{†}\mathcal{M})\equiv ⟨\mathcal{M}|\mathcal{M}⟩={\int }_V\mathrm{tr}({\mathcal{M}}^{†}\mathcal{M})dV=M,$$

$$\rho =\mathrm{tr}({\mathcal{M}}^{†}\mathcal{M}),$$

$$\mathcal{M}=\mathcal{R}\exp (i\mathcal{S}),$$

$$\rho =\mathrm{tr}({\mathcal{R}}^2),\text{u}=\tilde{D}\mathrm{tr}(\nabla \mathcal{S}),$$

$$\displaystyle\begin{array}[]{l}\partial_{t}{\cal R}+\boldsymbol{\nabla}{\cal R}\cdot\boldsymbol{\cal U}+\frac{1}{2}{\cal R}\boldsymbol{\nabla}\cdot\boldsymbol{\cal U}=0,\\ \partial_{t}{\cal S}+\frac{1}{2{\tilde{D}}}\boldsymbol{\cal U}^{2}-\frac{{\tilde{D}}}{2}{\cal R}^{-1}\boldsymbol{\nabla}^{2}{\cal R}-b\ln\left(\rho_{0}^{-1}{\cal R}^{2}\right)=0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\partial_{t}({\cal R}^{2})+\boldsymbol{\nabla}\cdot({\cal R}^{2}\boldsymbol{\cal U})=0,\\ \partial_{t}\boldsymbol{\cal U}+(\boldsymbol{\cal U}\cdot\boldsymbol{\nabla})\,\boldsymbol{\cal U}-\frac{{\tilde{D}}^{2}}{2}\boldsymbol{\nabla}\left({\cal R}^{-1}\boldsymbol{\nabla}^{2}{\cal R}\right)\\ \qquad-b{\tilde{D}}{\cal R}^{-2}\boldsymbol{\nabla}\left({\cal R}^{2}\right)=0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\partial_{t}\rho+\boldsymbol{\nabla}\cdot{\rm tr}({\cal R}^{2}\boldsymbol{\cal U})=0,\\ \frac{D\textbf{u}}{Dt}\equiv\partial_{t}\textbf{u}+{\rm tr}\left[(\boldsymbol{\cal U}\cdot\boldsymbol{\nabla})\,\boldsymbol{\cal U}\right]=\textbf{g},\end{array}$$

$$\text{g}=\mathrm{tr}\mathcal{G}$$

$$\mathcal{G}=\frac{{\tilde{D}}^2}{2}\nabla \left({\mathcal{R}}^{-1}{\nabla }^2\mathcal{R}\right)+b\tilde{D}{\mathcal{R}}^{-2}\nabla \left({\mathcal{R}}^2\right)$$

$$i{\partial }_t\mathcal{M}+\frac{1}{2}\tilde{D}{\partial }_{xx}\mathcal{M}+b\ln ({\rho }_0^{-1}\mathcal{M}{\mathcal{M}}^{†})\mathcal{M}=0,$$

$$\underset{x_L}{\overset{x_R}{\int }}\mathrm{tr}({\mathcal{M}}^{†}\mathcal{M})dx=M,$$

$$\mathcal{M}$$

Full text

††preprint: J. Theor. Appl. Mech. 57, 843-852 (2019)  [DOI: 10.15632/jtam-pl/112063]

I Introduction

Wave equations with logarithmic nonlinearity began to gain a considerable interest among physicists since works by Rosen (1968) and Bialynicki-Birula and Mycielski (1976). The corresponding models have been proven to be instrumental in dealing with physics of quantum fields and particles (Rosen, 1968; Rosen, 1969; Bialynicki-Birula and Mycielski, 1979; Zloshchastiev, 2010; Zloshchastiev, 2011) nonlinear generalizations of quantum mechanics (Bialynicki-Birula and Mycielski, 1976) optics and transport or diffusion phenomena, nuclear physics, theory of dissipative systems and quantum information (Yasue, 1978; Brasher, 1991; Znojil et al., 2017, Zloshchastiev, 2018a), theory of quantum liquids and superfluidity (Zloshchastiev, 2011; Avdeenkov and Zloshchastiev, 2011; Zloshchastiev, 2012; Bouharia, 2015; Zloshchastiev, 2017) and theory of physical vacuum and classical and quantum gravity. (Zloshchastiev, 2010; Zloshchastiev, 2011; Scott et al., 2016).

Special attention was paid to the applications of these equations in fluid mechanics, including the classical hydrodynamics of Korteweg-type materials, such as the barotropic Korteweg capillary fluids (De Martino et al., 2003; Lauro, 2008; Zloshchastiev, 2018b). In such a class of models, often referred as diffuse interface models, the capillary interface is viewed as a diffusion transition domain of rapid smooth variation of density, while surface tension is intrinsically incorporated (Dunn and Serrin, 1985; Dell’Isola F and Kosiński, 1993; Anderson et al., 1998). This allows us to describe flows with a spontaneous nucleation, coalescence and breakdown of density inhomogeneities in two-phase systems (Antanovskii, 1996). In this paper, we continue our studies in this direction, focusing on constructing models which can be used to describe a variety of multi-channel phenomena occuring in mechanics, such as multi-channel scattering, multi-component fluids and capillary microfluidic arrays and chips.

The paper is structured as follows. In section II we give a brief introduction to the theory of logarithmic wave equations, assuming a single-channel case for simplicity. In section IV we consider multi-channel systems, introduce the matrix generalization of the logarithmic wave equation, and formulate the basic notions of the approach.

II Logarithmic wave equation: Original formulation

Let us consider a material (fluid being a special case), which is in thermal contact with a reservoir of infinitely large heat capacity so as to maintain constant temperature. We assume that its microscopic structure can be regarded as a many-body system of particles, atoms and molecules, whose interaction energy is larger than its kinetic one; we shall refer to such materials as condensate-like. This class includes not only low-temperature strongly coupled systems, but also any matter in which density or geometrical constraints make interparticle interaction potential energies dominate over kinetic energies.

From the viewpoint of statistical mechanics, in such materials the microstates are naturally represented by a canonical ensemble; therefore their probability $P$ is given by a standard formula, in which kinetic energy is neglected compared to potential energy: $P\propto\exp{(-{\cal E}/T)}\approx\exp{(-{\rm U}/T)},$ where $T$, ${\cal E}$ and ${\rm U}$ are the absolute temperature, energy and potential energy of a many-body system, respectively; we shall work in the units where the Boltzmann constant $k_{B}=1$.

Let us also assume that the state of such a material can be described by a single complex-valued function written in a Madelung form (Rylov, 1999):

[TABLE]

where $\rho=\rho(\textbf{x},t)$ is a fluid density, and $S=S(\textbf{x},t)$ is a phase which is related to fluid velocity: $\textbf{u}\propto\boldsymbol{\nabla}S$. This function naturally obeys a normalization condition

[TABLE]

where $M$ and $V$ are the total mass and volume of the material. This poses restrictions which are similar to a quantum mechanical case: the set of all normalizable functions must constitute a Hilbert space, such as $L^{2}$. Assuming that the dynamics in such space is of the Hamiltonian type, one can derive the operator of the potential $\hat{\rm U}\sim T\ln{(A|\Psi|^{2})}=-K(T-T_{0})\ln{(\rho/\rho_{0})}$, where $A=1/\rho_{0}$ and $K$ are proportionality constants, and $\rho_{0}$ and $T_{0}$ are some critical values of, respectively, density and temperature at which the statistical effect vanishes (Zloshchastiev, 2018a).

Then the energy conservation law can be recast in the form of a logarithmic Schrödinger-like equation:

[TABLE]

where $b=K(T-T_{0})/\eta$, $\eta$ and ${\tilde{D}}$ are positive constants. Here, $\eta$ has a dimensionality of a Planck constant but not necessarily its value; from now on, we assume that $K$ is absorbed into $\eta$ for brevity. The parameter ${\tilde{D}}$ can be also interpreted as ${\tilde{D}}=\eta/m_{\text{eff}}$ where $m_{\text{eff}}$ is a mass of a particle.

It follows also that the nonlinear coupling must be linearly related to the thermal temperature $T$ and to the wave-mechanical temperature $T_{\Psi}$

[TABLE]

where the latter is defined as a thermodynamical conjugate of the entropy function $S_{\rho}=-\int_{V}\rho\ln{\!(\rho/\rho_{0})}\,dV$. Entropy functions of this kind directly emerge from the logarithmic term when one averages (II.3) using the inner product in a Hilbert space of functions $\Psi$ (Brasher, 1991; Zloshchastiev, 2010).

III Fluid-Schrödinger analogy

In this section, we derive a formal mapping between wave equations of type (II.3) and flow equations for a certain class of materials.

III.1 General case

Let us begin with a scalar wavefunction case (i.e., $\Psi$ is a matrix of rank one), and consider a general $U(1)$-symmetric Schrödinger equation of the form

[TABLE]

where $F(\rho)$ is a differentiable function on a positive semi-axis $\rho$, and $\Psi$ obeys the normalization condition (II.2).

By substituting the Madelung ansatz (II.1) into Eq. (III.1), and separating real and imaginary parts, one obtains

[TABLE]

where a dot means a scalar product. Taking a gradient of the latter equation and introducing the fluid velocity

[TABLE]

we obtain hydrodynamic laws for mass and momentum conservation for a two-phase compressible inviscid fluid with internal capillarity whose flow is irrotational and isothermal:

[TABLE]

with the stress tensor $T$ being of the Korteweg form with capillary effects (Dunn and Serrin, 1985):

[TABLE]

where $I$ is the identity matrix, and $\tilde{p}$ is a capillary pressure:

[TABLE]

where

[TABLE]

is a barotropic equation of state for our fluid, where prime means a derivative with respect to the value in brackets. Furthermore, the speed of sound can be derived as

[TABLE]

where the second term is non-zero only if $\boldsymbol{\nabla}^{2}\rho$ can be presented as a function of density.

Note that in the leading-order approximation with respect to ${\tilde{D}}$, one can neglect the capillary term and obtain

[TABLE]

In other words, Eq. (III.1) is a concise form of writing two fluid-mechanical Eqs. (III.8). This makes Schrödinger-like equations useful for studies of the Korteweg-type materials, as well as for certain generalizations to be discussed below.

III.2 Logarithmic fluid

Let us apply the formulae above for a special case of a logarithmic nonlinearity. We therefore impose

[TABLE]

which corresponds to Eq. (II.3). Then formulae (III.8)-(III.10) remain intact, but Eq. (III.11) simplifies to a linear equation:

[TABLE]

which belongs to a perfect-fluid class.

Finally, in the leading-order approximation with respect to ${\tilde{D}}$, one obtains

[TABLE]

which means that the logarithmic fluid is perfect, and its speed of sound is independent of density (Zloshchastiev, 2011) These two properties are useful for applying to models of those strongly-interacting fluids for which the perfect-fluid approximation can serve as a robust, or leading-order approximation.

IV Multi-channel systems and matrix logarithmic wave equation

Historically, studies of multi-channel phenomena started from multi-channel scattering in particle physics. For instance, in the reaction between two nuclei such as $~{}^{12}$C + d, a variety of mass rearrangements may be possible, such as $~{}^{13}$C + p, $~{}^{14}$N + $\gamma$, called the mass partition, but in all cases the sum of the particle masses in each partition is almost the same for each of the partitions. Therefore, one needs to describe the whole process, taking into account not only the existence of multiple channels but also their possible interactions (couplings) with each other.

In a broader context of condensed matter physics, multi-channel phenomena can be considered through the coupled-cluster method (CCM) as described in the work by Čížek (1966). According to this method, the search for a wave function is performed by means of the pre-multiplication ansatz, sometimes called also a “preconditioning” of a wavefunction (Acton, 1997):

[TABLE]

where $\Psi_{0}$ is a simple and often time-independent Slater determinant reference function, and $\hat{{\cal M}}$ is an operator written in an exponential form: $\hat{{\cal M}}=\exp\hat{\cal S}$. In quantum chemical problems, the fast convergence of the results was reported (McClain et al., 2016; Hagen et al., 2016).

In essence, the coupled-cluster method is about reconstructing the operator $\hat{\cal S}$ for a specific physical state $\Psi$. In other words, one replaces constructing of a Hilbert-space state vector by constructing of a Hilbert-space operator. The latter is usually an $N\times N$ matrix $\hat{{\cal M}}$, $N\to\infty$, but generally one also applies the truncations, i.e., finite-dimensional approximations $M<\infty$.

For our purposes, one can assume that we deal with multiple flows of the Korteweg-type with similar properties. It implies that in a free logarithmic wave equation (II.3) one should upgrade the function $\Psi(\textbf{x})=\langle\textbf{x}|{{\cal M}}\rangle$ from a scalar function to an $N\times N$ matrix function ${{\cal M}}$ (Znojil et al., 2017):

[TABLE]

whereas the normalization condition (II.2) is replaced with

[TABLE]

where a symbol “tr” denotes a conventional matrix trace. Consequently, the total mass density of such a system can be defined as

[TABLE]

then the normalization condition (IV.3) takes a habitual form.

Furthermore, the Madelung ansatz can be chosen in the form

[TABLE]

where ${\cal R}={\cal R}(\textbf{x},t)$ and ${\cal S}={\cal S}(\textbf{x},t)$ are self-adjoint invertible square matrices of rank $N$. Adjoint and invertible matrix conditions are imposed here to satisfy the correspondence principle: we must recover a conventional nonlinear Schrödinger equation at $N=1$. Note that ${\cal R}$ and ${\cal S}$ do not have to commute in general, and their order in Eq. (IV.5) is chosen to eliminate ${\cal S}$ in the logarithmic term in (IV.2).

Furthermore, matrices ${\cal R}$ and ${\cal S}$ are connected to the observables via the formulae

[TABLE]

for which reason we call $\boldsymbol{\cal U}={\tilde{D}}\boldsymbol{\nabla}{\cal S}$ the velocity vector matrix, and ${\cal R}^{2}$ the mass density matrix.

Substituting the ansatz (IV.5) into Eq. (IV.2), we obtain

[TABLE]

where ${\cal R}^{-1}$ denotes the inverse matrix of ${\cal R}$. From these one can derive the following equations

[TABLE]

where ${\cal R}^{-2}$ denotes the inverse of the mass density matrix. Here it is noticeable that the dimensional parameter $\rho_{0}$ has disappeared from the equations, indicating that the mass density and flow velocity of our material do not depend on a scaling constant $\rho_{0}$.

Applying the matrix trace operation to these equations and recalling the definitions (IV.6), we obtain flow equations for the observables

[TABLE]

where

[TABLE]

is a cumulative acceleration acting on the system, where

[TABLE]

is an acceleration vector matrix.

Equations (IV.16) can, of course, be recognized as matrix analogues of the continuity and Cauchy momentum equations for our system, respectively. These equations can be difficult to use for actual derivation of the observables $\rho$ and u; due to their complex structure involving products of non-commuting matrices. However, one should bear in mind that they are equivalent to one matrix equation (IV.2), which can be solved with considerably more ease. Once its solution is found, equations (IV.16) can be used for a fluid-mechanical interpretation of results.

V Examples

As an example, in this section we consider solutions for a case when ${\cal M}$ is a $2\times 2$ matrix, for reasons of simplification. For the same reasons, we also limit ourselves to a single spatial dimension: $\textbf{x}\to x$. Then the wave equation (IV.2) becomes

[TABLE]

while the normalization condition (IV.3) becomes

[TABLE]

assuming that the system is localized in the interval of $x\in[x_{L},\,x_{R}]$; this formula comes from the normalization condition $\int\limits d\vec{x}\,{\rm tr}({\cal M}^{\dagger}{\cal M})\equiv\langle{\cal M}|{\cal M}\rangle=M,$ where integration is taken over a spatial volume occupied by the system, $M$ being a constant usually interpreted as a number of particles inside such a volume.

Furthermore, below we consider the cases when the exact analytical solutions of Eq. (V.1) are known.

V.1 Diagonal case

Let us assume the simplest case when ${\cal M}$ is imposed to be a diagonal $2\times 2$ matrix. The imposed diagonality indicates that the cross-channel coupling terms can be neglected, i.e., our configuration or material is a mixture of two independent components. The ground-state solution of Eq. (V.1) can be found analytically. After some algebra, we obtain

[TABLE]

where $\psi_{a}({x},t)$ are complex-valued functions. In the fluid rest frame of reference, the solution has the form of the gausson, i.e., the Gaussian parcel modulated by de Broglie plane wave,

[TABLE]

where $a=1,\,2$, $\gamma=2b/{\tilde{D}}>0$, and

[TABLE]

is energy of a wave in the $a$th channel. The integration constants $C_{a}$ and $\nu_{a}$ are related, together with $b$, to a mean and variance of the solution. Note that this solution exists for models with a positive value of $\gamma$, as discussed in the concluding section.

If one imposes also a normalization condition (V.2) here, then one obtains an additional constraint for the integration constants

[TABLE]

where we denoted $\Delta F_{a}=F_{a}(x_{R})-F_{a}(x_{L})$ and $F_{a}(x)=C_{a}^{2}\,\exp{\left(\frac{\nu_{a}^{2}}{\gamma}\right)}\,\text{erf}\left(\frac{\nu_{a}-\gamma x}{\sqrt{\gamma}}\right)$. For instance, when $x_{R}$ and $x_{L}$ are set to, respectively, plus and minus infinity, then this constraint takes a simple form $\sum\limits_{a=1}^{2}C_{a}^{2}\,\exp{\left(\frac{\nu_{a}^{2}}{\gamma}\right)}=\sqrt{\frac{\gamma}{\pi}}\frac{M}{\rho_{0}}$.

Note that due to the Galilean symmetry of a logarithmic Schrödinger equation, one can always generate solitary wave solutions from a solution (V.4), whose centers of mass propagate with velocity $v_{a}$, independently for each channel. One can check also that the solution (V.3) naturally has a commutative property, $\left[{\cal M}^{\dagger},\,{\cal M}\right]=0$, thus eliminating the ordering uncertainty inside a logarithmic term.

Finally, one can verify that the cumulative acceleration (IV.17) and acceleration matrix (IV.18) both vanish on the solution (V.3), (V.4):

[TABLE]

which means that the material derivative of flow velocity is zero for this solution.

V.2 Off-diagonal case

Now let us assume that the channel-coupling terms are dominating. Therefore, one can impose that ${\cal M}$ is an off-diagonal matrix

[TABLE]

where $\bar{\psi}_{a}({x},t)$ are complex-valued functions.

The ground-state solution of Eq. (V.1) can be found analytically. In a rest frame, it also has the form of the Gaussian modulated by de Broglie plane wave,

[TABLE]

where $a=1,\,2$, and

[TABLE]

is energy of a wave for the $a$th channel, $C_{a}$ and $\nu_{a}$ are integration constants. As in the diagonal case above, this solution exists for models with a positive value of $\gamma$, as discussed in the concluding section.

If one imposes also a normalization condition (V.2) then one obtains an additional constraint for the integration constants, which is identical to Eq. (V.6).

Finally, acceleration matrix (IV.18) takes a non-zero value on the solution (V.8), (V.9), given by:

[TABLE]

but the cumulative acceleration (IV.17) still vanishes, thus indicating that the material derivative of flow velocity is zero for this solution.

VI Conclusion

In this paper, we formulated a formal mapping between a large class of nonlinear wave equations and flow equations of a barotropic fluid with internal surface tension and capillary effects. Then we generalized quantum wave equations of a logarithmic type to a matrix form. Furthermore, we mapped this equation to matrix flow equations of multi-channel or multi-component Korteweg-type materials including fluids with internal surface tension and capillary effects.

Aside from being a concise way of writing flow equations, such mapping reveals the existence of a Hilbert space of states represented by different solutions, for the same set of boundary conditions. This means that the material can randomly “choose” either of these states, so that a statistical ensemble occurs.

Furthermore, in Section V, we studied two simplest examples of the model, namely, a case of $2\times 2$ matrices in one-dimensional space, assuming a positive value of $\gamma$ and either a diagonal or off-diagonal form. We analytically derived Gaussian-type matrix solutions, and studied them in the context of fluid mechanics. It should be emphasized that the solutions presented above exist only for models with positive values of $\gamma$ where they usually represent states with lowest energy - an analogue of ground states in quantum mechanics. For negative values of $\gamma$, the topology of the problem changes crucially, and the topological solitons emerge, as was shown for the scalar (rank-1 matrix) case (Zloshchastiev, 2018b).

We expect that these matrix fluid models will find application not only in the theory of low-temperature superfluids, but also for mechanical studies of arrays of microfluids and other multi-channel flows where capillary effects are substantial.

Acknowledgements.

I am grateful to participants of the XXIII Fluid Mechanics Conference KKMP-2018 in Zawiercie, Poland (9-12 September 2018), where preliminaries of this work were reported and discussed. Proofreading of the manuscript by P. Stannard is greatly appreciated as well. This work is based on research supported by the Department of Higher Education and Training of South Africa, as well as by the National Research Foundation of South Africa under Grant No. 95965.