CC BY
0Russian Journal of Nonlinear Dynamics, 2022, vol. 18, no. 3, pp. 441-464. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd220308
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 76Dxx, 35Q35, 65F22
Construction of Inhomogeneous Velocity Fields Using Expansions in Terms of Eigenfunctions of the Laplace Operator
E. V. Vetchanin, E. A. Portnov
In this paper we present a method for constructing inhomogeneous velocity fields of an incompressible fluid using expansions in terms of eigenfunctions of the Laplace operator whose weight coefficients are determined from the problem of minimizing the integral of the squared divergence. A number of examples of constructing the velocity fields of plane-parallel and ax-isymmetric flows are considered. It is shown that the problem of minimizing the integral value of divergence is incorrect and requires regularization. In particular, we apply Tikhonov's regu-larization method. The method proposed in this paper can be used to generate different initial conditions in investigating the nonuniqueness of the solution to the Navier-Stokes equations.
Keywords: inhomogeneous velocity field, expansion in terms of eigenfunctions, ill-conditioned system of linear algebraic equations
1. Introduction
1. One of the special features of the equations governing the dynamics of an incompressible viscous fluid which attracts the attention of researchers is the possibility of emergence of a nonunique solution. We mention some of the studies dealing with the problem of nonuniqueness of a solution. The authors of [5] construct self-similar solutions for the plane-parallel and ax-isymmetric flow in a half-infinite channel with the fluid performing uniformly accelerated motion on one of its boundaries. It is shown numerically that, depending on the Reynolds number, the
Received August 04, 2022 Accepted September 22, 2022
This work was carried out within the framework of the state assignment of the Ministry of Science and Higher Education of Russia (FEWS-2020-0009).
Evgenii V. Vetchanin eugene186@mail.ru Evgenii A. Portnov pereliv00@mail.ru
Ural Mathematical Center, Udmurt State University ul. Universitetskaya 1, Izhevsk, 426034 Russia
equation for the self-similar solution can have several solutions. In [11], a similar problem for the three-dimensional flow is addressed. The paper [15] studies the transition to chaos in the problem of the flow in a planar infinite channel with injection on the boundaries and demonstrates the onset of various solutions as the Reynolds number changes. In [6, 7], the multiplicity of solutions is shown in the problem of the flow between two rotating disks. In [8], the nonunique-ness of the solution is found in calculations of the three-dimensional flow past a wing using the Navier-Stokes equations averaged over the Reynolds number.
A nonunique solution can arise within the framework of the model of an ideal incompressible fluid as well. For example, in the problem of the motion of point vortices on a sphere, two models corresponding to two different solutions of the equations of motion of an ideal fluid on a sphere [4] have been constructed. The first (classical) solution possesses one singular point corresponding to a vortex of strength r, and is given by the following stream function [1, 2]:
^ = -rin(1 - cos 9), 9 e [0, (1.1)
The second solution possesses two singular points corresponding to a vortex and its antipode
with strengths r and —r, and is given by the following stream function [3]:
9
0 = — 2r lntan -, 6 e [0, 7r]. (1.2)
In the first case the flow possesses uniform vorticity r, and in the second case, zero vorticity.
2. Of interest is the question posed by Professor A.V. Borisov concerning two described models of the motion of point vortices on a sphere, namely, which of them is more physical and more realistic? To solve this problem, it was proposed to numerically model the evolution of random initial vorticity distribution on a sphere. In this case, the initial flow can be given by expansion in terms of spherical harmonics (eigenfunctions of the Laplace operator in spherical coordinates) for vorticity1:
N n
mi
U = Yj Yj COS mV + Bn,m sin m^P^ (COS d), (1.3)
n=1 m=0
where Pm(cos 9) is the associated Legendre function of the first kind, N is the number of spherical harmonics involved in the construction of the vorticity field, and the expansion coefficients An m and Bn m can be chosen in an arbitrary way.
The known vorticity distribution allows us to find the stream function ^ from the equation
1 d2ip 1 d . dip R2 sin2 0 dip'2 + R2 sin 0 ~d9 Sm
The solution to Eq. (1.4) for a sphere of radius R and the vorticity distribution (1.3) will have the form
Nn
9) = YY (A'n,m COS mV + B'n ,m sin m&Pn? (cOS 9),
n=l m=0 (15)
A R2 B R2
t/ __p/ _ JJ,n,mri'
-->--) nryrt / . \ i irt m
n,m
n(n + 1)' nm n(n + 1)
1 This approach was proposed by Professor S. M.Ramodanov.
By varying the number of harmonics N and the coefficients A'n,m, B'n m in the expansion (1.5), one can obtain a great variety of inhomogeneous flows on the sphere with different scales of inhomogeneities (see Fig. 1). Further calculation of the evolution of vortex structures of such a flow can be performed using some numerical algorithm. However, the question whether specific models of the motion of vortices on a sphere are realistic still remains an open problem.
(a) N = 2 (b) iV = 4
Fig. 1. Examples of the flow on a sphere, defined by the expansion (1.5), for different N
3. The approach to constructing the vorticity field on a sphere with subsequent explicit recovery of the stream function, as described in the preceding section, turns out to be rather simple to use since there are no solid boundaries on the sphere. In the presence of solid boundaries on which the no-slip condition is given, the vorticity distribution must be kinematically consistent with boundary conditions for the velocity [9]. From a physical point of view, this requirement is a consequence of vorticity generation on solid boundaries. For an incompressible fluid in which there are no sources or sinks, the condition of kinematic consistency has the form
- DH(r, = / ^ - /
D dD
f (n(r') ■ Vh(r'))(r — r') , , x
-J I r-rY- dS> reS> (L6)
dD
where a is the internal angle at the point of the boundary related to 2n(d — 1), d = 2, 3 is the dimensionality of the problem, u is the vorticity field, Vh is the boundary velocity distribution, n is the external normal to the boundary of the region, D is the region occupied by the fluid, and &D is the boundary of region D.
Following the ideas described in [9], the construction of an inhomogeneous two-dimensional or three-dimensional flow can be performed in three stages.
At the first stage, the vorticity field satisfying homogeneous boundary conditions of the second kind is given. For example, for the rectangular region the components of such a field can
be written as
N M P
. . V^ V^ tk nnx nmy npz
Uk(x, n,m,p COS — COS -g- COS
n=0 m=0 p=0
x e [0, L], y e [0,H], z e [0,W]. At the second stage, following [9], it is necessary to solve the boundary integral equation
a2n(d - 1) (Vb(r) - 7(r) x n(r)) = f j.* ^ r° dQ'+
\r — r'\d D
r (^-n^x^^xtr-/) _ r (nW- W-rQ
I \r — r'\d J \r — r'\d
dD dD
to determine the strength of the vortex layer 7 on the boundary of the region responsible for vorticity generation (the flow of vorticity on the solid boundary).
At the third stage, the velocity field is calculated using the Biot - Savart law:
V
J \r-r'\d 2-/r(d — 1) J \r-r'\d v ;
2n(d — 1) J \r — r'\d 2n(d — 1) J \r — r'\d
D dD
Here the boundary distribution 7 defines indirectly the corrective additional term to the vorticity field and matches it with the boundary conditions for the velocity.
4. The approach described in the preceding subsection proceeds from some initial vorticity field which is matched with boundary conditions using the solution to the boundary integral equation (1.8). In this paper we propose an approach using as a starting point the representation of the velocity field components in the form of expansions in terms of eigenfunctions of the Laplace operator. The coefficients of the above-mentioned expansions are chosen from the condition of minimization of divergence, and the inhomogeneities in the velocity field are obtained by fixing some of the expansion coefficients. It is shown that the divergence minimization problem reduces to an incorrect system of linear algebraic equations whose solution can be found using Tikhonov's regularization method2 [13, 14].
One of the potential applications of the proposed method is to numerically investigate the multiplicity of solutions to the Navier-Stokes equations in simple regions admitting the construction of expansions in terms of eigenfunctions of the Laplace operator.
2. The main relations of the method for constructing an inhomogeneous velocity field
The method presented in this paper for constructing the inhomogeneous velocity field of an incompressible fluid is based on the following identity [10] known from vector analysis:
Vx(VxV ) = V(V-V) — AV, (2.1)
where V denotes the velocity field of the fluid.
2In a number of publications, the term "Tikhonov - Phillips regularization" is used since a similar regularization method was proposed in [12] independently of A. N. Tikhonov's results.
Given that the fluid is incompressible (V-V = 0), we write relation (2.1) as
AV = -V x u, (2.2)
where u = Vx V is the vorticity field.
In the case of known vorticity distribution u (kinematically consistent with boundary conditions for velocity) the velocity field V can be found using an analytical or numerical solution of the Poisson vector equation (2.2). If in the region occupied by the fluid it is possible to introduce an orthogonal curvilinear coordinate system q1, q2, q3 in such a way that the boundaries of the region coincide with its coordinate surfaces (qk = const, k = 1, 2, 3), then the solution to Eq. (2.2) can be constructed in the form of expansions in terms of the eigenfunctions of the Laplace operator:
Vk^ Q2, ?3) = Q2, ?3) + An,m,p^n,m,p(q^ Q2, Q3), k = 1 2 3 (2.3)
n=l m=1 p=l
where $k is a function (usually represented as a series) satisfying inhomogeneous boundary conditions for the velocity component Vk, is the eigenfunction of the kth component of
the Laplace vector operator, and An,mp are the expansion coefficients to be defined.
In the case where the vorticity field is unknown, the inhomogeneous velocity field can be obtained by choosing the coefficients An,mp in a fairly arbitrary way. However, such a velocity field does generally not satisfy the incompressibility equation
V-V = 0. (2.4)
Therefore, to find the coefficients An,mp, we suggest solving the problem of minimizing the integral of the squared divergence of velocity:
F (An, m ,p, An, m, p, An, m,p) = \VV |2|JI dq1 dq2 dq3 ^ ^^ (2.5)
D
where J is the Jacobi matrix of transformation to the curvilinear coordinate system.
Remark 1. In a sense, the problem of minimizing the integral of the squared divergence (2.5) is similar to the well-known method of least squares.
Since the coefficients An,mp appear in the expansions (2.3) linearly and the divergence operator is also linear, the problem (2.5) of finding the coefficients An,m,p reduces to minimizing the quadratic function:
F (An, m,p, An, m,p, Aim mp) = (X, QX) + (b, X) + C ^ min, (2.6)
where X is the vector composed of the coefficients An,m p, Q is a positive definite matrix whose form depends on the form of the eigenfunctions ^ m p, b is the vector whose form depends on
the form of the eigenfunctions ^kn, m p and the functions $k, and C is the constant depending on the form of the functions $k. Thus, to find the coefficients An , it is necessary to solve the following system of linear algebraic equations:
QX = -l-b. (2.7)
We point out some details related to the solution of the problem (2.6):
1. For practical calculations we will use, instead of complete expansions (2.3), their partial sums.
2. In the examples given below, we show that, even if the number of terms in the partial sums corresponding to the expansions is relatively small, the matrix Q possesses condition numbers of order 1016. Thus, the minimization problem (2.5) can become incorrect, and its solution requires regularization.
3. To obtain inhomogeneities in the velocity field, we will fix the values of some of the coefficients m, p.
Next, we consider examples of constructing inhomogeneous velocity fields in the following cases:
1. Plane-parallel flow in a rectangular region.
2. Axisymmetric flow in a circular tube.
3. Plane-parallel flow in a rectangular region
For the plane-parallel flow in a rectangular region the vector equation (2.2) and the minimization problem (2.5) take the form
dx2 dy2 dy' dx2 dy2 dx' F = + (32)
D
The form of the expansions defining the solution to Eqs. (3.1) and the explicit form of the functional (3.2) depend on the boundary conditions for the velocity components.
3.1. Example 1. The velocity profiles are given on all boundaries
Consider the flow in a rectangular region D = {(x, y) | x G [0, L], y G [0, H]} on whose boundaries the conditions of the first kind for the components of the velocity vector V = (Vx, Vy) are given:
V(0, y) = (A(y), 0), V(L, y) = f2(y), 0), y G [0, H], V(x, 0) = (0, gx(x)), V(x, H) = (0, <72(x)), X G [0, L]. We note that any of the functions fk and gk can be identically zero.
The partial sums that describe approximately the solutions to Eqs. (3.1) with the boundary conditions (3.3) can be written as
TTf x K1 ak sinh(h(L - x)) + fk sinh(Xkx) . , ,
vx(x, y) = Y——v ■ J: ^ -—sin{\ky +
k=1 sinh( xkL)
, ^ ra,m Vlh '
n=1 m=1 v
V) = ± 'J" M^H
' t=! slnH HH)
n=1 m=1 v
where
H H
2
ak = jj J /i(y)sin(Afcy)dy, Pk = jj J ¡-¿iv) sin(Afcy) dy, 0 0
L L . .
2 f 2 f (3.6)
lk = L / 9i(x)sm(fj,kx)dx, Sk = - / g2(x) sin(^kx) dx,
00 nk nk
IF' Atfc = T
Remark 2. When one calculates terms with hyperbolic functions in the expressions (3.4) and (3.5), machine infinities Inf can simultaneously arise in the numerator and the denominator of these expressions. The operation of dividing machine infinities leads to the onset of values NaN (Not a Number). In order to avoid the appearance of machine infinity in the numerator, it is convenient to represent the terms with hyperbolic functions in the expressions (3.4) and (3.5) in the form
1 I ak tanh(Afc(L1 - x)) ¡3k tanh(Afcx) ^
Afc — "77) — ~F) fc — 1, 2, 3, ...
tanh(Ak(L1 — x)) + tanh(Akx) \ cosh(Akx) cosh(Ak (L1 — x)) I
With the terms written in this way, the division of the finite quantity by machine infinity Inf will be performed. The result of this operation will be 0.
Using the expressions (3.4) and (3.5), we calculate the quadratic, linear and constant parts of the functional (3.2), respectively:
(X, QX) =
Nx N Nx N Nx N
v M L__
7r2 (a2 — n2)(m2 — p'2)
m
Nx N Nx N Nx N
7T2 (a2 — n2)(m2 — p2)
(3.7)
VVU' \ V V V A' B> ((-!)»+»-l)((-l)"^-l)
n=1 m=1 n=1 m=1 a=1, a=n ¡3=1, fi=m
+ V sr E> A> 4n[] ((-l)"+"-l)((-l)m+/3-l)
Z^ Z^n,m\ + Z^ Z^ Z^ Z^ nn,mn-a,fi _ „2^
n=1 m=1 n=1 m=1 a=1,a=nf=1,f=m
Nx min(K1,Ny)
(b, X) = £ E An
2 A:
n=1 m=1
IH
n,m /,2 I \2 V T + Am ' L
—am + (—1)n PJ +
min(K2,Nx) Ny
+ E E^
2/n
L
n=1 m
~n,m / 2 + a^ V H
= 1 Fn, + Am ' H
-Yn + (—1)m^J +
K1 Nx Ny
+E E E ^n,m
k=1 n=1 m=1, m=k
4fcM„Afc [H (-l)fc+m-l cosh(AfcL) - (-1)
n(p?n + A2k) V L k2 — m2
sin
h( Ak L
(ak + (—1)n fa) +
K Nx Ny
+ V" 4fcAtfcAm [L (-j)fc+ra -1 cosh(^feg) - (~ir , ,
+ ¿^ An,mzTvrmrV 77 —¿,2 _ „2----w*- + (-1) °k)>
k=1 n=1, n=k m=1
n(Am + /k) V H k2 — n2
sin
h(H
(3.8)
n
K
1 X2 H
c = Y(4 + f3ï) k k=l
L
+
coth (Xk L)
4 Uinh2 {XkL Xk
X,,
Ki
X2k H
k=l
2 sinh (XkL)
L coth (XkL) +
X
+
H
4 y sinh2 H)
+
coth (l~tkH)
H
K
Yk^k n •
k=l
K1 K2 \ 2p2X2
+ E E ("/■ -/ + (-1 )k«kh + (-1)^7/ + C D" ? vv)
k=l / =1
2
/ + Xk)
(cosh(AfcL) - (-l)-Q (costing) - (-1)* sinh(Xk L) sinh(H)
(3.9)
where An,m = An,m^n, B'n m = Bn,mXm, X is the vector composed of the coefficients
, m, Q is the matrix of the quadratic form (3.7), and b is the vector composed of the coefficients appearing in the expression (3.8).
Remark 3. The scaling of the expansion coefficients A'n m = An and B'n m = Bn mXm, which was performed in (3.7) and (3.8), allows one to make the matrix Q independent of the dimensions of the computational region.
Analysis of the expression (3.7) shows that the system of linear algebraic equations to which the minimization problem (3.2) leads in the case at hand breaks down into four independent subsystems of the form
Qkxk —
k = 1, 2, 3, 4,
(3.10)
where Ej1 and E| are identity matrices, Sk is in the general case a rectilinear matrix, and the vectors Xk are composed as follows:
Xl = (Al,l, Al,3, Al,5,
X2 = (Al,2, Al,4, Al,6,
X3 = (A2,l, A2,3, a2,5 ,
X4 = (A2,2, A2,4, A2,6,
A' A' A'
A3,l, A3,3, a3,5,
A' A' A'
a3,2, a3,4, a3,6,
A' A' A'
A4,l, a4,3, a4,5,
A' A' A'
A4,2, A4,4, A4,6,
B2,2, B2,4, B2,6,
B2,l, B2,3, B2,5,
Bl ,2, Bl ,4, Bl ,6,
Bl,l, Bl,3, Bl,5,
, B4 2, ,4, B4,6, ■ ■■), , B4, 1, B4,3, B4,5, ■ ■ ■), , B3,2, B3,4, B3,6, ■ ■ ■ ),
, B3, 1, B3,3, B3,5, ■ ■ ■ )■
^ 11)
As the number of harmonics Nx, Ny, K1 and K2 in the expressions (3 ■ 4) and (3 ■ 5) increases, the spectral properties of the matrix Qk deteriorate. In particular, when Nx = Ny = K1 = K2 = = N, their largest absolute eigenvalues |Amax(Qk)| tend to 2 as N increases, but the smallest absolute eigenvalues |Amin(Qk)| decrease fast, see Figs. 2a, 2b.
Remark 4. According to Fig. 2b, it may be assumed that Amin ~ exp(— pN), where p = const. The emergence of Amin ~ 10"16 with large N in numerical calculations is evidently due to limitations of the machine representation of real numbers.
1
x
x
PPPPPPPPPPP
—©—Qil
cond(Qfe)
10
-v-Q
(a)
20
JTV10-2O
Fig. 2. (a) The largest absolute eigenvalues, (b) the smallest absolute eigenvalues, (c) the condition numbers of the matrices Qk depending on the number N of harmonics retained in the expansions (3^4) and (3^5). It is assumed that Nx = Ny = K1 = K2 = N
The decrease in the value |Amin(Qk)| leads to an increase in the condition numbers
cond(Qfc) = '^^J (3.12)
|Amin(Qk )|
of the matrices Qk, see Fig. 2c.
Thus, we see that the systems of equations (3-10) become incorrect as the number of harmonics increases, and their solution requires regularization. Following the work of Tikhonov [13, 14], we will seek an approximate solution to the subsystems (3T0) by minimizing the smoothing functional:
1
Fa =
Q/.*;; + r2h
+ a \\X%||2 ^ min, (3.13)
where a > 0 is the regularization parameter and \\ • \\ is the Euclidean norm defined in a standard way.
The minimization problem (3.13) reduces to solving a system of linear algebraic equations of the form
(Q%Qk + <rE)XZ = ~QTbk, (3.14)
whose condition number has the following estimate for Nx, Ny, Kl, K2 ^ 3:
cond (Qfc Qfc + o-E) w (3.15)
To show how the proposed method works, we fix the dimensions of the region L = 5, H = 1, define explicit expressions for the boundary velocity profiles
/i(y) = i-(f-i)16, = -1)2)' 32^ = 0 (3.16)
and fix the following coefficients of the partial sums (3.4) and (3.5):
A34 = 0.5, B34 = -0.2. (3.17)
2y H
- 1
y e [0, H]
Remark 5. A function of the form
fn(y) = 1 -
with large n allows one to define a velocity profile close to a uniform one, see Fig. 3.
V
1 0.8 0.6 0.4 0.2 0
" -J______ „ _ Ifi ■
n = 4
n = 8
n = 2\|
fn(y)
0 0.2 0.4 0.6 0.8 1
Fig. 3. Examples of velocity profiles given by the function (3.18)
F(Xa)
10
-2
10
-3
Nz = 100, Ny = Ki = 10
o o o o o o o
Nx = 250, Ny = Kx = 20
Nx = 250, Ny = Ki = 40
□ □ □ □ □ □ D
v V V V V
10
M
Nx = 100, Ny — K\ — 10 Nx = 250, Ny = Ki = 20
3•10-3
10-15
Fig. 4. (a) The value of the functional (3.2) and (b) the value of the norm of the residual of the system (2.7), which have been calculated from the solution to the regularized system (3.14) for different values of a
Figure 4 shows the dependences, on the parameter a, of the integral value of divergence, i.e., the functional F(Xa), and the value of the residual R(Xa) = ||QXCT + , which have been calculated from the regularized solution Xa, for different values of Nx, Ny, K1. It can be seen from Fig. 4b that, as the regularization parameter increases, so does the value of the residual R(Xa). However, an essential growth of the integral divergence occurs only when a ~ 10_6, see Fig. 4a. Thus, the solution to the regularized problem can ensure small values of the integral divergence.
The vorticity field and the divergence distribution for the flow satisfying the boundary conditions (3^3), which have been constructed for a = 10_6 and for different values of Nx, Ny and K1, are shown in Figs. 5 and 6. A comparison of Figs. 5b and 5e shows that the increase in the number of longitudinal harmonics Nx leads to a decrease in the level of divergence whose largest values arise at the entry to the computational region (see Figs. 5c and 5f). In this case, the vorticity field and the form of the streamlines do not qualitatively change (see Figs. 5a and 5c). Figure 5a corresponds to the value of the functional (3^2) F ~ 4• 10_3, and Figure 5d corresponds to the value of the functional (3^2) F « 10"
^-3
Nx = 100, Ny = K1= 20, K2 = 0
Nx = 250, Ny = K1 = 20, K2=0
0.04 0.06
x
0.04 0.06
x
(c)
(f)
Fig. 5. (a) and (d) Vorticity field, (b) and (e) divergence distribution, (c) and (f) the initial area of divergence distribution. Regularization parameter a = 10_6
As the number of transversal harmonics Ny and Kl increases, so does the level of divergence at the entry to the computational region (compare Fig. 6c to Fig. 6f). However, the integral value of divergence becomes smaller than that in the preceding cases: F « 4.5 • 10"4. It should be noted that the flow pattern does not qualitatively change (see Fig. 6a), but oscillations arise in the vorticity field in the direction of the axis Y. Suppression of such oscillations remains an open problem.
3.2. Example 2. Conditions of the second kind are given on one of the boundaries
Consider the flow in a rectangular region D = {(x, y) | x e [0, L], y e [0, H]} on whose four boundaries the following conditions are given: on three boundaries, the conditions of the first kind for the components of the velocity vector V = (VX, Vy), and on the fourth boundary, homogeneous conditions of the second kind:
V(0, y) = f (y), 0), V(x, 0) = V(x, H) = 0, 3V (L,y)
dx
= 0.
(3.19)
This version of the problem is interesting for constructing the velocity field in a planar channel with the pressure given on one of its boundaries.
Nx = 250, Ny = K\ = 40, #2 = 0
x
(c)
Fig. 6. (a) Vorticity field, (b) divergence distribution, (c) the initial area of divergence distribution. Regularization parameter a = 10~6
The partial sums that describe approximately the solutions to Eqs. (3.1) with the boundary conditions (3.19) can be written as
,) = f;afe^h(Afe(L-,))sin(Afey) + g £ 2sin(^n(Am^ cosh(AfcL) ^„^ VLH
^ Y^V^h 2sin(^ra.T) sin(Amy)
k=1 \ K J n=1 m=1
>nX
! ! y/LH
n=1 m=1 v
where
H
2 f . \7 , nk n(2k - 1) .
= jj / f(v)sm(\ky)dy, h = »k =-2I---^ ^
0
Remark 6. For machine calculations it is convenient to represent the term with hyperbolic functions in the expression (3.20) in the form
ak sin(Xkv)
cosh(Xkx) ( 1 + tanh(Xk(L — x)) tanh(Xkx)
Let us retain in the expansions (3.20) and (3.21) the finite number of terms and calculate the quadratic, linear and constant parts of the functional (3.2):
(X, QX) =
Nx Ny Nx Ny Nx Ny
2
8m(l - (-l)m+/3)
n=1 m=1 n=1 m=1 a=1 ¡3=1, 3=m
E E K,m]2 + EEE E KmK/s^p a - 1 + (-l)"+"(2n - 1)) (m2 - p2) +
NX Ny N Ny Nx Ny ^ _ (_1)m+^
+ E E fin,™]2 + E E E E B'n,mA'a<^2 ,2n _ 1 {_1)a+n{2a _ _ m2\ '
n=1 m=1 n=1 m=1 a=1 p=1, p=m v v/v //v /
(3.23)
Nx min(A',Wy) , v
= E E
n=1 m=1 V ^n + L J
, VV V B' 4Ag(Afc(-l)"+1-/iwsinh(AfcL))(l-(-l)^) ¿¿Î^t^fc ^ (Ag + /4)(Ag-A^)V^coBh(AfcL) ' j
'' E'^,"-!^,) (3.25)
where An,m = An m^n, B'n m = Bn mXm, X is the vector composed of the coefficients An,m, B'nm, Q is the matrix of the quadratic form (3.23), and b is the vector composed of the coefficients appearing in the expression (3.24).
In the expressions (3.23) and (3.24) the scaling of the coefficients An,m = An, m^n and B'n,m = = Bn mXm has been performed. As in the preceding example, this makes the matrix Q independent of the dimensions of the computational region.
Analysis of the expression (3.23) shows that the system of equations (2.7) breaks down into two independent subsystems of the form
1. „ № Sk
QkXk = Qa = Li E2 ) ' k = 2> (3-26)
Jk ^k , ^ Ek/
where Ek and Ek are identity matrices, Sk is in the general case a rectangular matrix, and the vectors Xk are composed as follows:
X1 = (A1 ,1, A1,3, A1 , 5, • • • , A2 ,1, A2 ,3, A2 ,5, • • • , B1,2, B1 , 4, B1 , 6, • B2 , 2, B2 , 4, B2 ,6, • • • , X2 = (A1 , 2, A1 , 4, A1 , 6, • • • , A2,2, A2,4, A2,6, • • • , B1,1, B1 ,3, B1,5, • • • , B2,1, B2 ,3, B2 ,5, • • •) •
(3 • 27)
As in the preceding example, we find that, as the number of harmonics in the partial sums (3 • 20) and (3 • 21) increases, the largest absolute eigenvalue of the matrices Qk tends to 2, and the smallest absolute eigenvalue decreases, see Figs. 7a, 7b. Accordingly, the condition numbers of the matrices Qk increase, see Fig. 7c.
Thus, to solve the subproblems (3 • 26), it is also necessary to perform regularization. Instead of the subsystems (3 • 26) we will solve subsystems of the form (3 -14).
To show how the proposed method works, we fix the dimensions of the region L = 5, H = 1, define an explicit expression for the velocity profile in the inlet cross-section of the channel:
/(y) = l"(f-l)16 (3-28)
and fix the following coefficients of the partial sums (3 • 4) and (3 • 5):
A34 = 0 • 5, B34 = -0 • 2 • (3 • 29)
\K
1.8
1.6
1.4
1.2
0
>00000000000
cond(Qfe)
o Q,
q2
10 20
(a)
10 20
(c)
Fig. 7. (a) The largest absolute eigenvalues, (b) the smallest absolute eigenvalues, (c) the condition numbers of the matrices Qk depending on the number N of harmonics in the expressions (3.20) and (3.21).
It is assumed that Nx = Ny = K = N
In the case dealt with in this subsection, we find that, as the value of the regularization parameter a increases, so do the integral value of divergence F(Xa) and the value of the residual of the system R(Xa) = ||QXCT + ^H^, which have been calculated from the solution to the regularized problem (see Fig. 8). The most significant increase in the integral value of divergence
arises for a ~ 10
-6
r» ~Va _i_ Il.ll
1 2"lloo
= 100, Ny = K = 10
o o o o o o o
N± =. .250, Ny = K = 20
Fig. 8. (a) The value of the functional (3.2) and (b) the value of the norm of the residual of the system (2.7), which have been calculated from the solution of the regularized system (3.14) for different values of a
The vorticity field and the divergence distribution for the flow satisfying the boundary conditions (3.19), which have been constructed for a = 10_6 and for different values of Nx, Ny and K, are shown in Figs. 9 and 10. A comparison of Figs. 9b and 9e shows that the increase in the number of longitudinal harmonics Nx leads to a decrease in the level of divergence whose largest values arise at the entry to the computational region (see Figs. 9c and 9f). In this case, the vorticity field and the form of the streamlines do not qualitatively change (see Figs. 9a and 9c).
Nx = 100, Ny = K = 20
Nx = 250, Ny = if = 20
0 0.02 0.04 x 0.06 0.08 0.1
Fig. 9. (a) and (d) vorticity field, (b) and (e) divergence distribution, (c) and (f) the initial area of divergence distribution. Regularization parameter a = 10_6
Nx = 250, Ny — K — AO
0 0.02 0.04 „. 0.06 0.08 0.1
(c)
Fig. 10. (a) vorticity field, (b) divergence distribution, (c) the initial area of divergence distribution. Regularization parameter a = 10_6
Figure 9a corresponds to the value of the functional (3 • 2) F & 2 • 7 • 10 3, and Fig. 9d corresponds to the value of the functional (3 • 2) F & 8 • 4 • 10"4.
As the number of transversal harmonics Ny and K increases, so does the level of divergence at the entry to the computational region (compare Fig. 10c to Fig. 9f). However, the integral value of divergence becomes smaller than that in the preceding cases: F & 4 • 2 • 10_4. It should be noted that the flow pattern does not qualitatively change (see Fig. 10a), but oscillations arise in the vorticity field in the direction of the axis Y.
4. The axisymmetric flow in a circular channel
The vector equation (2^2) and the minimization problem (2^5) take the form
d2 VZ d2VZ 1 dVZ + -77^ +
dz2 dr2 r dr F =
1 dru d2Vr d2 Vr 1 dVr Vr du
+ -rr^f + ~ = T-, (4.1)
r dr
dz2 dr2 r dr r2 dz
( dVZ dVr Vr \2 , ,
r —+ —--1--) azar —>• mm .
V dz dr r
(4.2)
D
As in the examples considered in Section 3, the form of expansions defining the solution to Eqs. (4T) and the explicit form of the functional (4^2) depend here on the boundary conditions for the velocity components.
4.1. Example 1. The velocity profile is given at the entry and the exit of the region
Let us consider the axisymmetric flow in a circular tube. We denote the computational region D = {(z, r) | z G [0, L], r G [0, R]} and define on its boundaries the following conditions for the velocity components V = (Vz, Vr) and their derivatives:
V(0, r) = (/i (r), 0), V(L, r) = (/2(r), 0),
dVfr' 0) = 0, Vr(z, 0) = 0, V(z, R) = 0.
(4.3)
The partial sums that describe approximately the solutions to Eqs. (4.1) with the boundary conditions (4.3) can be written as
r, - + E E A"
!sin( ßnz)J0( X
k=i
Vr (z, r) = ££ Bn
n=1 m=1
n=1 m=1
a/2 sin(yinz) J1 {ymr)
(4.4)
(4.5)
where
ak
1
R R
zV2 J fi(r)Jo(hr) dr, ßk = jj^2 J f2(r)J0{\kr)dr,
0 k 0
(4.6)
nk
" ■ " \ ^k nk jZ ßk = T' = "77) vk = -j^i -Jk =
L
R
R\Ji(ik )l
jr
) Jk
R\Jo(nk )\
, k = 1, 2, 3, ...
Here £k and nk denote nontrivial zeroes of Bessel functions of orders zero and one, respectively:
Jo (£k )=0, Ji (nk )=0.
Using the expressions (4.4) and (4.5), we calculate the quadratic, linear and constant parts of the functional (4.2), respectively:
(X, QX)
E E [An,m]2 + E E E E ^mB<
, 4aim(—1)m+3 (1 — (—1)n+a)
a-,/3
n=1 m=1
n=1 m=1 a=1, a=n 3=1
K(n2 — a2) (cm — vD
+
, 4nC3 (—1)m+3 (1 — (—1)n+a)
+ E E + E E E E Bra,mAa,/3 ' ( o 2\772 '
i i ,i 1 1 , ,3 M a2 — n2) (42 — nm)
n=1 m=1 n=1 m=1 a=1,a=n3=1 v / v»p <mJ
Nz min(K,Nr)
(b, X) = E E An,
n=1 m=1
K Nz Nr
+ E EE Bn,
k=1 n=1 m=1
V2 J"2
+xm)
(—am + (—1)n m +
nXk (—1)k+m(cosh( Xk L) — (—1)n) JZ
— vm) —n+xk) sinh( xkL)
a+(—1)nPk), (4.8)
K X2
k=1
L
2 \ sinh2 (XkL)
K
Sin
xk
k=1
M xk L
Lcoth{\kL) + Y)[J~}2, (4.9)
where A'n,m = An,m^n, B'n m = Bn,mvm, X is the vector composed of the coefficients An,m and B'n m, Q is the matrix of the quadratic form (4.7), and b is the vector composed of the coefficients appearing in the expression (4.8).
We note that, as in the planar case, after scaling the unknowns A'n,m = An,m^n and B'n m = = Bn,mvm the matrix Q ceases to depend on the dimensions of the computational region.
Analysis of the expression (4.7) shows that the system of linear algebraic equations corresponding to the minimization problem (4.2) breaks down into two independent subsystems of the form
Qkxk — ~2&fc'
k = 1, 2,
(4.10)
where and E| are identity matrices, Sk is in the general case a rectangular matrix, and the vectors Xk are composed as follows:
X1 = (A1,1, A1,2, A1,3,
A3,1, A3,2, A3,3,
X2 = (A2,1, A2,2, A2,3, . .., A4,1, A4,2, A4,3,
B2,1, B2,2, B2,3, B1,1, B1,2, B1,3,
B4,1, B4,2, B4,3, B3,1, B3,2, B3,3,
(4.11)
As in the examples above, we calculate the dependences of the largest absolute eigenvalues |Amax(Qk)|, the smallest absolute eigenvalues |Amin(Qk)| of the matrices Qk and their condition numbers cond(Qk) on the number of harmonics N in the expressions (4.4) and (4.5), see Fig. 11. Here it is assumed that Nz = Nr = K = N. It can be seen from Figs. 11a, 11b that for each of the matrices Qk, k = 1, 2, with increasing N the value |Amax| tends to 2, and the value |Amin| decreases fast, which leads to an increase in the condition numbers of the matrices Qk (see Fig. 7c).
Nz Nr
Nz Nr
Nz Nr
k
k
\K
1.8
1.6
1.4
1.2^
>00000000000
o Q,
q2
10 20
(a)
cond(Qfe)
10 20
(c)
Fig. 11. (a) The largest absolute eigenvalues, (b) the smallest absolute eigenvalues, (c) the condition numbers of the matrices Qk depending on the number N of harmonics in the expressions (4.4) and (4.5). It is assumed that Nz = Nr = K = N
Thus, to solve the subproblems (4.10), it is necessary to perform regularization. Instead of the subsystems (4.10), we will solve subsystems of the form (3.14).
To show how the proposed method works in the axisymmetric case, we fix the dimensions of the region L = 5, R = 1, define an explicit expression for the velocity profile in the inlet and outlet cross-sections of the channel:
/i(r) = /2(r) = l--^ (4.12)
and fix the following coefficient of the partial sums (4.4) and (4.5):
B6>l = 0.2. (4.13)
As in the examples given above, in the case considered in this subsection we find that, as the value of the regularization parameter a increases, so do the integral value of divergence F(Xa) and the value of the residual of the system R(Xa) = ||QX°" + , which have been calculated from the solution to the regularized problem (see Fig. 12). The most significant increase in the integral value of divergence occurs for a ~ 10_s.
The vorticity field and the divergence distribution for the flow satisfying the boundary conditions (4.3), which have been constructed for a = 10_6 and for different values of Nz, Nr and K, are shown in Fig. 13. A comparison of Figs. 13b and 13e shows that the increase in the number of longitudinal and transversal harmonics Nz leads to a decrease in the level of divergence whose largest values arise at the entry to the computational region (see Figs. 13c and 13f). In this case, the vorticity distribution in the flow core changes significantly and the regions with vortices diminish. In particular, it can be seen from Fig. 13f that the largest value of vorticity arises on the wall of the channel. Figure 13a corresponds to the value of the functional (3.2) F & 1.4 • 10_4, and Fig. 13d corresponds to the value of the functional (3.2) F & 6.1 • 10_5.
We note that the vorticity field exhibits pronounced oscillations in the longitudinal direction. The study of the causes of their onset and the ways of suppressing them remains an open problem.
10
-2
F(XC
10
-3
10
-4
10
10-
-5
Nz - 100, Nr — K — 5 1
Nz = 200, Nr = K = io A/
OOQOQO
□ □ u-zr
|Qx<T + èbIL
10
-15
10-io
(a)
10
-5
Fig. 12. (a) The value of the functional (3.2) and (b) the value of the norm of the residual of the system (2.7), which have been calculated from the solution to the regularized system (3.14) for different values of a
Nz = 100, Nr = K = 5
Nz = 200, Nr = K = 10
0.05
0 V 0.5 -0.05
0.05
0 V 0.5 -0.05
0 0.02 0.04 x 0.06 0.08 0.1
(f)
Fig. 13. (a) and (d) vorticity field, (b) and (e) divergence distribution, (c) and (f) the initial area of divergence distribution. Regularization parameter a = 10_6
4.2. Example 2. Conditions of the second kind are given on one of the boundaries
Let us consider the axisymmetric flow in a circular tube. We denote the computational region D = {(z, r) | z £ [0, H], r G [0, and define on its boundaries the following conditions for the components of the velocity vector V = (Vz, Vr):
V(0, r) = f (r), 0), V(z, R) = 0,
dV(L,r) dVz(z, 0) ,
v ' = 0, —' = 0, Vr(z, 0) = 0.
(4.14)
dz
dr
This case is interesting for constructing the velocity field in a channel with the pressure given on one of its boundaries.
The partial sums that describe approximately the solutions to Eq. (4.1) with the boundary conditions (4.14) can be represented as
cosh(Afc(L — z))J0(Xkr) V2sm—nz)J0(\mr)
vz(z, r)=E°* fcV ov *;+E E ° ^r™', ^
k=l COSHAkL) n=l m=l vLJm
VM r) = E E B^^md, (4.16)
where
n=l m=l
R
ak = JJF^ J /(r)^o( V) dr, fj,k = K = = %
k 0
(4.17)
jz _ "I" 1 vsfc/ I jr _ ruv't'/l Z- — 1 O Q ■Jk ~ 7k ! Jk ~ ! h —
R\MZk)\ jr R\MVk)\
V2
Here £k and r/k denote nontrivial zeroes of Bessel functions of orders zero and one, respectively:
Jo (4) = 0, Ji (nk) = 0.
Using the expressions (4.15) and (4.16), we calculate the quadratic, linear and constant parts of the functional (4.2), respectively:
(X, QX) =
Nz Nr Nz Nr Nz Nr 8(m (_1)m+3 + 1
=S £Km]2+S 5 S § - 4) (2«m-1+(-D"+"(2n - D)+
+ E E Km]2 + E E E E KmKß^^2 _ ^ (2n- 1 + (-l)*+"(2a - 1)) ' (4'18)
n=l m=l n=l m=l a=l ß=l
Nz mm(K,Nr) n ^ 2
\ 2 jz
(b,X) = J2 E Km rr 2 +
n=l m=l v
,VVVR 472;/mA2JA:(-l)fc+-(A,rasinh(AfcL) - (—l)"Afc) hkh n,m VLE(Al+A4)(Al-^)cosh(AfcL) ' j
C = E^Manh(AfcL) - (4-20)
where Anm = An,m^n, B'nm = Bn,mvm, X is the vector composed of the coefficients A'nm and B'n m, Q is the matrix of the quadratic form (4.18), and b is the vector composed of the coefficients appearing in the expression (4.19).
In this case, unlike in the examples above, the system of linear algebraic equations corresponding to the minimization problem (4.2) does not break down into several independent
Nz Nr
subsystems. However, this system of linear algebraic equations has a form similar to that which arose in the preceding examples:
1 /E1 S \
QX = --b, Q=(^e2J, (4-21)
where E1 and E2 are identity matrices, S is a square matrix, and the vector X is composed as follows:
X = (A1,1) A1,2, A1,3, • • • , A2,1, A2,2, A2,3, • • • , B1,1, B1,2, B1,3, • • • , B2,1, B2,2, B2,3, • • •)•
(4.22)
Figure 14 shows the dependences of the largest absolute eigenvalue |Amax(Q)|, the smallest absolute eigenvalue |Amin(Q)| of the matrix Q and their condition number cond(Q) on the number of harmonics N in the partial sums (4 • 15) and (4 T6). It is assumed that Nz = Nr = K = N. In this example, the behavior of the spectral properties of the matrix does not qualitatively differ from that considered in the preceding examples.
Fig. 14. (a) The largest absolute eigenvalue, (b) the smallest absolute eigenvalue, (c) the condition number of the matrix Q depending on the number N of harmonics in the expressions (4 T5) and (4 • 16)
To solve the system (4 • 21), we apply the regularization method discussed above.
To show how the proposed method works in the axisymmetric case, we fix the dimensions of the region L = 5, R =1, define an explicit expression for the velocity profile in the inlet cross-section of the channel:
m = 1-^2 (4-23)
and fix the following coefficient of the partial sums (4 T5) and (4 • 16):
B41 = 0 • 1 • (4 • 24)
As in the examples given above, in the case considered in this subsection we find that, as the value of the regularization parameter a increases, so do the integral value of divergence F(Xa) and the value of the residual of the system R(Xa) = ||QX°" + , which have been calculated from the solution to the regularized problem (see Fig. 15). The most significant increase in the integral value of divergence occurs for a ~ 10_s.
The vorticity field and the divergence distribution for the flow satisfying the boundary conditions (4 T4), which have been constructed for a = 10_6 and for different values of Nz,
10
-3
F(Xa)
10
-4
10
-5
Nz = 100, Nr = K =i 5 Mz. - 200, Nr ~ K ^: 10
O O O O O
10-6 L
ío-15
□ □ □
10"2 10"3 10-4 10~5 10-6
|QXCT + ¿bll
I ^ ¿ lloo
10-10 (a)
10"
10
-7
Nz = 100, ?Vr = = 5 Nz = 200, Nr = ft" = 10
10
-15
10-io (b)
10"
Fig. 15. (a) The value of the functional (3 . 2) and (b) the value of the uniform norm of the residual of the system (2 . 7), which have been calculated from the solution to the regularized system (3 . 14) for different values of a
Nr and K, are shown in Fig. 16. A comparison of Figs. 16b and 16e shows that the increase in the number of longitudinal and transversal harmonics Nz leads to a decrease in the level of divergence; the largest values of divergence arise at the entry to the computational region (see Figs. 16c and 16f). In this case, the vorticity distribution in the flow core changes significantly and the regions with vortices diminish. In particular, it can be seen from Fig. 16f that the largest value of vorticity arises on the wall of the channel. Figure 16a corresponds to the value of the functional (4 . 2) F & 4 . 0-10-5, and Figure 16d corresponds to the value of the functional (4 . 2) F & « 1. 8 • 10"5.
Nz = 100, Nr = K = 5
Nz = 200, Nr = K = 10
4 1
0 V 0.5
-4 0
0.02 1
D
-0.02 V 0.5
-0.04
-0.06 0
_ H 0.02 1 0
-0.02 y 0.5 -0.04
-0.06 o
Fig. 16. (a) and (d) vorticity field, (b) and (e) divergence distribution, (c) and (f) the initial area of divergence distribution. Regularization parameter a = 10_6
5. Conclusion
In this paper we present a method for constructing inhomogeneous velocity fields using expansions in terms of eigenfunctions of the Laplace operator. The components of the velocity field of the flow are approximated by the partial sums of the corresponding expansions, part of the expansion coefficients is fixed and the other part is determined from the problem of minimizing the integral of the squared divergence. It is shown that the minimization problem is incorrect and its solution requires regularization. Using the solution to the regularized problem, examples of plane-parallel and axisymmetric inhomogeneous flows are constructed. The following problems remain open:
1. As the number of harmonics in the partial sums approximating the velocity field components increases, high-frequency oscillations arise in the vorticity field. It is of interest to study the mechanism of these oscillations and the ways of suppressing them.
2. As the number of harmonics in the partial sums approximating the velocity field components increases, the condition number of the system that arises in solving the minimization problem increases fast to ~ 1016. It is of interest to determine the scaling of the coefficients of partial sums that ensures moderate values of the condition number.
3. Application of the proposed method to constructing two-dimensional and three-dimensional velocity fields in various curvilinear orthogonal coordinate systems.
The authors express their gratitude to E. M. Artemova and I. K. Marchevsky for the interest shown in this work and for discussing the results.
Conflict of interest
The authors declare that they have no conflicts of interest.
References
[1] Bogomolov, V.A., Dynamics of the Vorticity on a Sphere, Fluid Dynam., 1977, vol. 12, no. 6, pp. 863-870; see also: Izv. Akad. Nauk SSSR Mekh. Zhidk. Gaza, 1977, no. 6, pp. 57-65.
[2] Bogomolov, V.A., On Two-Dimensional Hydrodynamics of a Sphere, Izv. Atmos. Ocean. Phys, 1979, vol. 15, no. 1, pp. 18-22; see also: Izv. Akad. Nauk SSSR Fiz. Atmos. Okeana, 1979, vol. 15, no. 1, pp. 29-35.
[3] Borisov, A. V., Kilin, A. A., and Mamaev, I. S., A New Integrable Problem of Motion of Point Vortices on the Sphere, Nelin. Dinam., 2007, vol. 3, no. 2, pp. 211-223 (Russian).
[4] Borisov, A. V. and Mamaev, I.S., Mathematical Methods in the Dynamics of Vortex Structures, Izhevsk: R&C Dynamics, Institute of Computer Science, 2005 (Russian).
[5] Brady, J.F. and Acrivos, A., Steady Flow in a Channel or Tube with an Accelerating Surface Velocity. An Exact Solution to the Navier - Stokes Equations with Reverse Flow, J. Fluid Mech., 1981, vol. 112, pp. 127-150.
[6] Holodniok, M., Kubícek, M., and Hlavácek, V., Computation of the Flow between Two Rotating Coaxial Disks, J. Fluid Mech., 1977, vol. 81, no. 4, pp. 689-699.
[7] Holodniok, M., Kubícek, M., and Hlavácek, V., Computation of the Flow between Two Rotating Coaxial Disks: Multiplicity of Steady-State Solutions, J. Fluid Mech., 1981, vol. 108, pp. 227-240.
[8] Kamenetskiy, D. S., Bussoletti, J.E., Hilmes, C.L., Venkatakrishnan, V., Wigton, L.B., and Johnson, F.T., Numerical Evidence of Multiple Solutions for the Reynolds-Averaged Navier - Stokes Equations, AIAA J, 2014, vol. 52, no. 8, pp. 1686-1698.
[9] Kempka, S.N., Glass, M.W., Peery, J.S., Strickland, J.H., and Ingber, M.S., Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations, Sandia Report No. SAND-96-0583, Albuquerque, N.M.: Sandia National Lab.(SNL-NM), 1996.
[10] Kochin, N.E., Vector Calculus and the Principles of Tensor Calculus, Moscow: Nauka, 1965 (Russian).
[11] Miklavcic, M. and Wang, C., Viscous Flow due to a Shrinking Sheet, Quart. Appl. Math., 2006, vol. 64, no. 2, pp. 283-290.
[12] Phillips, D.L., A Technique for the Numerical Solution of Certain Integral Equations of the First Kind, J. Assoc. Comput. Mach, 1962, vol. 9, pp. 84-97.
[13] Tikhonov, A. N., On Incorrect Problems of Linear Algebra and a Stable Method for Their Solution, Dokl. Akad. Nauk SSSR, 1965, vol. 163, no. 3, pp. 591-594 (Russian).
[14] Tikhonov, A.N. and Arsenin, V. Ya., Methods for Solving Incorrectly Posed Problems, 2nd ed., Moscow: Nauka, 1979 (Russian).
[15] Watson, E. B.B., Banks, W.H. H., Zaturska, M.B., and Drazin, P. G., On Transition to Chaos in Two-Dimensional Channel Flow Symmetrically Driven by Accelerating Walls, J. Fluid Mech., 1990, vol. 212, pp. 451-485.
