Global Galerkin method for stability studies in incompressible CFD and other possible applications
SSchool of Mechanical Engineering, Faculty of Engineering, Tel-Aviv University, Tel-Aviv, Israel e-mail: [email protected]
To appear in “
Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics ”, ed. A. Gelfgat, Springer, 2018.
Global Galerkin method for stability studies in incompressible CFD and other possible applications Alexander Gelfgat
Abstract
In this paper the author reviews a version of the global Galerkin that was developed and applied in a series of earlier publications. The method is based on divergence-free basis functions satisfying all the linear and homogeneous boundary conditions. The functions are defined as linear superpositions of the Chebyshev polynomials of the first and second klinds that are combined into divergence free vectors. The description and explanations of treatment of boundary conditions inhomogeneities and singularities are given. Possible implementation for steady state solvers, path-continuation, stability solvers and straight-forward integration in time are discussed. The most important results obtained using the approach are briefly reviewed and possible future applications are discussed.
Keywords: weighted residual methods, spectral methods, Galerkin method, instability, path-continuation
1 Introduction
This paper revisits a version of the global Galerkin method whose development started in the author’s PhD thesis (Gelfgat, 1988) and resulted in an effective approach for stability analysis of model incompressible flows, so that well known and less known results had been published continuously between 1994 and 2005. Several years ago I was asked to present this approach again, and it triggered some interest among several young colleagues. It is mainly their interest that inspired me to review all that was done and to make a consistent description of all the technical details, which were distributed over several papers, and sometimes omitted. Naturally, I have added some comments on my personal current opinion of what was done and what possibly can be done in the future. The development of this approach started in the era when computer memory was measured in kilobytes, so that even storage of a large matrix was a problem. The only way to reduce the required memory was to reduce the number of degrees of freedom. Since discretization of the flow region (domain) requires fine grids, the decrease of number of degrees of freedom was sought in different versions of global weighted residual methods, called also spectral methods. These methods do not discretize the domain, and approximate solutions as truncated series of basis functions. It was shown by many authors that weighted residuals methods allow for a very significant reduction of degrees of freedom (see, e.g., Canuto et al., shown in Gelfgat & Tanasawa (1994). Later, as computational power grew, so did the possibilities of applications of the method. A small number of degrees of freedom resulted in Jacobian matrices of small size, which could be treated numerically for the computation of steady flows, as well as for the solution of the eigenvalue problems associated with linear stability (Gelfgat et al., 1997, 1999a,b). Divergence free functions were extended to cylindrical coordinates, which allowed us to consider flow in the rotating disk – cylinder system and to obtain the first stability results (Gelfgat et al., 1996) for this configuration. Since then a number of different problems addressing steady flows, multiplicity of states, and stability were solved for different flows in two-dimensional rectangular cavities and cylinders. The periodic circumferential coordinate in the cylindrical geometry allows one to study stability of an axisymmetric base flow with respect to three-dimensional perturbations. The linearized problem for the 3D disturbances separates for each circumferential Fourier mode, so that the final answer is obtained by consideration of several 2D-like problems. Subsequently, a considerable number of results were obtained for stability of flows in cylinders, which were driven by the rotation of boundaries, by buoyancy convection, and by magnetic field. The first attempt to study 3D instability in Cartesian coordinates was carried out in Gelfgat (1999) for the Rayleigh-Bénard problem in a rectangular box, i.e., the stability of a quiescent fluid heated from below. No attempts were made by the author to study stability of fully developed 3D flow. At the same time, the three-dimensional bases found an unexpected application for visualization of incompressible 3D flows (Gelfgat, 2014). In what follows we start from a brief description of the weighted residual and Galerkin method formulated for an incompressible fluid dynamics problem. Then the proposed way of constructing basis functions is explained, starting from a one-dimensional two-point problem. The treatment of inhomogeneities and boundary conditions singularities is explained. This follows by discussion of the resulting dynamical system and explanations of how it was treated in the cited works, as well as how it can be treated for fully 3D problems. After that, several illustrative examples are presented. Finally, a discussion of possible future studies is given.
2 The problem and the numerical method
We consider flow of an incompressible fluid in a two-dimensional rectangle ≤ 𝑥𝑥 ≤𝐴𝐴 𝑥𝑥 , 0 ≤ 𝑦𝑦 ≤ 𝐴𝐴 𝑦𝑦 , or in a three-dimensional box ≤ 𝑥𝑥 ≤ 𝐴𝐴 𝑥𝑥 , 0 ≤ 𝑦𝑦 ≤ 𝐴𝐴 𝑦𝑦 , 0 ≤ 𝑧𝑧 ≤ 𝐴𝐴 𝑧𝑧 . The rectangular shape of the domain is the main restriction for all the following. Below we discuss how this restriction can be relaxed to a canonical shape, i.e., the domain bounded by the coordinate surfaces belonging to a system of orthogonal curvilinear coordinates. The momentum and continuity equations for velocity 𝒗𝒗 and pressure 𝑝𝑝 read 𝜕𝜕𝒗𝒗𝜕𝜕𝜕𝜕 + ( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 = −∇𝑝𝑝 + ∆𝒗𝒗 + 𝒇𝒇 , (2.1) ∇ ∙ 𝒗𝒗 = 0 , (2.2) where 𝑅𝑅𝑅𝑅 is the Reynolds number. We assume also that boundary conditions for all three velocity components are linear and homogeneous, e.g., the no-slip conditions on all boundaries. Equations (2.1)-(2.2) can be considered together with other scalar transport equations for, e.g., temperature and/or concentration, an example of which will be given below. To formulate the Galerkin method we assume that the solution 𝒗𝒗 belongs to a space V of divergence-free vectors satisfying all the (linear and homogeneous) boundary conditions, LHBC. Assume that vectors { 𝝋𝝋 𝐾𝐾 } 𝐾𝐾=1∞ form a basis in this space. Then the solution 𝒗𝒗 can be represented as 𝒗𝒗 = ∑ 𝑐𝑐 𝐾𝐾 𝝋𝝋 𝐾𝐾∞𝐾𝐾=1 . (2.3) The coefficients 𝑐𝑐 𝐾𝐾 can be obtained by evaluation of inner products of Eq. (2.3) with a basis vector 𝝋𝝋 𝐿𝐿 : 〈𝒗𝒗 , 𝝋𝝋 𝐿𝐿 〉 = ∑ 𝑐𝑐 𝐾𝐾 〈𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝐿𝐿 〉 ∞𝐾𝐾=1 , 𝐿𝐿 = 1,2,3, … (2.4) Here we assumed that the space V is supplied with an inner product 〈∙ , ∙〉 , which is yet to be defined. If the functions { 𝝋𝝋 𝐾𝐾 } 𝑘𝑘=1∞ form an orthogonal basis, then the coefficients 𝑐𝑐 𝐾𝐾 are obtained as 𝑐𝑐 𝐾𝐾 = 〈𝒗𝒗 , 𝝋𝝋 𝐾𝐾 〉〈𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝐾𝐾 〉 , 𝐾𝐾 = 1,2,3, … (2.5) However, if the basis functions are not orthogonal, the expressions (2.4) form an infinite system of linear algebraic equations, which can be solved only up to a certain truncation. Keeping the first 𝑁𝑁 terms in the series (2.3) and defining a vector of coefficients 𝒄𝒄 = { 𝑐𝑐 , 𝑐𝑐 , 𝑐𝑐 , … , 𝑐𝑐 𝑁𝑁 } , the first 𝑁𝑁 coefficients can be obtained as 𝒄𝒄 = 𝐺𝐺 −1 𝒇𝒇 , 𝐺𝐺 𝐾𝐾𝐿𝐿 = 〈 𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝐿𝐿 〉 , 𝑓𝑓 𝐿𝐿 = 〈𝒗𝒗 , 𝝋𝝋 𝐿𝐿 〉 . (2.6) Here 𝐺𝐺 is the Gram matrix, and 𝑁𝑁 is the truncation number. In the relations (2.4)-(2.6) we assumed that the vector 𝒗𝒗 is known. If it is unknown, the coefficients 𝑐𝑐 𝐾𝐾 can be obtained only approximately by minimization of the residual of the momentum equation. To show how they can be calculated, we first assume that the coefficients 𝑐𝑐 𝐾𝐾 are time-dependent and the basis functions { 𝝋𝝋 𝐾𝐾 } 𝑘𝑘=1∞ depend only on coordinate values. Then the representation (2.3) defines a time- and space-dependent function 𝒗𝒗 . Note that the continuity equation (2.2), as well as all the boundary conditions are already satisfied because 𝒗𝒗 belongs to the space V . Clearly, the solution we are looking for also belongs to this space, so that the momentum equation is the only one to be solved. Now we rewrite the momentum equation (2.1) in two additional and equivalent forms 𝜕𝜕𝒗𝒗𝜕𝜕𝜕𝜕 + ( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 + ∇𝑝𝑝 − ∆𝒗𝒗 − 𝒇𝒇 = 𝑹𝑹 = (2.7) 𝜕𝜕𝒗𝒗𝜕𝜕𝜕𝜕 = −∇𝑝𝑝 + ∆𝒗𝒗 − ( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 + 𝒇𝒇 (2.8) Here 𝑹𝑹 is the residual of the momentum equation. If 𝒗𝒗 is the solution then 𝑹𝑹 ≡ . For a general case we assume that all possible residuals belong to a certain functional space W , which is also supplied with an inner product 〈∙ , ∙〉 𝑊𝑊 . Assume also that { 𝝓𝝓 𝐾𝐾 } 𝐾𝐾=1∞ is a basis in W . Then the requirement that the residual be zero is equivalent to the requirement that the residual 𝑹𝑹 is orthogonal to all the basis functions 𝝓𝝓 𝐾𝐾 , namely 〈𝑹𝑹 , 𝝓𝝓 𝐾𝐾 〉 𝑊𝑊 = 0, 𝐾𝐾 = 1, 2, 3, … (2.9) The equations (2.9) form an infinite set of non-linear time-dependent ODEs that must be solved to find the time-dependent coefficients 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) . However, this system of equations is not closed yet, since nothing is said about the pressure. To proceed, we assume that the pressure belongs to a space S of scalar time-dependent functions, differentiable at least twice in the domain, with the basis { 𝑠𝑠 𝐾𝐾 } 𝐾𝐾=1∞ . The pressure is represented as 𝑝𝑝 = ∑ 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) 𝑠𝑠 𝐾𝐾∞𝐾𝐾=1 . (2.10) The equation for pressure is formed in the standard way, by applying the divergence operator to both sides of the momentum equation (2.1), which yields ∆𝑝𝑝 = −𝑑𝑑𝑑𝑑𝑑𝑑 [( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 − 𝒇𝒇 ] . (2.11) The boundary conditions required this equation will be discussed later. Note, that the velocity representation (2.3), even in the truncated form used below, guarantees zero velocity divergence. Since the Laplacian and divergence operators are evaluated analytically, they commute in every approximate (truncated) formulation. We introduce the residual 𝐷𝐷 of the pressure equation 𝐷𝐷 = ∆𝑝𝑝 + 𝑑𝑑𝑑𝑑𝑑𝑑 [( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 − 𝒇𝒇 ] . (2.12) For the residual 𝐷𝐷 we demand only piecewise continuity in all spatial directions, so that it formally belongs to another space of scalar functions. This space is denoted as D , its basis as { 𝑞𝑞 𝐾𝐾 } 𝐾𝐾=1∞ , and the scalar product as 〈∙ , ∙〉 𝐷𝐷 . Now, for the solution of the problem represented by 𝒗𝒗 and 𝑝𝑝 , the scalar non-linear algebraic equations 〈𝐷𝐷 , 𝑞𝑞 𝐿𝐿 〉 𝐷𝐷 = 0, 𝐿𝐿 = 1, 2, 3, … (2.13) must be satisfied. Equations (2.13) define the coefficients 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) and must be satisfied together with the equations (2.9). Note that these equations do not contain the time derivative and, therefore, are algebraic constraints for the ODEs (2.9). Obviously, V ⊆ W and S ⊆ D . To build a numerical procedure for obtaining the first 𝑁𝑁 𝑣𝑣 and 𝑁𝑁 𝑝𝑝 coefficients of the velocity and pressure series (2.3) and (2.10), we truncate the series together with the residual projection equations (2.9) and (2.13). This results in 𝒗𝒗 ≈ ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 𝝋𝝋 𝐾𝐾𝑁𝑁 𝑣𝑣 𝐾𝐾=1 , 〈𝑹𝑹 , 𝝓𝝓 𝐿𝐿 〉 𝑊𝑊 = 0, 𝐿𝐿 = 1, 2, 3, … , 𝑁𝑁 𝑣𝑣 , (2.14) 𝑝𝑝 ≈ ∑ 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) 𝑠𝑠 𝐾𝐾𝑁𝑁 𝑝𝑝 𝐾𝐾=1 , 〈𝐷𝐷 , 𝑞𝑞 𝑀𝑀 〉 𝐷𝐷 = 0, 𝑀𝑀 = 1, 2, 3, … , 𝑁𝑁 𝑝𝑝 , (2.15) which defines 𝑁𝑁 𝑣𝑣 non-linear ODEs and 𝑁𝑁 𝑝𝑝 non-linear algebraic constraints for calculation of the time-dependent coefficients 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) and 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) . This is called method of weighted residuals (Fletcher (1984), Boyd (2000), Canuto et al. (2006)). The ODEs and the algebraic constraints resulting from projections of the momentum and pressure equation residuals onto the basis functions are ∑ 𝑐𝑐̇ 𝐾𝐾 ( 𝑡𝑡 ) 〈𝝋𝝋 𝐾𝐾 , 𝝓𝝓 𝑀𝑀 〉 𝑊𝑊𝑁𝑁 𝑣𝑣 𝐾𝐾=1 = − ∑ ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 𝑐𝑐 𝐿𝐿 ( 𝑡𝑡 ) 〈 ( 𝝋𝝋 𝐿𝐿 ∙ ∇ ) 𝝋𝝋 𝐾𝐾 , 𝝓𝝓 𝑀𝑀 〉 𝑊𝑊 − 𝑁𝑁 𝑣𝑣 𝐿𝐿=1𝑁𝑁 𝑣𝑣 𝐾𝐾=1 ∑ 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) 𝑁𝑁 𝑝𝑝 𝐾𝐾=1 〈∇𝑠𝑠 𝐾𝐾 , 𝝓𝝓 𝑀𝑀 〉 𝑊𝑊 + ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 〈∆𝝋𝝋 𝐾𝐾 , 𝝓𝝓 𝑀𝑀 〉 𝑊𝑊 + 𝑁𝑁 𝑣𝑣 𝐾𝐾=1 〈𝒇𝒇 , 𝝓𝝓 𝑀𝑀 〉 𝑊𝑊 , 𝑀𝑀 = 1, 2, 3, … , 𝑁𝑁 𝑣𝑣 (2.16) ∑ 𝑑𝑑 𝐾𝐾 ( 𝑡𝑡 ) 𝑁𝑁 𝑝𝑝 𝐾𝐾=1 〈∆𝑠𝑠 𝐾𝐾 , 𝑞𝑞 𝐽𝐽 〉 𝐷𝐷 = − ∑ ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 𝑐𝑐 𝐿𝐿 ( 𝑡𝑡 ) 〈∇ ∙ ( 𝝋𝝋 𝐿𝐿 ∙ ∇ ) 𝝋𝝋 𝐾𝐾 , 𝑞𝑞 𝐽𝐽 〉 𝐷𝐷 + 〈∇ ∙ 𝒇𝒇 , 𝑞𝑞 𝐽𝐽 〉 𝐷𝐷𝑁𝑁 𝑣𝑣 𝐿𝐿=1𝑁𝑁 𝑣𝑣 𝐾𝐾=1 (2.17) 𝐽𝐽 = 1, 2, 3, … , 𝑁𝑁 𝑝𝑝 Following common definitions, { 𝝋𝝋 𝐾𝐾 } 𝑘𝑘=1∞ and { 𝑠𝑠 𝐾𝐾 } 𝐾𝐾=1∞ are called coordinate basis systems, while { 𝝓𝝓 𝐾𝐾 } 𝑘𝑘=1∞ and { 𝑞𝑞 𝐾𝐾 } 𝐾𝐾=1∞ - are called projection basis systems.
Note that the weighted residuals method can be formulated also for coordinate basis functions that do not satisfy some or all boundary conditions. Fletcher (1984) distinguishes between coordinate functions that satisfy only a (linear) differential equation, satisfy only boundary conditions, and do not satisfy anything, calling these three cases boundary, interior, and mixed, respectively. Within this classification, and noting that the equations are non-linear, we discuss mainly interior methods. To arrive at the Galerkin formulation we assume that the boundary conditions do not explicitly involve time. Then the l.h.s. of Eq. (2.8) satisfies the LHBC of the velocity and is divergence-free. Then, the same can be said about the r.h.s. of Eq. (2.8), and finally about the residual of the momentum equation 𝑹𝑹 defined in Eq. (2.7). Thus, for time-independent boundary conditions, the residual 𝑹𝑹 belongs to the space V , so that we can choose 𝝓𝝓 𝐾𝐾 = 𝝋𝝋 𝐾𝐾 for all 𝐾𝐾 . Also assuming that the coordinate systems { 𝝋𝝋 𝐾𝐾 } 𝑘𝑘=1∞ and { 𝑠𝑠 𝐾𝐾 } 𝐾𝐾=1∞ are usually constructed from trigonometric functions or polynomials that are infinitely differentiable, the residual 𝐷𝐷 will also be differentiable infinite number of times. Since the physical problem does not specify any pressure boundary conditions, we can assume that both 𝑝𝑝 and 𝐷𝐷 belong to the same space S , and we can choose 𝑞𝑞 𝐾𝐾 = 𝑠𝑠 𝐾𝐾 for all 𝐾𝐾 . Thus, we arrive at a particular version of the weighted residuals method, in which the coordinate and projection basis systems coincide. This version is known as the Galerkin or Boubnov – Galerkin method. An obvious reason to choose the Galerkin method among all possible weighted residual formulations follows from the fact that V ⊆ W . By projecting onto a smaller space V , we expect that with increase of the truncation number convergence will be faster. Another, less obvious and more profound reason follows from the definition of the inner products via volume integrals. For scalar and vector functions defined in the domain 𝑉𝑉 and on its boundary we define 〈𝑓𝑓 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ), 𝑔𝑔 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 〉 𝜌𝜌 = ∫ 𝜌𝜌 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 𝑓𝑓𝑔𝑔𝑑𝑑𝑉𝑉 𝑉𝑉 , (2.18) 〈𝒖𝒖 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ), 𝒗𝒗 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 〉 𝜌𝜌 = ∫ 𝜌𝜌 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 𝒖𝒖 ∙ 𝒗𝒗 𝑑𝑑𝑉𝑉 𝑉𝑉 , (2.19) where 𝜌𝜌 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) > 0 is the weight function. The simplest and the most robust formulation is obtained with the unit weight 𝜌𝜌 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) = 1 , for which the inner products (2,18), (2.19) may also have some physical meaning. For example, the norm produced by (2.19) becomes twice the dimensionless kinetic energy. An important additional advantage of the unit weight follows from consideration of the inner product of the gradient of a scalar field 𝑓𝑓 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) with a divergence free vector field 𝒖𝒖 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) . Assume that 𝛤𝛤 is the domain boundary, 𝒏𝒏 its normal, and that the component of 𝒖𝒖 normal to the boundary vanishes. Keeping in mind that ∇ ∙ 𝒖𝒖 = 0 and 𝒖𝒖 ∙ 𝒏𝒏 = 0 , we have 〈∇𝑓𝑓 , 𝒖𝒖〉 = ∫ ∇𝑓𝑓 ∙ 𝒖𝒖𝑑𝑑𝑉𝑉 𝑉𝑉 = ∫ [ ∇ ∙ ( 𝑓𝑓𝒖𝒖 ) − 𝑓𝑓∇ ∙ 𝒖𝒖 ] 𝑉𝑉 𝑑𝑑𝑉𝑉 = ∫ ∇ ∙ ( 𝑓𝑓𝒖𝒖 ) 𝑉𝑉 𝑑𝑑𝑉𝑉 = ∫ 𝑓𝑓𝒖𝒖 ∙ 𝒏𝒏 𝛤𝛤 𝑑𝑑𝛤𝛤 = 0 (2.20) Thus, if the velocity does not penetrate the boundary, and the inner product is chosen as in (2.19), the projection of the pressure gradient on all the velocity basis functions vanishes. This means that equation systems (2.16) and (2.17) separate, so that velocity can be calculated from (2.14) without any knowledge about pressure. Note that this calculation involves both minimization of the residual and summation of the truncated velocity series. The pressure then can be computed from (2.15) using the previously found velocity field. In this case Eqs. (2.16), which minimizes the residual by computing the time-dependent coefficients 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) , become ∑ 𝑐𝑐̇ 𝐾𝐾 ( 𝑡𝑡 ) 〈𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝑀𝑀 〉 𝑣𝑣 𝐾𝐾=1 = − ∑ ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 𝑐𝑐 𝐿𝐿 ( 𝑡𝑡 ) 〈 ( 𝝋𝝋 𝐿𝐿 ∙ ∇ ) 𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝑀𝑀 〉 𝑣𝑣 𝐿𝐿=1𝑁𝑁 𝑣𝑣 𝐾𝐾=1 + ∑ 𝑐𝑐 𝐾𝐾 ( 𝑡𝑡 ) 〈∆𝝋𝝋 𝐾𝐾 , 𝝋𝝋 𝑀𝑀 〉 + 𝑁𝑁 𝑣𝑣 𝐾𝐾=1 〈𝒇𝒇 , 𝝋𝝋 𝑀𝑀 〉 (2.21) 𝑀𝑀 = 1, 2, 3, … , 𝑁𝑁 𝑣𝑣 This is the formal Galerkin formulation for calculating an approximate solution of a problem. Note that the exclusion of pressure by the Galerkin projection (2.20) is not restricted to closed regions with non-penetrative boundaries. The inhomogeneity in the velocity boundary conditions can be removed by change of variables, which is discussed below in more detail. To proceed we need to explain how to build basis functions, which are divergence-free and satisfy the whole set of LHBC. Then we will discuss how to handle inhomogeneous boundary conditions, curvilinear coordinates, and weight functions others than unity.
3 Basis functions
We start from the question of how to satisfy all the LHBC for a scalar unknown function. We assume that the basis functions for all three-dimensional time-dependent scalar variables, e.g., temperature and/or concentration, are represented as products of some one-dimensional bases. Thus for a function 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 , 𝑡𝑡 ) defined in a box ≤ 𝑥𝑥 ≤ 𝐴𝐴 𝑥𝑥 , 0 ≤ 𝑦𝑦 ≤ 𝐴𝐴 𝑦𝑦 , 0 ≤ 𝑧𝑧 ≤ 𝐴𝐴 𝑧𝑧 we seek a representation in the form of the tensor (Kronecker) product of one-dimensional bases 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 , 𝑡𝑡 ) = ∑ ∑ ∑ 𝑑𝑑 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) 𝑓𝑓 𝑖𝑖 ( 𝑥𝑥 ) 𝑔𝑔 𝑖𝑖 ( 𝑦𝑦 ) ℎ 𝑘𝑘 ( 𝑧𝑧 ) 𝑁𝑁 𝑥𝑥 𝑖𝑖=1𝑁𝑁 𝑦𝑦 𝑖𝑖=1𝑁𝑁 𝑧𝑧 𝑘𝑘=1 , (3.1) Here 𝑑𝑑 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) are unknown time-dependent coefficients, 𝑁𝑁 𝑥𝑥 , 𝑁𝑁 𝑦𝑦 , 𝑁𝑁 𝑧𝑧 are the truncation numbers specified in each spatial direction separately, and 𝑓𝑓 𝑖𝑖 ( 𝑥𝑥 ), 𝑔𝑔 𝑖𝑖 ( 𝑦𝑦 ), ℎ 𝑘𝑘 ( 𝑧𝑧 ) are one-dimensional bases that must be defined in each direction. The three-dimensional basis corresponding to those defined in the previous chapter is 𝐹𝐹 𝐽𝐽 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) = 𝑓𝑓 𝑖𝑖 ( 𝑥𝑥 ) 𝑔𝑔 𝑖𝑖 ( 𝑦𝑦 ) ℎ 𝑘𝑘 ( 𝑧𝑧 ), 𝐽𝐽 = 𝑁𝑁 𝑥𝑥 �𝑁𝑁 𝑦𝑦 ( 𝑘𝑘 −
1) + 𝑗𝑗 − � + 𝑑𝑑 (3.2) Starting from here we will use capital letters for the global indices, and small letters for the one-dimensional indices. If all the boundary conditions for 𝜃𝜃 are linear and homogeneous, the functions 𝑓𝑓 𝑖𝑖 ( 𝑥𝑥 ), 𝑔𝑔 𝑖𝑖 ( 𝑦𝑦 ) and ℎ 𝑘𝑘 ( 𝑧𝑧 ) must satisfy the boundary conditions posed in the x- , y- , and z- directions, respectively. Assume that there are M boundary conditions in the x- direction posed in the form ∑ 𝛼𝛼 𝑚𝑚𝑚𝑚 𝑓𝑓 ( 𝑚𝑚 ) ( 𝑥𝑥 𝑚𝑚 ) 𝐿𝐿𝑚𝑚=0 = 0, 𝑚𝑚 = 1, 2, … , 𝑀𝑀 . (3.3) Here 𝛼𝛼 𝑚𝑚𝑚𝑚 are known coefficients, 𝑙𝑙 is the derivative number, and 𝑥𝑥 𝑚𝑚 are the borders (that are parts of coordinate surfaces 𝑥𝑥 = 𝑥𝑥 𝑚𝑚 ) where the boundary conditions are posed. Note that m is the number of a boundary condition and not a number of the point, since several boundary conditions can be defined at the same point, so that sometimes 𝑥𝑥 𝑚𝑚 = 𝑥𝑥 𝑚𝑚 for 𝑚𝑚 ≠ 𝑚𝑚 . For the following we can also extend (3.3) by assuming that negative 𝑙𝑙 correspond to integrals, and that surfaces 𝑥𝑥 = 𝑥𝑥 𝑚𝑚 are not necessarily the boundary points, but also can lie inside the flow region, as it happens in a two-fluid example below. Now assume that a set of functions { 𝑠𝑠 𝑘𝑘 ( 𝑥𝑥 )} 𝑘𝑘=1∞ forms a basis in, say, 𝐶𝐶 ∞ ([0, 𝐴𝐴 𝑥𝑥 ]) . This can be, for example, a trigonometric Fourier basis, or a set of orthogonal polynomials defined on [0, 𝐴𝐴 𝑥𝑥 ] . For some reasons we choose this basis for representation of the solution, but the functions 𝑠𝑠 𝑘𝑘 ( 𝑥𝑥 ) do not satisfy the boundary conditions (3.3). We build an alternative basis by considering superpositions of 𝑀𝑀 + 1 consequent basis functions as 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 ) = ∑ 𝛽𝛽 𝑛𝑛𝑘𝑘 𝑠𝑠 𝑘𝑘+𝑛𝑛 ( 𝑥𝑥 ) 𝑀𝑀𝑛𝑛=0 (3.4) Substituting 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 ) into the boundary conditions (3.3) we obtain ∑ 𝛽𝛽 𝑛𝑛𝑘𝑘 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 𝑚𝑚 ) 𝑀𝑀𝑛𝑛=0 = ∑ 𝛼𝛼 𝑚𝑚𝑚𝑚 ∑ 𝛽𝛽 𝑛𝑛𝑘𝑘 𝑠𝑠 𝑘𝑘+𝑛𝑛 ( 𝑚𝑚 ) ( 𝑥𝑥 𝑚𝑚 ) 𝑀𝑀𝑛𝑛=0𝐿𝐿𝑚𝑚=0 = 0, 𝑚𝑚 = 1, 2, … , 𝑀𝑀 (3.5) For each 𝑘𝑘 the relations (3.5) form a system of M equations that can be used to find 𝑀𝑀 + 1 coefficients 𝛽𝛽 𝑛𝑛𝑘𝑘 . To make the equations solvable we choose 𝛽𝛽 , 𝑘𝑘 = 1 , that leaves us with M linear algebraic equations for M remaining coefficients 𝛽𝛽 𝑛𝑛𝑘𝑘 for every fixed 𝑘𝑘 . The matrix of this system will be regular if the boundary conditions (3.3) are independent. In this way, the functions 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 ) are fully defined and satisfy all the boundary conditions (3.3). Obviously they form a basis in the subspace 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠 { 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 )} ⊂ 𝐶𝐶 ∞ ([0, 𝐴𝐴 𝑥𝑥 ]) , which contains functions that satisfy conditions (3.3) and are differentiable infinite number of times. In this way we form bases 𝑓𝑓 𝑖𝑖 ( 𝑥𝑥 ), 𝑔𝑔 𝑖𝑖 ( 𝑦𝑦 ) and ℎ 𝑘𝑘 ( 𝑧𝑧 ) for the representation (3.1). In the following all the examples will be based on the Chebyshev polynomials of the first and the second type, 𝑇𝑇 𝑘𝑘 ( 𝑥𝑥 ) and 𝑈𝑈 𝑘𝑘 ( 𝑥𝑥 ) , briefly described in the Appendix A. In the further text these polynomials will always be chosen as the basis { 𝑠𝑠 𝑘𝑘 ( 𝑥𝑥 )} 𝑘𝑘=1∞ . Example: Two-point boundary value problem
Consider a two-point boundary value problem for 𝜃𝜃 ( 𝑥𝑥 ) posed on the interval ≤ 𝑥𝑥 ≤ with the two ( 𝑀𝑀 = 2 ) boundary conditions 𝑠𝑠 𝜃𝜃 ′ (0) + 𝑏𝑏 𝜃𝜃 (0) = 0, 𝑠𝑠 𝜃𝜃 ′ (1) + 𝑏𝑏 𝜃𝜃 (1) = 0, (3.1.1) where 𝑠𝑠 , 𝑏𝑏 , 𝑠𝑠 , 𝑏𝑏 are known coefficients. Recalling the Chebyshev polynomials and taking into account 𝛽𝛽 , 𝑘𝑘 = 1 , our new basis functions are defined as 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 ) = ∑ 𝛽𝛽 𝑛𝑛𝑘𝑘 𝑇𝑇 𝑘𝑘+𝑛𝑛 ( 𝑥𝑥 ) = 𝑇𝑇 𝑘𝑘 ( 𝑥𝑥 ) + 𝛽𝛽 𝑇𝑇 𝑘𝑘+1 ( 𝑥𝑥 ) + 𝛽𝛽 𝑇𝑇 𝑘𝑘+2 ( 𝑥𝑥 ) . (3.1.2) The boundary values of the Chebyshev polynomials and their derivatives are given in Appendix A. Substituting 𝑟𝑟 𝑘𝑘 ( 𝑥𝑥 ) into the boundary conditions (3.1.1) and using Eqs. (A3) and (A5), we obtain two equations for the coefficients 𝛽𝛽 and 𝛽𝛽 : −𝛽𝛽 [2 𝑠𝑠 ( 𝑘𝑘 + 1) + 𝑏𝑏 ] + 𝛽𝛽 [2 𝑠𝑠 ( 𝑘𝑘 + 2) + 𝑏𝑏 ] = − 𝑠𝑠 𝑘𝑘 − 𝑏𝑏 (3.1.3) 𝛽𝛽 [2 𝑠𝑠 ( 𝑘𝑘 + 1) + 𝑏𝑏 ] + 𝛽𝛽 [2 𝑠𝑠 ( 𝑘𝑘 + 2) + 𝑏𝑏 ] = − 𝑠𝑠 𝑘𝑘 − 𝑏𝑏 (3.1.4) These equation are easily solved analytically, which yields 𝛽𝛽 = �2𝑎𝑎 𝑘𝑘 +𝑏𝑏 ��2𝑎𝑎 ( 𝑘𝑘+2 ) +𝑏𝑏 �−�2𝑎𝑎 𝑘𝑘 +𝑏𝑏 ��2𝑎𝑎 ( 𝑘𝑘+2 ) +𝑏𝑏 � [ ( 𝑘𝑘+1 ) +𝑏𝑏 ][ ( 𝑘𝑘+2 ) +𝑏𝑏 ] + [ ( 𝑘𝑘+1 ) +𝑏𝑏 ][ ( 𝑘𝑘+2 ) +𝑏𝑏 ] , (3.1.5) 𝛽𝛽 = �2𝑎𝑎 𝑘𝑘 +𝑏𝑏 ��2𝑎𝑎 ( 𝑘𝑘+1 ) +𝑏𝑏 �+�2𝑎𝑎 𝑘𝑘 +𝑏𝑏 ��2𝑎𝑎 ( 𝑘𝑘+1 ) +𝑏𝑏 � [ ( 𝑘𝑘+1 ) +𝑏𝑏 ][ ( 𝑘𝑘+2 ) +𝑏𝑏 ] + [ ( 𝑘𝑘+1 ) +𝑏𝑏 ][ ( 𝑘𝑘+2 ) +𝑏𝑏 ] , (3.1.6) and defines a new basis, whose satisfy both boundary conditions. This idea was introduced in Orszag (1971a) for two-point homogeneous Dirichlet boundary conditions, and in Orszag (1971b) for boundary conditions of the Orr-Sommerfeld equation. It was formalized for an arbitrary set of boundary conditions (3.1.1) in Gelfgat (1988). Note that for a similar problem defined on the interval ≤ 𝑥𝑥 ≤ 𝐴𝐴 𝑥𝑥 we need only to replace 𝑥𝑥 by 𝑥𝑥 𝐴𝐴 𝑥𝑥 ⁄ in Eq. (3.1.2), which will not affect the expressions (3.1.5) and (3.1.6). To construct basis functions, which are two-dimensional divergence-free vectors satisfying all the boundary conditions, we start by constructing a divergence-free basis that does not yet involve any boundary conditions. Using the Chebyshev polynomials, and assuming the flow region is a rectangle ≤ 𝑥𝑥 ≤ 𝐴𝐴 𝑥𝑥 , 0 ≤ 𝑦𝑦 ≤ 𝐴𝐴 𝑦𝑦 , we define 𝒘𝒘 𝑖𝑖𝑖𝑖 = �𝑤𝑤 𝑖𝑖𝑖𝑖 ( 𝑥𝑥 ) 𝑤𝑤 𝑖𝑖𝑖𝑖 ( 𝑦𝑦 ) � = � 𝐴𝐴 𝑥𝑥 𝑇𝑇 𝑖𝑖 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � 𝑈𝑈 𝑖𝑖−1 � 𝑦𝑦𝐴𝐴 𝑦𝑦 �− 𝐴𝐴 𝑦𝑦 𝑈𝑈 𝑖𝑖−1 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � 𝑇𝑇 𝑖𝑖 � 𝑦𝑦𝐴𝐴 𝑦𝑦 �� , 𝑑𝑑 , 𝑗𝑗 = 1,2,3, … (3.2.1) 𝒘𝒘 = � 𝐴𝐴 𝑥𝑥 𝑈𝑈 𝑖𝑖−1 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � � , 𝒘𝒘 𝑖𝑖0 = � 𝐴𝐴 𝑦𝑦 𝑈𝑈 𝑖𝑖−1 � 𝑥𝑥𝐴𝐴 𝑥𝑥 �� (3.2.2) Applying the relation (A2), it is easily seen that ∇ ∙ 𝒘𝒘 𝑖𝑖𝑖𝑖 = 0 . Now, to implement the boundary conditions and to keep the divergence zero, we keep the x- and y- components of the vector dependent on each other, as in (3.2.1). To do so we implement all the velocity boundary conditions in the x- direction in the x- dependent part of the vector, and do the same in the y- direction. In the x- direction the boundary conditions are posed at 𝑥𝑥 = 0 and 𝑥𝑥 = 𝐴𝐴 𝑥𝑥 for the two velocity components, so that we have four boundary conditions in total. The same is done in the y- direction. Thus, we construct the basis using linear superpositions of five (4+1) consecutive Chebyshev polynomials. This results in 𝝋𝝋 𝑖𝑖𝑖𝑖 ( 𝑥𝑥 , 𝑦𝑦 ) = � 𝐴𝐴 𝑥𝑥 ∑ 𝜎𝜎 𝑖𝑖𝑖𝑖 𝑖𝑖+𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝜏𝜏 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � − 𝐴𝐴 𝑦𝑦 ∑ 𝜎𝜎 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝜏𝜏 𝑗𝑗𝑖𝑖 𝑖𝑖+𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � � (3.2.3) It is easy to check that ∇ ∙ 𝝋𝝋 𝑖𝑖𝑖𝑖 = 0 . The coefficients 𝜎𝜎 𝑖𝑖𝑚𝑚 and 𝜏𝜏 𝑖𝑖𝑚𝑚 must be found by substituting 𝝋𝝋 𝑖𝑖𝑖𝑖 into the boundary conditions. For example, assume that the rectangle has a stress free boundary at 𝑦𝑦 = 𝐴𝐴 𝑦𝑦 , while all the other boundaries are no-slip. This leads to the following boundary conditions for the velocity 𝒗𝒗 : 𝒗𝒗 ( 𝑥𝑥 = 0, 𝑦𝑦 ) = 𝒗𝒗 ( 𝑥𝑥 = 𝐴𝐴 𝑥𝑥 , 𝑦𝑦 ) = 𝒗𝒗 ( 𝑥𝑥 , 𝑦𝑦 = 0) = 0, 𝑑𝑑 𝑦𝑦 �𝑥𝑥 , 𝑦𝑦 = 𝐴𝐴 𝑦𝑦 � = 𝜕𝜕𝑣𝑣 𝑥𝑥 𝜕𝜕𝑦𝑦 �𝑥𝑥 , 𝑦𝑦 = 𝐴𝐴 𝑦𝑦 � = 0 (3.2.4) The assignment 𝜎𝜎 𝑖𝑖0 = 𝜏𝜏 𝑑𝑑 = 1 and substitution of (3.2.3) into (3.2.4) yields 𝜎𝜎 𝑖𝑖1 = 𝜎𝜎 𝑖𝑖3 = 0, 𝜎𝜎 = − , 𝜎𝜎 = (3.2.5) 𝜎𝜎 𝑖𝑖2 = − 𝑖𝑖𝑖𝑖+2 − ( 𝑖𝑖+1 )( 𝑖𝑖+4 ) 𝑖𝑖 ( 𝑖𝑖+2 )( 𝑖𝑖+3 ) , 𝜎𝜎 𝑖𝑖4 = ( 𝑖𝑖+1 )( 𝑖𝑖+4 ) 𝑖𝑖 ( 𝑖𝑖+3 ) , 𝑑𝑑 > 0 (3.2.6) 𝜏𝜏 = , 𝜏𝜏 𝑖𝑖1 = 2 𝑖𝑖 +2𝑖𝑖+1𝑖𝑖 +5𝑖𝑖 +7𝑖𝑖 , 𝑑𝑑 > 0 (3.2.7) 𝜏𝜏 = − , 𝜏𝜏 𝑖𝑖2 = − 𝑖𝑖 +8𝑖𝑖 +26𝑖𝑖 +40𝑖𝑖+24𝑖𝑖 +8𝑖𝑖 +22𝑖𝑖 +21𝑖𝑖 , 𝑑𝑑 > 0 (3.2.8) 𝜏𝜏 = − , 𝜏𝜏 𝑖𝑖3 = − 𝑖𝑖 +4𝑖𝑖+3𝑖𝑖 +5𝑖𝑖 +7𝑖𝑖 , 𝑑𝑑 > 0 (3.2.9) 𝜏𝜏 = , 𝜏𝜏 𝑖𝑖4 = 𝑖𝑖 +8𝑖𝑖 +22𝑖𝑖 +27+12𝑖𝑖 +8𝑖𝑖 +22𝑖𝑖 +21𝑖𝑖 , 𝑑𝑑 > 0 (3.2.10) Generalization of the 2D basis functions for a three-dimensional case is not straight-forward, since it is unclear how to produce divergence-free vectors, similar to those defined in (3.2.3), that will form a complete set of basis functions. To use a similar approach, we need to carry out several additional calculations. Assume that 𝒗𝒗 = ( 𝑢𝑢 , 𝑑𝑑 , 𝑤𝑤 ) is a divergence-free three-dimensional vector that satisfies the no-slip conditions on all the boundaries of the rectangular box ≤ 𝑥𝑥 ≤ 𝐴𝐴 𝑥𝑥 , 0 ≤ 𝑦𝑦 ≤ 𝐴𝐴 𝑦𝑦 , 0 ≤ 𝑧𝑧 ≤𝐴𝐴 𝑧𝑧 . Since ∇ ∙ 𝒗𝒗 = 𝜕𝜕𝑢𝑢 𝜕𝜕𝑥𝑥⁄ + 𝜕𝜕𝑑𝑑 𝜕𝜕𝑦𝑦⁄ + 𝜕𝜕𝑤𝑤 𝜕𝜕𝑧𝑧⁄ = 0 , so that 𝑤𝑤 = − ∫ ( 𝜕𝜕𝑢𝑢 𝜕𝜕𝑥𝑥⁄ + 𝜕𝜕𝑑𝑑 𝜕𝜕𝑦𝑦⁄ ) 𝑑𝑑𝑧𝑧 , a 3D incompressible vector field can be decomposed as 𝒗𝒗 = � 𝑢𝑢𝑑𝑑𝑤𝑤� = � 𝑢𝑢 𝑤𝑤 � + � 𝑑𝑑𝑤𝑤 � , 𝑤𝑤 = − ∫ 𝜕𝜕𝜕𝜕𝜕𝜕𝑥𝑥 𝑑𝑑𝑧𝑧 , 𝑧𝑧0 𝑤𝑤 = − ∫ 𝜕𝜕𝑣𝑣𝜕𝜕𝑦𝑦 𝑑𝑑𝑧𝑧 𝑧𝑧0 (3.3.1) Since both vectors 𝒗𝒗 ( 𝑥𝑥 , 𝑧𝑧 ) = ( 𝑢𝑢 , 0, 𝑤𝑤 ) and 𝒗𝒗 ( 𝑦𝑦 , 𝑧𝑧 ) = (0, 𝑑𝑑 , 𝑤𝑤 ) are divergence-free, this decomposition shows that a divergence-free velocity field can be represented as superposition of two fields which have components only in the ( x,z ) or ( y,z ) planes. For each of the two vectors we can construct basis functions similar to (3.2.3) 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) = ⎣⎢⎢⎢⎡ 𝐴𝐴 𝑥𝑥 ∑ 𝑎𝑎� 𝑖𝑖𝑖𝑖 ( 𝑖𝑖+𝑚𝑚 ) 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̂ 𝑘𝑘𝑛𝑛 𝑈𝑈 𝑘𝑘+𝑛𝑛−1 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � 𝑛𝑛�𝑛𝑛=0 − 𝐴𝐴 𝑧𝑧 ∑ 𝑠𝑠� 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−14𝑚𝑚=0 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̂ 𝑘𝑘𝑘𝑘 ( 𝑘𝑘+𝑛𝑛 ) 𝑇𝑇 𝑘𝑘+𝑛𝑛 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � 𝑛𝑛�𝑛𝑛=0 ⎦⎥⎥⎥⎤ (3.3.2) 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) = ⎣⎢⎢⎢⎢⎡ 𝐴𝐴 𝑦𝑦 ∑ 𝑠𝑠� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑗𝑗𝑖𝑖 ( 𝑖𝑖+𝑚𝑚 ) 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̃ 𝑘𝑘𝑛𝑛 𝑈𝑈 𝑘𝑘+𝑛𝑛−1 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � 𝑛𝑛�𝑛𝑛=0 − 𝐴𝐴 𝑧𝑧 ∑ 𝑠𝑠� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̃ 𝑘𝑘𝑘𝑘 ( 𝑘𝑘+𝑛𝑛 ) 𝑇𝑇 𝑘𝑘+𝑛𝑛 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � 𝑛𝑛�𝑛𝑛=0 ⎦⎥⎥⎥⎥⎤ (3.3.3) As above, the coefficients 𝑠𝑠� 𝑖𝑖𝑚𝑚 , 𝑏𝑏� 𝑖𝑖𝑚𝑚 , 𝑐𝑐̂ 𝑘𝑘𝑛𝑛 , 𝑠𝑠� 𝑖𝑖𝑚𝑚 , 𝑏𝑏� 𝑖𝑖𝑚𝑚 , and 𝑐𝑐̃ 𝑘𝑘𝑛𝑛 are defined after substitution of the functions in the boundary conditions. Note that the number of polynomials included in the linear superpositions in z-direction, 𝑠𝑠� + 1 , is not yet defined. This is because at 𝑧𝑧 = 𝐴𝐴 𝑧𝑧 only the sum 𝑤𝑤 ( 𝑥𝑥 , 𝑦𝑦 , 𝐴𝐴 𝑧𝑧 ) + 𝑤𝑤 ( 𝑥𝑥 , 𝑦𝑦 , 𝐴𝐴 𝑧𝑧 ) is zero, while the boundary values of 𝑤𝑤 ( 𝑥𝑥 , 𝑦𝑦 , 𝐴𝐴 𝑧𝑧 ) and 𝑤𝑤 ( 𝑥𝑥 , 𝑦𝑦 , 𝐴𝐴 𝑧𝑧 ) are not known. An example of such a decomposition can be found in Gelfgat (2014). Therefore there are only three boundary conditions in the z- direction to be satisfied by the basis functions 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) and 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) , so that 𝑠𝑠� = 3 . It is still possible to use these bases, but to satisfy the boundary conditions for 𝑤𝑤 at 𝑧𝑧 = 𝐴𝐴 𝑧𝑧 , additional algebraic constraints will be needed. Note that there is no such a problem if boundary conditions in the z- direction are periodic. To remove the above algebraic constraint at 𝑧𝑧 = 𝐴𝐴 𝑧𝑧 , we assume that 𝑠𝑠� = 4 in (3.3.2) and (3.3.3), so that the functions 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) and 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) are divergence-free and satisfy all the boundary conditions. In this case, using (3.2.5) and (3.2.6), 𝑠𝑠� 𝑖𝑖𝑚𝑚 = 𝑠𝑠� 𝑖𝑖𝑚𝑚 = 𝑏𝑏� 𝑖𝑖𝑚𝑚 = 𝑏𝑏� 𝑖𝑖𝑚𝑚 = 𝑐𝑐̂ 𝑖𝑖𝑚𝑚 = 𝑐𝑐̃ 𝑖𝑖𝑚𝑚 = 𝜎𝜎 𝑖𝑖𝑚𝑚 , 𝑏𝑏� 𝑖𝑖1 = 𝑏𝑏� 𝑖𝑖1 = 𝑏𝑏� 𝑖𝑖3 = 𝑏𝑏� 𝑖𝑖3 = 𝑏𝑏� 𝑖𝑖4 = 𝑏𝑏� 𝑖𝑖4 = 0 , 𝑏𝑏� 𝑖𝑖2 = 𝑏𝑏� 𝑖𝑖2 = − ( 𝑖𝑖+2 ) 𝑖𝑖 (3.3.4) Projection of the solution 𝒗𝒗 onto 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠�𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) � and 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠 �𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) � results in a vector 𝒗𝒗� 𝒗𝒗� = � 𝑢𝑢� 𝑤𝑤� � + � 𝑑𝑑�𝑤𝑤� � , (3.3.5) which is divergence-free, satisfies all the boundary conditions, but may not be a good approximation of 𝒗𝒗 because the set of basis functions still is not complete. To complete the basis we notice that 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠�𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) � and 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠 �𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) � project the velocity on the ( 𝑥𝑥 , 𝑧𝑧 ) and ( 𝑦𝑦 , 𝑧𝑧 ) planes, so that it is straight-forward to add projections on the ( 𝑥𝑥 , 𝑦𝑦 ) planes as well. Thus, similarly to (3.3.2) and (3.3.3) we add another set of divergence-free basis functions satisfying all the boundary conditions 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) = ⎣⎢⎢⎢⎢⎡ 𝐴𝐴 𝑥𝑥 ∑ 𝑎𝑎� 𝑖𝑖𝑖𝑖 ( 𝑖𝑖+𝑚𝑚 ) 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̅ 𝑘𝑘𝑛𝑛 𝑇𝑇 𝑘𝑘+𝑛𝑛 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � 𝑦𝑦 ∑ 𝑠𝑠� 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑥𝑥𝐴𝐴 𝑥𝑥 � ∑ 𝑏𝑏� 𝑗𝑗𝑖𝑖 ( 𝑖𝑖+𝑚𝑚 ) 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑦𝑦𝐴𝐴 𝑦𝑦 � ∑ 𝑐𝑐̅ 𝑘𝑘𝑛𝑛 𝑇𝑇 𝑘𝑘+𝑛𝑛 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � ⎦⎥⎥⎥⎥⎤ (3.3.6) Similarly to previous functions, for the no-slip boundary conditions posed on all boundaries, 𝑠𝑠� 𝑖𝑖𝑚𝑚 = 𝑏𝑏� 𝑖𝑖𝑚𝑚 = 𝜎𝜎 𝑖𝑖𝑚𝑚 , 𝑐𝑐̅ 𝑘𝑘1 = 𝑐𝑐̅ 𝑘𝑘3 = 𝑐𝑐̅ 𝑘𝑘4 = 0 , 𝑐𝑐̅ 𝑘𝑘2 = − ( 𝑘𝑘+2 ) 𝑘𝑘 (3.3.7) Finally, we seeek a three-dimensional solution of the form 𝒗𝒗 ≈ � � � 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) 𝑁𝑁 ( 𝑥𝑥 , 𝑦𝑦 ) 𝑘𝑘=0𝑀𝑀 ( 𝑥𝑥 , 𝑦𝑦 ) 𝑖𝑖=0 ( 𝑡𝑡 ) 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 𝐿𝐿 ( 𝑥𝑥 , 𝑦𝑦 ) 𝑖𝑖=0 + + � � � 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) 𝑁𝑁 ( 𝑥𝑥 , 𝑧𝑧 ) 𝑘𝑘=0𝑀𝑀 ( 𝑥𝑥 , 𝑧𝑧 ) 𝑖𝑖=0 ( 𝑡𝑡 ) 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 𝐿𝐿 ( 𝑥𝑥 , 𝑧𝑧 ) 𝑖𝑖=0 + (3.3.8) + � � � 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) 𝑘𝑘=0𝑀𝑀 ( 𝑦𝑦 , 𝑧𝑧 ) 𝑖𝑖=0 ( 𝑡𝑡 ) 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) 𝐿𝐿 ( 𝑦𝑦 , 𝑧𝑧 ) 𝑖𝑖=0 Projection of the residual of the momentum equation on all three sets of basis vectors yields three sets of ODEs for calculation of the three sets of time-dependent coefficients 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) , and 𝐶𝐶 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) . Since the basis functions are not orthogonal it will be necessary to invert the Gram matrix that is formed as 𝐺𝐺 = ⎣⎢⎢⎢⎢⎡〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑦𝑦 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑧𝑧 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑦𝑦 , 𝑧𝑧 ) 〉〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑦𝑦 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑧𝑧 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑦𝑦 , 𝑧𝑧 ) 〉〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑦𝑦 , 𝑥𝑥 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑦𝑦 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑦𝑦 , 𝑥𝑥 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑥𝑥 , 𝑧𝑧 ) 〉 〈 𝝋𝝋 𝑑𝑑𝑗𝑗𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) , 𝝋𝝋 𝑝𝑝𝑞𝑞𝑟𝑟 ( 𝑦𝑦 , 𝑧𝑧 ) 〉⎦⎥⎥⎥⎥⎤ (3.3.9) A simple numerical test for the no-slip boundary conditions and equal truncation numbers (starting from 4 and larger) in each direction and for each set of the functions shows that the Gram matrix is singular. This means that some of the functions are linearly dependent and must be excluded. Based on the above discussion, we see that in the case of periodic boundary conditions in the z- direction, all the set (3.3.6) can be omitted. However, some functions of this set must be added in the case of no-slip boundaries. This shows that the complete set of linearly independent basis functions differs for different boundary conditions. Unfortunately, the author could not arrive at a rigorous mathematical procedure that would determine which functions must be excluded at certain boundary conditions. At the same time, a simple numerical experiment can be helpful. Considering no-slip conditions at all the boundaries, we varied 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) in the last sum of Eq. (3.3.7). In other words, we varied the truncation number in the z- direction for the functions defined in (3.3.6) only. We found that by taking 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) = 0 , i.e., one basis function in the z- direction, we obtain a regular Gram matrix. Increase of 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) to 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) ≥ , makes the Gram function singular. Furthermore, taking a single basis function in the z- direction, with the third index 𝑘𝑘 ≥ , which means polynomials of larger degrees in (3.3.6), also results in a singular Gram matrix. This shows that the addition of the first polynomial in the z- direction, corresponding to 𝑁𝑁 ( 𝑦𝑦 , 𝑧𝑧 ) = 0 , is essential, while the others can be omitted. Using (3.3.7), this first polynomial is 𝑧𝑧 − 𝑧𝑧 ) . Returning to the sets (3.3.2) and (3.3.3) with the coefficients defined in (3.3.4), we observe that the polynomials corresponding to the x- and y- vector components start form the second degree, while those corresponding to the z-component start from the third one. Thus, the missing second-order polynomial must be added with the help of another set (3.3.6). Consider flow in a cylinder with radius R and height H . The whole problem is defined now in the cylindrical coordinates ( 𝑟𝑟 , 𝜃𝜃 , 𝑧𝑧 ), 0 ≤ 𝑟𝑟 ≤ 𝑅𝑅 , 0 ≤ 𝜃𝜃 ≤ 𝜋𝜋 , 0 ≤ 𝑧𝑧 ≤ 𝐻𝐻 . Using 𝜋𝜋 -periodicity in the azimuthal direction we represent the flow as a Fourier series 𝒗𝒗 = ∑ 𝒗𝒗 𝑘𝑘 ( 𝑟𝑟 , 𝑧𝑧 ) 𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝑑𝑑𝑘𝑘𝜃𝜃 ) ∞𝑘𝑘=−∞ , (3.4.1) so that the basis functions in the 𝜃𝜃 -direction are complex exponents 𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝑑𝑑𝑘𝑘𝜃𝜃 ) . The continuity equation for 𝒗𝒗 𝑘𝑘 ( 𝑟𝑟 , 𝑧𝑧 ) = �𝑢𝑢 𝑘𝑘 ( 𝑟𝑟 , 𝑧𝑧 ), 𝑑𝑑 𝑘𝑘 ( 𝑟𝑟 , 𝑧𝑧 ), 𝑤𝑤 𝑘𝑘 ( 𝑟𝑟 , 𝑧𝑧 ) � is ( 𝑟𝑟𝜕𝜕 𝑘𝑘 ) 𝜕𝜕𝑟𝑟 + 𝑖𝑖𝑘𝑘𝑟𝑟 𝑑𝑑 𝑘𝑘 + 𝜕𝜕𝑤𝑤 𝑘𝑘 𝜕𝜕𝑧𝑧 = 0 (3.4.2) Here we must distinguish between the axisymmetric case and axisymmetric Fourier mode, for which 𝑘𝑘 = 0 , and all the other 𝑘𝑘 ≠ modes. Axisymmetric flow (or the axisymmetric mode) is represented by a single set of the basis functions, which is built similarly to the above 2D Cartesian case, but taking into account the continuity equation (3.4.2). Note that if an axisymmetric flow has also a non-zero azimuthal component, the latter can be treated as a scalar function. Then the axisymmetric vector basis should include only the radial and axial components. An example of such basis, successfully used in several studies (see below) is 𝑼𝑼 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) = �𝑢𝑢 𝑤𝑤 � = ⎩⎨⎧
12 𝑟𝑟𝑅𝑅 ∑ 𝜎𝜎 𝑖𝑖𝑖𝑖 𝑖𝑖+𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � − 𝐴𝐴 𝑧𝑧 ∑ 𝜎𝜎 𝑖𝑖𝑚𝑚 𝑈𝑈� 𝑖𝑖+𝑚𝑚−1 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇 𝑗𝑗𝑖𝑖 𝑖𝑖+𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � ⎭⎬⎫ , (3.4.3) where 𝑈𝑈� 𝑛𝑛 � 𝑟𝑟𝑅𝑅 � = 𝑇𝑇 𝑛𝑛+1 � 𝑟𝑟𝑅𝑅 � + 𝑟𝑟 ( 𝑠𝑠 + 1) 𝑈𝑈 𝑛𝑛 � 𝑟𝑟𝑅𝑅 � . As above, the zero divergence of 𝑼𝑼 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) follows from Eq. (A2), and the coefficients 𝜎𝜎 𝑖𝑖𝑚𝑚 and 𝜇𝜇 𝑖𝑖𝑚𝑚 are obtained by substitution of 𝑼𝑼 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) in the boundary conditions. Note also that the r -component of velocity vanishes at the polar axis 𝑟𝑟 = 0 , so that for flow in a full cylinder only three conditions in the radial direction must be additionally satisfied. If the domain is a cylindrical layer (e.g., Taylor-Couette flow) then the polar axis is cut out and one remains with the four boundary conditions, as in the Cartesian case. Now consider Fourier modes of (3.4.1) corresponding to 𝑘𝑘 ≠ . Using the same idea as in Eq. (3.1.1) we observe that 𝒗𝒗 𝑘𝑘 = �𝑢𝑢 𝑘𝑘 𝑑𝑑 𝑘𝑘 𝑤𝑤 𝑘𝑘 � = �− 𝑢𝑢 𝑘𝑘1𝑖𝑖𝑘𝑘 𝜕𝜕 ( 𝑟𝑟𝜕𝜕 𝑘𝑘 ) 𝜕𝜕𝑟𝑟 � + �− 𝑟𝑟𝑖𝑖𝑘𝑘 𝜕𝜕𝑤𝑤 𝑘𝑘 𝜕𝜕𝑧𝑧 𝑤𝑤 𝑘𝑘 � . (3.4.4) Obviously, the sum of two azimuthal components of this relation satisfies the boundary conditions for 𝒗𝒗 𝑘𝑘 . It is not clear, however, whether each of them satisfies the boundary conditions separately. The author is not sure that this can be proved for a general case, but it can be easily examined for no-slip conditions at 𝑟𝑟 = 𝑅𝑅 and 𝑧𝑧 = 0, 𝐻𝐻 . Since 𝑤𝑤 𝑘𝑘 ( 𝑟𝑟 = 𝑅𝑅 ) = 0 the derivative 𝜕𝜕𝑤𝑤 𝑘𝑘 𝜕𝜕𝑧𝑧⁄ also vanishes at 𝑟𝑟 = 𝑅𝑅 . Since the sum of two azimuthal components vanishes at 𝑟𝑟 = 𝑅𝑅 , the azimuthal component of the first vector of the r.h.s also vanishes there, so that each component satisfies boundary conditions in the radial direction. Similarly, we consider 𝜕𝜕 ( 𝑟𝑟𝑢𝑢 𝑘𝑘 ) 𝜕𝜕𝑟𝑟⁄ at 𝑧𝑧 = 0, 𝐻𝐻 and conclude that it vanishes there because 𝑢𝑢 𝑘𝑘 ( 𝑟𝑟 , 0) = 𝑢𝑢 𝑘𝑘 ( 𝑟𝑟 , 𝐻𝐻 ) = 0 . Then also the second azimuthal component vanishes at 𝑧𝑧 = 0, 𝐻𝐻 , and all the no-slip boundary conditions for the azimuthal velocity, as well as the axis condition, are satisfied by each r.h.s. vector of (3.4.4) separately. Thus, considering flows with no-slip cylindrical boundaries, we can decompose 𝒗𝒗 𝑘𝑘≠0 into two independent bases, so that the whole flow will be represented as 𝒗𝒗 = ∑ ∑ 𝐴𝐴 𝑖𝑖𝑖𝑖 ( 𝑡𝑡 ) 𝑀𝑀 𝑧𝑧 𝑖𝑖=0 𝑈𝑈 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) 𝑀𝑀 𝑟𝑟 𝑖𝑖=0 + ∑ �∑ ∑ �𝐵𝐵 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) 𝑉𝑉 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) + 𝐶𝐶 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) 𝑊𝑊 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) � 𝑁𝑁 𝑧𝑧 𝑖𝑖=0𝑁𝑁 𝑟𝑟 𝑖𝑖=0 � 𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝑑𝑑𝑘𝑘𝜃𝜃 ) 𝑘𝑘=+∞𝑘𝑘=−∞𝑘𝑘≠0 (3.4.5) Here the functions 𝑈𝑈 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) represent the axisymmetric part of the flow and are defined in (3.4.3).The functions 𝑉𝑉 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) and 𝑊𝑊 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) represent the two vectors in r.h.s. of (3.4.4) and are defined as 𝑽𝑽 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) = ⎩⎪⎨⎪⎧−𝑑𝑑𝑘𝑘 � 𝑟𝑟𝑅𝑅 � 𝑞𝑞 ∑ 𝜎𝜎� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇̅ 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � ∑ 𝜎𝜎� 𝑖𝑖𝑚𝑚 𝑈𝑈� 𝑖𝑖+𝑚𝑚−1 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇̅ 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � ⎭⎪⎬⎪⎫ , (3.4.5) 𝑾𝑾 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) = ⎩⎪⎪⎨⎪⎪⎧ � 𝑟𝑟𝑅𝑅 � ∑ 𝜎𝜎� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇̿ 𝑖𝑖𝑚𝑚 𝑈𝑈 𝑖𝑖+𝑚𝑚−1 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � − 𝑖𝑖𝑘𝑘𝐴𝐴 𝑧𝑧 𝑟𝑟 ∑ 𝜎𝜎� 𝑖𝑖𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚 � 𝑟𝑟𝑅𝑅 � ∑ 𝜇𝜇� 𝑗𝑗𝑖𝑖 𝑖𝑖+𝑚𝑚 𝑇𝑇 𝑖𝑖+𝑚𝑚𝑚𝑚 � 𝑧𝑧𝐴𝐴 𝑧𝑧 � ⎭⎪⎪⎬⎪⎪⎫ , (3.4.6) Here 𝑞𝑞 = 0 for | 𝑘𝑘 | = 1 and 𝑞𝑞 = 1 for | 𝑘𝑘 | > 1 , 𝑈𝑈� 𝑛𝑛 ( 𝑥𝑥 ) = ( 𝑞𝑞 + 1) 𝑟𝑟 𝑞𝑞 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ) + 2 𝑠𝑠𝑟𝑟 𝑞𝑞+1 𝑈𝑈 𝑛𝑛−1 ( 𝑥𝑥 ) . Again, the zero divergence of the functions 𝑽𝑽 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) 𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝑑𝑑𝑘𝑘𝜃𝜃 ) and 𝑾𝑾 𝑖𝑖𝑖𝑖 ( 𝑟𝑟 , 𝑧𝑧 ) 𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝑑𝑑𝑘𝑘𝜃𝜃 ) follow from Eqs. (3.4.2) and (A2), and the coefficients 𝜎𝜎� 𝑖𝑖𝑚𝑚 , 𝜇𝜇̅ 𝑖𝑖𝑚𝑚 , 𝜎𝜎� 𝑖𝑖𝑚𝑚 , and 𝜇𝜇̿ 𝑖𝑖𝑚𝑚 are defined after substitution of (3.4.5) and (3.4.6) in the boundary conditions. The integer parameter q appears because of different boundary conditions posed for | 𝑘𝑘 | = 1 and | 𝑘𝑘 | ≠ at the polar axis (Canuto, et al., 2006; Gelfgat et al. 1999), which are At 𝑟𝑟 = 0 : 𝑢𝑢 𝑘𝑘=0 = 0, 𝑑𝑑 𝑘𝑘=0 = 0, 𝜕𝜕𝑤𝑤 𝑘𝑘=0 𝜕𝜕𝑟𝑟 = 0 𝑢𝑢 𝑘𝑘= ± ≠ 𝑑𝑑 𝑘𝑘= ± ≠ 𝑤𝑤 𝑘𝑘= ± = 0 (3.4.7) 𝑢𝑢 | 𝑘𝑘 | >1 = 0, 𝑑𝑑 | 𝑘𝑘 | >1 = 0, 𝑤𝑤 | 𝑘𝑘 | >1 = 0 For no-slip conditions at the top, bottom and sidewall of the cylinder the coefficients 𝜎𝜎 𝑖𝑖𝑚𝑚 , 𝜇𝜇 𝑖𝑖𝑚𝑚 , 𝜎𝜎� 𝑖𝑖𝑚𝑚 , 𝜇𝜇̅ 𝑖𝑖𝑚𝑚 , 𝜎𝜎� 𝑖𝑖𝑚𝑚 , and 𝜇𝜇̿ 𝑖𝑖𝑚𝑚 are defined as 𝜎𝜎 𝑖𝑖1 = − 𝑖𝑖 +7𝑖𝑖 +15𝑖𝑖+9𝑖𝑖 +6𝑖𝑖 +12𝑖𝑖+8 , 𝜎𝜎 𝑖𝑖2 = − 𝑖𝑖 ( 𝑖𝑖+2 ) , 𝜎𝜎 𝑖𝑖3 = 𝑖𝑖 +3𝑖𝑖 +3𝑖𝑖+1𝑖𝑖 +6𝑖𝑖 +12𝑖𝑖+8 , 𝜎𝜎 𝑖𝑖4 = 0 (3.4.8) 𝜎𝜎� 𝑖𝑖1 = − ( 𝑖𝑖+1 ) , 𝜎𝜎� 𝑖𝑖2 = , 𝜎𝜎� 𝑖𝑖3 = ( 𝑖𝑖+2 ) , 𝜎𝜎� 𝑖𝑖4 = 0 (3.4.9) 𝜎𝜎� 𝑖𝑖1 = − 𝜎𝜎� 𝑖𝑖2 = 0, 𝜎𝜎� 𝑖𝑖3 = 0, 𝜎𝜎� 𝑖𝑖4 = 0 (3.4.10) 𝜇𝜇 𝑖𝑖1 = 𝜇𝜇 𝑖𝑖3 = 𝜇𝜇̿ 𝑖𝑖1 = 𝜇𝜇̿ 𝑖𝑖3 = 0, 𝜇𝜇 = 𝜇𝜇̿ = − , 𝜇𝜇 = 𝜇𝜇̿ = ; (3.4.11) 𝜇𝜇 𝑖𝑖2 = 𝜇𝜇 � 𝑖𝑖2 = − 𝑖𝑖𝑖𝑖+2 − ( 𝑖𝑖+1 )( 𝑖𝑖+4 ) 𝑖𝑖 ( 𝑖𝑖+2 )( 𝑖𝑖+3 ) , 𝜇𝜇 𝑖𝑖4 = 𝜇𝜇 � 𝑖𝑖4 = ( 𝑖𝑖+1 )( 𝑖𝑖+4 ) 𝑖𝑖 ( 𝑖𝑖+3 ) , 𝑑𝑑 > 0 (3.4.12) 𝜇𝜇̅ 𝑖𝑖1 = 𝜇𝜇̅ 𝑖𝑖3 = 𝜇𝜇̅ 𝑖𝑖4 = 0, 𝜇𝜇̅ 𝑖𝑖2 = − (3.4.13) This example of divergence-free basis functions built for cylindrical geometries also shows how construction of a divergence free basis satisfying all the boundary conditions can be approached in other orthogonal coordinate systems. The process can be noticeably simplified if two periodic spatial coordinates are involved in the formulation of the problem. In these cases only a one-dimensional basis will be needed. Alternatively, the divergence-free basis in cylindrical and spherical coordinates with two periodic directions can be defined as in Dumas & Leonard (1994) or Meseguer & Melibovsky (2007).
4. Inhomogeneous boundary conditions
If the problem has inhomogeneous boundary conditions, they can be included as algebraic constraints. A better method is a change of variables so that all inhomogeneities are moved from the boundary conditions to the equations. Then the boundary conditions become homogeneous, and the corresponding basis functions can be built as in Section 3.4. We assume that all the boundary conditions are linear. Such change of variables can use analytic, as well as numerically calculated functions. Several examples of this are briefly described below. The simplest example for the change of variables is convection in a box, which has constant temperatures at the opposite sides, while all the other boundaries are perfectly thermally conducting or perfectly insulated. In the first case the temperature varies linearly along these boundaries, while in the second case temperature derivative normal o the boundary must vanish. These boundary conditions are satisfied by the linear temperature profile, which corresponds to the temperature distribution in a purely conducting case. For example, if for the dimensionless temperature 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) , the boundary conditions in the x- direction are 𝜃𝜃 ( 𝑥𝑥 = 0, 𝑦𝑦 , 𝑧𝑧 ) = 1, 𝜃𝜃 ( 𝑥𝑥 = 1, 𝑦𝑦 , 𝑧𝑧 ) = 0 , the change of variable is 𝜃𝜃 = (1 − 𝑥𝑥 ) + 𝜃𝜃� ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) . The function (1 − 𝑥𝑥 ) satisfies all homogeneous and inhomogeneous boundary conditions, so that all the boundary conditions for new unknown function 𝜃𝜃� ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) are homogeneous. This change of variables was applied in all cited works of Gelfgat that treated convection in rectangular cavities starting from Gelfgat & Tanasawa (1994). The inhomogeneities can be excluded from the boundary conditions by extracting a known analytical function from the solution only if the boundary conditions are continuous, including continuity at the corners of the computational domain. Another example of this is parabolic heating of a vertical cylinder that was considered in Gelfgat et al. (2000, 2001b). The dimensionless temperature at the cylindrical sidewall was prescribed as 𝜃𝜃 ( 𝑟𝑟 = 1, 𝑧𝑧 ) = 𝑓𝑓 ( 𝑧𝑧 ) = 𝑧𝑧 (1 − 𝑧𝑧 𝐴𝐴⁄ ), 0 ≤ 𝑧𝑧 ≤ 𝐴𝐴 , and was zero at the top and bottom, 𝑧𝑧 = 0, 𝐴𝐴 . Since the function 𝑓𝑓 ( 𝑧𝑧 ) vanishes at the top and the bottom, the simplest change of variables in this case is 𝜃𝜃 = 𝑓𝑓 ( 𝑧𝑧 ) + 𝜃𝜃� ( 𝑟𝑟 , 𝑧𝑧 ) , which was applied in the mentioned studies. Clearly, when the boundary conditions are discontinuous, use of a simple analytic function for the change of variables becomes problematic. Such a function can be built, for example, as a solution of Laplace equation with discontinuous boundary conditions. This analytic solution will suffer from Gibbs phenomenon which may destroy the convergence of the whole process. On the other hand, a low-order numerical solution can smooth the discontinuity up to an acceptable level, as in many calculations of lid-driven cavity flow, however this will take us too far from our Chebyshev-polynomial-based Galerkin approach. Thus, for calculation of the lid-driven cavity flow in Gelfgat (2005) we used analytically smoothed boundary condition, then solved the Stokes problem with the smoothed boundary conditions, and then used the Stokes problem solution for a change of variables. The Stokes problem was solved using the same Galerkin approach. In the studies of swirling flows in a cylinder with rotating lid, as well as independently rotating top, bottom and sidewall of the cylinder (Gelfgat et al., 1996a,b, 2001; Marques et al., 2003) we solved the Stokes problem for the azimuthal velocity component. A similar Galerkin method in scalar formulation was applied. The boundary conditions with discontinuities in the corners were included as additional algebraic constraints by adding collocation points along the boundaries. A more complicated case was treated in Erenburg et al. (2003). There we considered convection in a rectangular cavity with partially heated and partially insulated sidewall as sketched in Fig. 1. All the boundaries are no-slip. The dimensionless boundary conditions for the temperature are (here 𝐴𝐴 = 𝐻𝐻 𝐿𝐿 , ⁄ 𝑠𝑠 = ℎ 𝐿𝐿⁄ , 𝑠𝑠 = ℎ 𝐿𝐿⁄ ) 𝜃𝜃 ( 𝑥𝑥 ; 𝑦𝑦 = 0, 𝐴𝐴 ) = 0 , (4.1) 𝜃𝜃 ( 𝑥𝑥 = 0,1; 𝑠𝑠 ≤ 𝑦𝑦 ≤ 𝑠𝑠 ) = 1 , (4.2) 𝜕𝜕𝜕𝜕𝜕𝜕𝑥𝑥 ( 𝑥𝑥 = 0,1; 𝑦𝑦 < 𝑠𝑠 𝑜𝑜𝑟𝑟 𝑦𝑦 > 𝑠𝑠 ) = 0 . (4.3) To arrive at a formulation with continuous and homogeneous boundary conditions on all the boundaries for a single unknown function, we decompose the temperature into two functions 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) = Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) + 𝜃𝜃� ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) (4.4) where 𝜃𝜃� ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is the new unknown function for which a continuous set of boundary conditions is required, i.e., 𝜃𝜃� ( 𝑥𝑥 ; 𝑦𝑦 = 0, 𝐴𝐴 ) = 0 , 𝜕𝜕𝜕𝜕�𝜕𝜕𝑥𝑥 ( 𝑥𝑥 = 0,1; 𝑦𝑦 ) = 0 (4.5) x H y h h θ = θ cold = ∂∂ r θ =∂ ∂ r θ = ∂∂ r θ =∂ ∂ r θ θ = θ cold θ = θ hot θ = θ hot L Fig. 1. Sketch of the temperature boundary conditions of the problem of Erenburg et al. (2003). The function Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is used to adjust the boundary conditions for 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) to (4.1)-4.3). Therefore, the boundary conditions for Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) are Θ ( 𝑥𝑥 ; 𝑦𝑦 = 0, 𝐴𝐴 ) = 0 , (4.6) Θ ( 𝑥𝑥 = 0,1; 𝑠𝑠 ≤ 𝑦𝑦 ≤ 𝑠𝑠 ) = 1 − 𝜃𝜃� , (4.7) 𝜕𝜕Θ𝜕𝜕𝑥𝑥 ( 𝑥𝑥 = 0,1; 𝑦𝑦 < 𝑠𝑠 𝑜𝑜𝑟𝑟 𝑦𝑦 > 𝑠𝑠 ) = 0 (4.8) To avoid the appearance of an additional source term in the energy equation we also require that Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) be a solution of the Laplace equation, ∆Θ = 0 . (4.9) The solution of problem (4.6)-(4.9) can be represented as Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) = Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) + Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) (4.10) where Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is the part of the solution of (4.6)-(4.9) corresponding to 𝜃𝜃� = 0 and Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is the part dependent on 𝜃𝜃� . Both functions Θ and Θ are calculated numerically by the global Galerkin method inside the cavity and collocation points at the sidewalls. Obviously, the part Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is defined by the geometry of the problem only, and is time-independent, so that it must be calculated only once. The problem formulation for Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) is straight-forward, and its solution can be presented as Θ = 𝐿𝐿 −1 𝜃𝜃� , where 𝐿𝐿 is the operator defining the problem, and approximated by its Galerkin/collocation projection. The representation of the temperature (4.4) now becomes 𝜃𝜃 ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) = Θ ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) + 𝜃𝜃� ( 𝑥𝑥 , 𝑦𝑦 , 𝑡𝑡 ) = Θ + ( 𝐿𝐿 −1 + 𝐼𝐼 ) 𝜃𝜃� , (4.11) where 𝐼𝐼 is the identity operator. The energy equation becomes ( 𝐿𝐿 −1 + 𝐼𝐼 ) 𝜕𝜕𝜕𝜕�𝜕𝜕𝜕𝜕 + ( 𝒗𝒗 ∙ ∇ )( 𝐿𝐿 −1 + 𝐼𝐼 ) 𝜃𝜃� = ∆ ( 𝐿𝐿 −1 + 𝐼𝐼 ) 𝜃𝜃� − ( 𝒗𝒗 ∙ ∇ ) Θ . (4.12) Thus, after the function Θ and the operator ( 𝐿𝐿 −1 + 𝐼𝐼 ) are calculated, the remaining problem for 𝜃𝜃� is defined with the homogeneous boundary conditions (4.5) only. Other details can be found in Erenburg et al. (2003). Basis functions for two-fluid flow and boundary conditions at liquid-liquid interface
Here we give an example of basis functions that were used to calculate a two-phase flow with a capillary interface separating two liquids. A two-fluid Dean flow sketched in Fig. 2 was considered in Gelfgat et al. (2001d). The two fluids occupy adjacent thin cylindrical layers, 𝑠𝑠 ≤ 𝑟𝑟 ≤ 𝑠𝑠 + 𝑏𝑏 and 𝑠𝑠 + 𝑏𝑏 ≤ 𝑟𝑟 ≤ 𝑠𝑠 + 𝑑𝑑 , respectively, with the assumption 𝑠𝑠� = 𝑠𝑠 𝑑𝑑⁄ ≫ . The two fluids are separated by the border 𝑟𝑟 = 𝑠𝑠 + 𝑏𝑏 . The base flow in both fluids is driven by a constant azimuthal pressure gradient. Note that while the pressure is then a multi-valued function of 𝜃𝜃 , 𝜕𝜕𝑝𝑝 𝜕𝜕𝜃𝜃⁄ is not and can be considered as an imposed bulk force This formulation is an extension of the classical Dean (1928) problem to two-fluid system and is a simplified model of flow in a curved channel. Here we leave all the details concerning the evaluation of the base flow and the formulation of the stability problem to Gelfgat et al. (2001d), and focus only on the boundary conditions and incorporation of them into the basis functions. The three-dimensional velocity perturbation in cylindrical coordinates is defined as 𝒗𝒗 = �𝑢𝑢 ( 𝑥𝑥 ), 𝑑𝑑 ( 𝑥𝑥 ), 𝑤𝑤 ( 𝑥𝑥 ) �𝑅𝑅𝑥𝑥𝑝𝑝 ( 𝜆𝜆𝑡𝑡 + 𝑑𝑑𝑠𝑠𝜃𝜃 + 𝑑𝑑𝑘𝑘𝑧𝑧 ) . For convenience, we define a new dimensionless coordinate 𝑥𝑥 = ( 𝑟𝑟 − 𝑠𝑠 ) 𝑑𝑑⁄ , and define 𝑏𝑏� = ( 𝑏𝑏 − 𝑠𝑠 ) 𝑑𝑑⁄ . The dimensionless amplitude of the interface perturbation is 𝛿𝛿̅ . Indices 1 and 2 denote the variables in each sublayer, 𝜌𝜌 = 𝜌𝜌 𝜌𝜌 ⁄ and 𝜇𝜇 = 𝜇𝜇 𝜇𝜇 ⁄ is the ratio of densities and viscosities, respectively. After the axial velocity 𝑤𝑤 and the pressure 𝑝𝑝 are eliminated, the conditions for the radial and azimuthal components at the bounding surfaces and the interface read 𝑥𝑥 = 0: 𝑢𝑢 = 𝑑𝑑 = 𝑑𝑑𝜕𝜕 𝑑𝑑𝑥𝑥 = 0 (5.1) 𝑥𝑥 = 1: 𝑢𝑢 = 𝑑𝑑 = 𝑑𝑑𝜕𝜕 𝑑𝑑𝑥𝑥 = 0 (5.2) 𝑥𝑥 = 𝑏𝑏� : 𝑢𝑢 = 𝑢𝑢 , 𝑑𝑑 = 𝑑𝑑 (5.3) 𝑑𝑑𝜕𝜕 𝑑𝑑𝑥𝑥 = 𝑑𝑑𝜕𝜕 𝑑𝑑𝑥𝑥 (5.4) 𝑑𝑑𝑣𝑣 𝑑𝑑𝑥𝑥 = 𝜇𝜇
12 𝑑𝑑𝑣𝑣 𝑑𝑑𝑥𝑥 (5.5) 𝑑𝑑 𝜕𝜕 𝑑𝑑𝑥𝑥 + 𝑘𝑘 𝑢𝑢 = 𝜇𝜇 � 𝑑𝑑 𝜕𝜕 𝑑𝑑𝑥𝑥 + 𝑘𝑘 𝑢𝑢 � (5.6) 𝜆𝜆 𝛿𝛿 � = 𝑢𝑢 − 𝑑𝑑𝑠𝑠 𝑏𝑏� 𝑉𝑉𝛿𝛿 � (5.7) 𝜆𝜆 � 𝜌𝜌 𝑑𝑑𝑢𝑢 𝑑𝑑𝑥𝑥 − 𝑑𝑑𝑢𝑢 𝑑𝑑𝑥𝑥 � = 𝛿𝛿 � 𝑠𝑠 � 𝑘𝑘 � 𝜌𝜌 − � 𝑉𝑉 + + � 𝑑𝑑 𝑑𝑑𝑥𝑥 − 𝑑𝑑𝑑𝑑𝑥𝑥 − 𝑘𝑘 � ( 𝜇𝜇 𝑢𝑢 − 𝑢𝑢 ) − 𝜌𝜌
12 𝑖𝑖𝑛𝑛𝑎𝑎� 𝑉𝑉 𝑑𝑑𝑑𝑑𝑥𝑥 ( 𝜌𝜌 𝑢𝑢 − 𝑢𝑢 ) − 𝑘𝑘 𝑊𝑊𝑅𝑅 � 𝑏𝑏� − 𝑘𝑘 � 𝛿𝛿̅ (5.8) To construct basis functions, we start from non-deformable interface. In this case we add 𝑢𝑢 = 𝑢𝑢 = 0 to the boundary condition (5.3) and omit (5.7) and (5.8). Then the unknowns 𝑢𝑢 ( 𝑥𝑥 ) and 𝑑𝑑 ( 𝑥𝑥 ) are approximated by truncated series Fig. 2. Sketch of the two-fluid Dean flow problem b a r θ d V θ 𝑢𝑢 ( 𝑥𝑥 ) = ∑ 𝑐𝑐 𝑚𝑚 ( 𝑡𝑡 ) 𝑁𝑁𝑚𝑚=1 𝜓𝜓 𝑚𝑚 ( 𝑥𝑥 ), 𝑑𝑑 ( 𝑥𝑥 ) = ∑ 𝑑𝑑 𝑚𝑚 ( 𝑡𝑡 ) 𝑁𝑁𝑚𝑚=1 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) , (5.9) where 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) = ⎩⎪⎨⎪⎧∑ 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) 𝑇𝑇 𝑚𝑚+𝑚𝑚 � 𝑥𝑥𝑏𝑏� � , 0 ≤ 𝑥𝑥 ≤ 𝑏𝑏�∑ 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) 𝑇𝑇 𝑚𝑚+𝑚𝑚 � 𝑥𝑥−𝑏𝑏�1−𝑏𝑏� � , 𝑏𝑏� ≤ 𝑥𝑥 ≤ , (5.10) 𝜓𝜓 𝑚𝑚 ( 𝑥𝑥 ) = ⎩⎪⎨⎪⎧∑ 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) 𝑇𝑇 𝑚𝑚+𝑚𝑚 � 𝑥𝑥𝑏𝑏� � , 0 ≤ 𝑥𝑥 ≤ 𝑏𝑏�∑ 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) 𝑇𝑇 𝑚𝑚+𝑚𝑚 � 𝑥𝑥−𝑏𝑏�1−𝑏𝑏� � , 𝑏𝑏� ≤ 𝑥𝑥 ≤ . (5.11) Here the superscripts (1) and (2) denote the sublayers. The coefficients 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) , 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) , 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) , and 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) are obtained after substitution of the basis functions (5.11) into the boundary conditions (5.1)-(5.6). The inner product is defined as 〈𝑓𝑓 , 𝑔𝑔〉 = ∫ 𝑓𝑓 ( 𝑥𝑥 ) 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 = ∫ 𝑓𝑓 ( 𝑥𝑥 ) 𝑏𝑏�0 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 + ∫ 𝑓𝑓 ( 𝑥𝑥 ) 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 , (5.12) and after application of the Galerkin method, the time-dependent coefficients 𝑐𝑐 𝑚𝑚 ( 𝑡𝑡 ) and 𝑑𝑑 𝑚𝑚 ( 𝑡𝑡 ) are the same for the whole domain. For the linear stability problem with deformable interface the solution is represented as 𝑢𝑢 ( 𝑥𝑥 ) = 𝑐𝑐 ( 𝑡𝑡 ) 𝜙𝜙 ( 𝑥𝑥 ) + ∑ 𝑐𝑐 𝑚𝑚 ( 𝑡𝑡 ) 𝑁𝑁𝑚𝑚=1 𝜓𝜓 𝑚𝑚 ( 𝑥𝑥 ), 𝑑𝑑 ( 𝑥𝑥 ) = ∑ 𝑑𝑑 𝑚𝑚 ( 𝑡𝑡 ) 𝑁𝑁𝑚𝑚=1 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) . (5.13) The bases 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) and 𝜓𝜓 𝑚𝑚 ( 𝑥𝑥 ) remain unchanged. An additional function 𝜙𝜙 ( 𝑥𝑥 ) is introduced to satisfy the boundary conditions for a deformable interface. Its choice is arbitrary. In our case we define it as 𝜙𝜙 ( 𝑥𝑥 ) = ⎩⎪⎨⎪⎧∑ 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) 𝑇𝑇 𝑞𝑞+𝑚𝑚 � 𝑥𝑥𝑏𝑏� � , 0 ≤ 𝑥𝑥 ≤ 𝑏𝑏�∑ 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) 𝑇𝑇 𝑞𝑞+𝑚𝑚 � 𝑥𝑥−𝑏𝑏�1−𝑏𝑏� � , 𝑏𝑏� ≤ 𝑥𝑥 ≤ . (5.14) The coefficients 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) and 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) are used to satisfy the boundary conditions (5.1), (5.2), (5.4) and (5.6) subject to the normalization condition 𝜙𝜙�𝑏𝑏�� = 1 . Choice of the index 𝑞𝑞 is arbitrary, e.g., 𝑞𝑞 = 0 . With the normalization condition 𝜙𝜙�𝑏𝑏�� = 1 applied, the coefficient 𝑐𝑐 can be interpreted as the amplitude of the radial velocity at the deformed interface. This coefficient, and the interface amplitude 𝛿𝛿̅ , are defined by the two remaining boundary conditions (5.7) and (5.8). Thus, the Galerkin projections of the governing equations together with the boundary conditions (5.7) and (5.8) form a closed algebraic system for calculation of the coefficients 𝑐𝑐 𝑚𝑚 and 𝑑𝑑 𝑚𝑚 together with two additional unknown scalars 𝑑𝑑 𝑚𝑚 and 𝛿𝛿̅ . The stability problem reduces to a generalized eigenvalue problem. In Gelfgat et al. (2001d) coefficients of the basis functions 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) , 𝛼𝛼 𝑚𝑚𝑚𝑚 ( ) , 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) , 𝛽𝛽 𝑚𝑚𝑚𝑚 ( ) , 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) , and 𝛾𝛾 𝑞𝑞𝑚𝑚 ( ) were obtained by means of computer algebra.
6. The dynamical ODEs system for time-dependent coefficients
In the following we assume that all the inner products are defined with unit weight. We also assume that all the necessary changes of variables are made, so that the boundary conditions of all the unknown vector and scalar fields are linear and homogeneous. Then, after the basis functions are constructed, and the Galerkin projections are made, and the pressure is eliminated by Eq. (2.20), we arrive at an ODE system (2.21) for the time-dependent coefficients. We store all the unknown time-dependent coefficients of the problem in the vector 𝑿𝑿 ( 𝑡𝑡 ) = { 𝑋𝑋 𝐼𝐼 ( 𝑡𝑡 )} 𝐼𝐼=1𝑁𝑁 , where 𝑁𝑁 is the total number of scalar unknowns (degrees of freedom). Note that the vector 𝑿𝑿 ( 𝑡𝑡 ) contains velocity coefficients, as well as coefficients of all the other possible unknowns, e.g., temperature and/or concentration, excluding pressure. Then the resulting dynamic ODE system has the following form (the Einstein summation rule is assumed) 𝐺𝐺 𝐼𝐼𝐽𝐽
𝑋𝑋̇ 𝐽𝐽 = 𝐿𝐿 𝐼𝐼𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹 𝐼𝐼 . (6.1.1) Here 𝐺𝐺 𝐼𝐼𝐽𝐽 is the Gram matrix, 𝐿𝐿 𝐼𝐼𝐽𝐽 , 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 , and 𝐹𝐹 𝐼𝐼 are projections of the linear, nonlinear and bulk force of the momentum equation (2.7), respectively. The transport equation for temperature, concentration, electric and magnetic fields, etc., arrive to similar ODE systems that can be complicated by non-linear terms not belonging to the material derivative. This dynamical system has several nice properties that follow from the fact that the basis functions satisfy all boundary conditions, and are divergence-free. From Green’s theorems for a scalar function and for a divergence-free velocity 〈∆𝜃𝜃 , 𝜃𝜃〉 = −〈∇𝜃𝜃 , ∇𝜃𝜃〉 , 〈∆𝒗𝒗 , 𝒗𝒗〉 = −〈∇ × 𝒗𝒗 , ∇ × 𝒗𝒗〉 . (6.1.2) It follows that the matrices 𝐿𝐿 𝐼𝐼𝐽𝐽 corresponding to the dissipative terms are symmetric and negative definite independently on the truncation number. Furthermore, from the conservation properties 〈 ( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 , 𝒗𝒗〉 = 0 , 〈 ( 𝒗𝒗 ∙ ∇ ) 𝜃𝜃 , 𝜃𝜃〉 = 0 (6.1.3) it follows that 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐼𝐼 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 = 0 (6.1.4) for any truncation number. This means that the non-linear term always conserves the momentum. Additionally, these properties yield a very convenient tool for code debugging. Computation of steady state flows reduces to an algebraic system of quadratic equations 𝐿𝐿 𝐼𝐼𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹 𝐼𝐼 = 0 , (6.1.5) for which we do not need to consider the Gram matrix. The application of Newton iteration is straightforward and requires computation and inversion of the Jacobian matrix ℑ 𝐼𝐼𝐽𝐽 = 𝐿𝐿 𝐼𝐼𝐽𝐽 + �𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 + 𝑁𝑁 𝐼𝐼𝐾𝐾𝐽𝐽 �𝑋𝑋 𝐾𝐾 . (6.1.6) Linear stability analysis of the calculated steady states reduces to computation of the eigenvalues of another Jacobian matrix, which includes the inverted Gram matrix ℑ�𝒀𝒀 = 𝜆𝜆𝒀𝒀 , ℑ� 𝐼𝐼𝐽𝐽 = 𝐺𝐺 𝐼𝐼𝑀𝑀−1 �𝐿𝐿
𝑀𝑀𝐽𝐽 + �𝑁𝑁 𝑀𝑀𝐽𝐽𝐾𝐾 + 𝑁𝑁 𝑀𝑀𝐾𝐾𝐽𝐽 �𝑋𝑋 𝐾𝐾 � = 𝐿𝐿�
𝐼𝐼𝐽𝐽 + �𝑁𝑁� 𝐼𝐼𝐽𝐽𝐾𝐾 + 𝑁𝑁�
𝐼𝐼𝐾𝐾𝐽𝐽 �𝑋𝑋 𝐾𝐾 . (6.1.7) The inverse of the Gram matrix is also needed for straightforward time integration, for which the dynamical system (6.1) has the form 𝑋𝑋̇ 𝐼𝐼 = 𝐺𝐺 𝐼𝐼𝑀𝑀−1 �𝐿𝐿
𝑀𝑀𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁 𝑀𝑀𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹 𝑀𝑀 � = 𝐿𝐿�
𝐼𝐼𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁�
𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹� 𝐼𝐼 , (6.1.8) where matrices multiplied by 𝐺𝐺 −1 are denoted by � . Explicit representation of the dynamical system (6.1.8) allows one to perform additional analytical evaluations required by weakly non-linear analysis of bifurcations. This was implemented for the Hopf bifurcation in Gelfgat et al. (1996) and Gelfgat (2004). Assume that with increasing Reynolds number, at 𝑅𝑅𝑅𝑅 = 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 , a complex conjugate pair Λ = ± 𝑑𝑑𝑖𝑖 of leading eigenvalues of the problem (6.1.7) crosses the imaginary axis. Then the instability sets in as a Hopf bifurcation if all the conditions of the Hopf theorem hold, which is the most common case. We are interested in an asymptotic approximation of the oscillation period and amplitude at small super-criticality. Assume that 𝑿𝑿 is the steady state at the critical point, and 𝑼𝑼 and 𝑽𝑽 are the left and right eigenvectors corresponding to the eigenvalue Λ = 𝑑𝑑𝑖𝑖 . We also denote the r.h.s. of the dynamic system (6.8) as 𝑭𝑭 ( 𝑿𝑿 ; 𝑅𝑅𝑅𝑅 ) =
𝑿𝑿̇ . Then, according to Hassard et al. (1981), the oscillating state, i.e., the limit cycle, is approximated as 𝑅𝑅𝑅𝑅 = 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 + 𝜇𝜇 𝜀𝜀 + 𝑂𝑂 ( 𝜀𝜀 ) (6.1.9) 𝜏𝜏 ( 𝑅𝑅𝑅𝑅 ) = [1 + 𝜏𝜏 𝜀𝜀 + 𝑂𝑂 ( 𝜀𝜀 )] (6.1.10) 𝑿𝑿 ( 𝑡𝑡 , 𝑅𝑅𝑅𝑅 ) = 𝑋𝑋 ( 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 ) + 𝜀𝜀𝑅𝑅𝑅𝑅𝑠𝑠𝑙𝑙 �𝑽𝑽𝑅𝑅𝑥𝑥𝑝𝑝 � �� + 𝑂𝑂 ( 𝜀𝜀 ) (6.1.11) Here 𝜀𝜀 is a formal positive parameter, 𝑅𝑅𝑅𝑅 − 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 is the super-criticality, 𝜏𝜏 is the oscillation period, and 𝑿𝑿 ( 𝑡𝑡 , 𝑅𝑅𝑅𝑅 ) is the asymptotic oscillatory solution of the ODE system (6.1.8) for the Reynolds number defined in (6.1.9). This asymptotic expansion is defined by two parameters 𝜇𝜇 and 𝜏𝜏 , which are calculated using the following procedu (Hassard et al., 1981) 𝜇𝜇 = − 𝑅𝑅𝑅𝑅𝑎𝑎𝑚𝑚 ( 𝜎𝜎 ) 𝛼𝛼 𝑟𝑟 , 𝜏𝜏 = − [ 𝐼𝐼𝑚𝑚 ( 𝜎𝜎 ) + 𝜇𝜇 𝛼𝛼 𝑖𝑖 ] , 𝛼𝛼 = 𝛼𝛼 𝑟𝑟 + 𝑑𝑑𝛼𝛼 𝑖𝑖 = � dΛ𝑑𝑑𝑅𝑅𝑅𝑅 � 𝑅𝑅𝑅𝑅=𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 , (6.1.12) The parameter 𝜎𝜎 is obtained as 𝜎𝜎 = 𝐻𝐻 + �𝑔𝑔 𝑔𝑔 − 𝑔𝑔 | − | 𝑔𝑔 | � (6.1.13) 𝑔𝑔 = 2 𝑼𝑼 𝑇𝑇 𝒇𝒇 , 𝑔𝑔 = 2 𝑼𝑼 𝑇𝑇 𝒇𝒇� , 𝑔𝑔 = 2 𝑼𝑼 𝑇𝑇 𝒇𝒇 . (6.1.14) The vectors 𝒇𝒇 and 𝒇𝒇 are the second derivatives of the r.h.s., and 𝐻𝐻 is the third derivative in the complex plane 𝐶𝐶 , 𝜉𝜉 ∈ 𝐶𝐶 : 𝒇𝒇 = � 𝜕𝜕 𝜕𝜕𝜉𝜉 𝑭𝑭𝜉𝜉 [ 𝑋𝑋 + 𝑅𝑅𝑅𝑅𝑠𝑠𝑙𝑙 ( 𝑽𝑽𝜉𝜉 ); 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 ] � 𝜉𝜉=0 , 𝒇𝒇 = � 𝜕𝜕 𝜕𝜕𝜉𝜉𝜕𝜕𝜉𝜉� 𝑭𝑭 [ 𝑋𝑋 + 𝑅𝑅𝑅𝑅𝑠𝑠𝑙𝑙 ( 𝑽𝑽𝜉𝜉 ); 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 ] � 𝜉𝜉=0 (6.1.15) 𝐻𝐻 = � 𝜕𝜕 𝜕𝜕𝜕𝜕𝜉𝜉� � 𝑼𝑼 𝑇𝑇 𝑭𝑭�𝑋𝑋 + 𝑅𝑅𝑅𝑅𝑠𝑠𝑙𝑙�𝑽𝑽𝜉𝜉 + 𝝔𝝔 𝜉𝜉 + 𝝔𝝔 𝜉𝜉𝜉𝜉̅� ; 𝑅𝑅𝑅𝑅 𝑐𝑐𝑟𝑟 ��� 𝜉𝜉=0 , (6.1.16) and the vectors 𝝔𝝔 and 𝝔𝝔 are solutions of the following systems of linear algebraic equations ( 𝐼𝐼 is the identity matrix and ⨂ denotes the Kronecker product) ℑ�𝝔𝝔 = −𝒉𝒉 , �ℑ� − 𝑑𝑑𝑖𝑖 𝐼𝐼�𝝔𝝔 = −𝒉𝒉 , 𝒉𝒉 𝑖𝑖𝑖𝑖 = �𝐼𝐼 − 𝑅𝑅𝑅𝑅𝑠𝑠𝑙𝑙 ( 𝑉𝑉⨂𝑈𝑈 𝑇𝑇 ) �𝒇𝒇 𝑖𝑖𝑖𝑖 (6.1.17) The eigenvalue derivative (6.1.12) can be easily evaluated numerically. However, the most difficult part of this calculation is evaluation of the second and the third derivatives of the right hand side of the ODE system (6.1.15) and (6.1.16). These derivatives must be evaluated in the complex plane. The number of degrees of freedom is large, so that accurate enough differentiation by finite differences that will involve small increments of 𝜉𝜉 can be problematic. However, the explicit form of (6.1.8) allows for analytical calculation of the derivatives. The result is 𝑓𝑓 , 𝑘𝑘 = �𝑁𝑁 𝐾𝐾𝐼𝐼𝐽𝐽 �𝑉𝑉 𝐼𝐼 ( 𝑟𝑟 ) 𝑉𝑉 𝐽𝐽 ( 𝑟𝑟 ) − 𝑉𝑉 ( 𝑖𝑖 ) 𝑉𝑉 𝐽𝐽 ( 𝑖𝑖 ) � + 𝑑𝑑𝑁𝑁 𝐾𝐾𝐼𝐼𝐽𝐽 �𝑉𝑉 𝐼𝐼 ( 𝑟𝑟 ) 𝑉𝑉 𝐽𝐽 ( 𝑖𝑖 ) − 𝑉𝑉 𝐼𝐼 ( 𝑖𝑖 ) 𝑉𝑉 𝐽𝐽 ( 𝑟𝑟 ) �� , (6.1.18) 𝑓𝑓 , 𝑘𝑘 = 𝑁𝑁 𝑘𝑘𝑖𝑖𝑖𝑖 �𝑉𝑉 𝐼𝐼 ( 𝑟𝑟 ) 𝑉𝑉 𝐽𝐽 ( 𝑟𝑟 ) + 𝑉𝑉 ( 𝑖𝑖 ) 𝑉𝑉 𝐽𝐽 ( 𝑖𝑖 ) � , (6.1.19) 𝐺𝐺 = 𝑈𝑈 𝑖𝑖 �𝑁𝑁 𝑖𝑖𝑖𝑖𝑚𝑚 + 𝑁𝑁 𝑖𝑖𝑚𝑚𝑖𝑖 �� 𝝔𝝔 , 𝑖𝑖 𝑉𝑉 𝑚𝑚 + 2 𝝔𝝔 , 𝑖𝑖 𝑉𝑉� 𝑚𝑚 � , (6.1.20) where superscripts ( 𝑟𝑟 ) and ( 𝑑𝑑 ) denote real and imaginary parts, respectively. These analytical expressions allow one to compute the asymptotic expansions (6.1.9)-(6.1.11) without significant loss of accuracy. The CPU-time requirements for such calculations are comparable with the calculation of two steady state solutions and their spectra. The sign of 𝜇𝜇 determines whether the bifurcation is sub- or super-critical. All the computations described in the previous section, are restricted by two main difficulties. To describe these, we note that in all the studies where this method was successfully applied, the number of basis functions in one direction was between 20 and 70. Therefore, to make some estimates, we will address three types of truncation in one direction with 30, 60 and 100 basis functions for both two- and three-dimensional problems. The first difficulty is connected with the Gram matrix 𝐺𝐺 . In two-dimensional and quasi two-dimensional cases, e.g., 3D instability of an axisymmetric base flow, there is no problem storing the Gram matrix in memory, nor in computing its inverse. In fact, even with the truncation number 100, the Gram matrix consists of blocks of order 100 =10 , which leads to an order of 10 non-zero entries, which fits in several Gb memory. This matrix, which is symmetric and positive definite, can be inverted by Choleski decomposition, which is much faster than Gaussian elimination. This inverse must be computed only once for each specific problem and can be stored on disk. However, for all the calculations, except the Newton iterations, the r.h.s. of dynamic system must be multiplied by 𝐺𝐺 −1 . If the truncation number is small, the matrices 𝐿𝐿� = 𝐺𝐺 −1 𝐿𝐿 and 𝑁𝑁� = 𝐺𝐺 −1 𝑁𝑁 in the dynamical system (6.1.7) can be computed and stored before other heavy computations begin. At larger truncation number the storage of matrix 𝑁𝑁� becomes impossible (see below), which makes the time-dependent calculations too slow. At the same time, the stability analysis, as well as the weakly non-linear analysis, require only a few, usually less than 10, evaluations of the Jacobian matrix ℑ� , thus making the computational process affordable. Treatment of the Gram matrix in a three-dimensional formulation is much more difficult. Here the largest block to be inverted has order of ((2 𝑀𝑀 ) ) elements, where M is the truncation number. Storage of such large matrices becomes problematic already at truncation numbers 𝑀𝑀 ≥ , and is unaffordable at 𝑀𝑀 = 100 . Thus, among all the possibilities described, only steady state calculations are possible. A solution for this can be orthonormalization of the set of the basis functions, which is discussed below. The second difficulty is the numerical evaluation of non-linear term in (6.1.5) and (6.1.7), which requires the order of ((2 𝑀𝑀 ) ) multiplications, assuming that the matrix 𝑁𝑁� is stored and evaluation of its terms does not require additional operations. Again, it can be affordable for the steady state, stability, and weakly non-linear calculations when 2D and quasi-2D problems are solved, because all these require very few evaluations of the r.h.s. In the fully 3D cases it seems to be not feasible, unless some additional evaluations of the non-linear terms are performed.
As a simplest example of handling the non-linear terms and avoiding too many multiplications, we consider the Burgers equation in the interval ≤ 𝑥𝑥 ≤ 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 + 𝑢𝑢 𝜕𝜕𝜕𝜕𝜕𝜕𝑥𝑥 = 𝜈𝜈 𝜕𝜕 𝜕𝜕𝜕𝜕𝑥𝑥 , 𝑢𝑢 (0, 𝑡𝑡 ) = 𝑢𝑢 (1, 𝑡𝑡 ) = 0 . (6.3.1) The initial condition used in Gelfgat (2006) was 𝑢𝑢 ( 𝑥𝑥 , 0) = 𝑠𝑠𝑑𝑑𝑠𝑠 (2 𝜋𝜋𝑥𝑥 ) + 𝑠𝑠𝑑𝑑𝑠𝑠 ( 𝜋𝜋𝑥𝑥 ) 2 ⁄ . We seek the solution as a truncated series 𝑢𝑢 ( 𝑥𝑥 , 𝑡𝑡 ) = ∑ 𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) 𝑀𝑀−1𝑖𝑖=0 , 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) = 𝑇𝑇 𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑖𝑖+2 ( 𝑥𝑥 ) , (6.3.2) where the basis functions 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) are constructed as linear superpositions of the Chebyshev polynomials to satisfy the boundary conditions (see Eqs. (A3)). After the Galerkin projections are applied we obtain the ODE system (6.1.1), where coefficients 𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) are stored in the vector 𝑿𝑿 . For this problem the matrices in (6.1.1) are defined as 𝐺𝐺 𝑖𝑖𝑖𝑖 = 〈𝜑𝜑 𝑖𝑖 , 𝜑𝜑 𝑖𝑖 〉 , 𝐿𝐿 𝑖𝑖𝑖𝑖 = 〈𝜑𝜑 𝑖𝑖′′ , 𝜑𝜑 𝑖𝑖 〉 , 𝐹𝐹 𝑖𝑖 = 0 , 𝑁𝑁 𝑖𝑖𝑖𝑖𝑘𝑘 = 〈𝜑𝜑 𝑖𝑖 𝜑𝜑 𝑘𝑘′ , 𝜑𝜑 𝑖𝑖 〉 , (6.3.3) and evaluation of the non-linear term requires 𝑀𝑀 operations. To reduce the number of multiplications, we notice that 𝜑𝜑 𝑖𝑖 𝜑𝜑 𝑘𝑘′ is a polynomial of order 𝑗𝑗 + 𝑘𝑘 + 1 that satisfies the boundary conditions of the problem, so that we can express it as a series of 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) : 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) 𝜑𝜑 𝑘𝑘′ ( 𝑥𝑥 ) = ∑ 𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘𝑖𝑖+𝑘𝑘−1𝑚𝑚=0 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) (6.3.4) The coefficients 𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘 can be evaluated analytically using Eqs. (A7) and (A12) in the following way 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) 𝜑𝜑 𝑘𝑘′ ( 𝑥𝑥 ) = �𝑇𝑇 𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑖𝑖+2 ( 𝑥𝑥 ) � [ 𝑇𝑇 𝑘𝑘′ ( 𝑥𝑥 ) − 𝑇𝑇 𝑘𝑘+2′ ( 𝑥𝑥 )] = (6.3.5) = �𝑇𝑇 𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑖𝑖+2 ( 𝑥𝑥 ) � � 𝑘𝑘 � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝 ( 𝑥𝑥 ) [( 𝑘𝑘−1 ) ] 𝑝𝑝=0 − 𝑘𝑘 + 2) � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝 ( 𝑥𝑥 ) [( 𝑘𝑘+1 ) ] 𝑝𝑝=0 � = 2 𝑘𝑘 � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 ��𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖 ( 𝑥𝑥 ) + 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � − �𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 ) + 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖+2 ( 𝑥𝑥 ) �� [( 𝑘𝑘−1 ) ] 𝑝𝑝=0 − 𝑘𝑘 + 2) � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 ��𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖 ( 𝑥𝑥 )+ 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � − �𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 )+ 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖+2 ( 𝑥𝑥 ) �� [( 𝑘𝑘+1 ) ] 𝑝𝑝=0 = 2 𝑘𝑘 � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 ) + 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖+2 ( 𝑥𝑥 ) � [( 𝑘𝑘−1 ) ] 𝑝𝑝=0 − 𝑘𝑘 + 2) � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 )+ 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) − 𝑇𝑇 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖+2 ( 𝑥𝑥 ) � [( 𝑘𝑘+1 ) ] 𝑝𝑝=0 = 2 𝑘𝑘 � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �−𝜑𝜑 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 ) + 𝜑𝜑 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � [( 𝑘𝑘−1 ) ] 𝑝𝑝=0 − 𝑘𝑘 + 2) � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �−𝜑𝜑 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 )+ 𝜑𝜑 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � [( 𝑘𝑘+1 ) ] 𝑝𝑝=0 = − � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �−𝜑𝜑 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 ) + 𝜑𝜑 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � [( 𝑘𝑘−1 ) ] 𝑝𝑝=0 − 𝑘𝑘 + 2) � 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 �−𝜑𝜑 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 ( 𝑥𝑥 )+ 𝜑𝜑 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 ( 𝑥𝑥 ) � [( 𝑘𝑘+1 ) ] 𝑝𝑝= [( 𝑘𝑘−1 ) ] +1 Defining additionally 𝜑𝜑 𝑘𝑘<0 = 0 , noticing that [( 𝑘𝑘 −
1) 2 ⁄ ] + 1 = [( 𝑘𝑘 + 1) 2 ⁄ ] , and comparing the above result with (6.3.4) we observe that the coefficients 𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘 can be assembled by the following procedure. Starting from 𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘 = 0 , 𝑏𝑏 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] −𝑖𝑖−2 = 𝑏𝑏 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] −𝑖𝑖−2 + 2( 𝑘𝑘 + 2) 𝑠𝑠 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] , (6.3.6) 𝑑𝑑𝑓𝑓 𝑘𝑘 − − 𝑘𝑘 + 1) 2 ⁄ ] − 𝑗𝑗 − ≥ 𝑏𝑏 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] +𝑖𝑖 = 𝑏𝑏 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] +𝑖𝑖 − 𝑘𝑘 + 2) 𝑠𝑠 𝑘𝑘−1−2 [( 𝑘𝑘+1 ) ] , 𝑑𝑑𝑓𝑓 𝑘𝑘 − − 𝑘𝑘 + 1) 2 ⁄ ] + 𝑗𝑗 ≥ For 𝑝𝑝 = 0 to 𝑝𝑝 = [( 𝑘𝑘 −
1) 2 ⁄ ]: 𝑏𝑏 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 = 𝑏𝑏 𝑘𝑘−1−2𝑝𝑝−𝑖𝑖−2 + 4 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 , 𝑑𝑑𝑓𝑓 𝑘𝑘 − − 𝑝𝑝 − 𝑗𝑗 − ≥ 𝑏𝑏 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 = 𝑏𝑏 𝑘𝑘−1−2𝑝𝑝+𝑖𝑖 − 𝑠𝑠 𝑘𝑘−1−2𝑝𝑝 , 𝑑𝑑𝑓𝑓 𝑘𝑘 − − 𝑝𝑝 + 𝑗𝑗 ≥ Now, using (6.3.4), we form a new set of time-dependent coefficients: ∑ 𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) 𝑐𝑐 𝑘𝑘 ( 𝑡𝑡 ) 𝑀𝑀−1𝑖𝑖 , 𝑘𝑘=0 𝜑𝜑 𝑖𝑖 ( 𝑥𝑥 ) 𝜑𝜑 𝑘𝑘′ ( 𝑥𝑥 ) = ∑ 𝐶𝐶 𝑚𝑚 ( 𝑡𝑡 ) ( 𝑀𝑀−1 ) 𝑚𝑚=0 𝜑𝜑 𝑚𝑚 ( 𝑥𝑥 ) , (6.3.7) 𝐶𝐶 𝑚𝑚 ( 𝑡𝑡 ) = 𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) 𝑐𝑐 𝑘𝑘 ( 𝑡𝑡 ) ∑ �𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘 + 𝑏𝑏 𝑚𝑚𝑘𝑘𝑖𝑖 � 𝑖𝑖+𝑘𝑘−1𝑚𝑚=0 , 𝑚𝑚 = 𝑗𝑗 + 𝑘𝑘 (6.3.8) And finally, ∑ 𝑁𝑁 𝑖𝑖𝑖𝑖𝑘𝑘 𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) 𝑐𝑐 𝑘𝑘 ( 𝑡𝑡 ) 𝑀𝑀−1𝑖𝑖 , 𝑘𝑘=0 = ∑ 〈𝜑𝜑 𝑖𝑖 𝜑𝜑 𝑘𝑘′ , 𝜑𝜑 𝑖𝑖 〉𝑐𝑐 𝑖𝑖 ( 𝑡𝑡 ) 𝑐𝑐 𝑘𝑘 ( 𝑡𝑡 ) 𝑀𝑀−1𝑖𝑖 , 𝑘𝑘=0 = ∑ 𝐶𝐶 𝑚𝑚 ( 𝑡𝑡 ) ( 𝑀𝑀−1 ) 𝑚𝑚=0 〈𝜑𝜑 𝑚𝑚 , 𝜑𝜑 𝑖𝑖 〉 (6.3.9) Now, we can estimate the number of multiplications required. The coefficients 𝑏𝑏 𝑚𝑚𝑖𝑖𝑘𝑘 , and the sums in (6.3.8) depend only on the basis functions and, therefore, can be computed only once in the beginning of the computational process. Computation of the coefficients 𝐶𝐶 𝑚𝑚 ( 𝑡𝑡 ) requires 𝑀𝑀 multiplications and is the most CPU-time consuming part. Then, evaluation of (6.3.9) requires 𝑀𝑀 − multiplications providing that all the inner products are pre-computed. Since the operations in (6.3.8) and (6.3.9) are easily scalable, vectorization and/or parallel computing can additionally speed up the calculations. Returning to the non-linear terms of momentum equation, we observe that in the case of no-slip conditions the non-linear terms 𝑢𝑢 𝜕𝜕𝑢𝑢 𝜕𝜕𝑥𝑥⁄ , 𝑑𝑑 𝜕𝜕𝑢𝑢 𝜕𝜕𝑦𝑦⁄ , 𝑤𝑤𝜕𝜕𝑢𝑢 𝜕𝜕𝑧𝑧⁄ , etc., satisfy the no-slip boundary conditions for 𝑢𝑢 , 𝑑𝑑 , and 𝑤𝑤 , respectively. Thus, these terms can also be decomposed into series of appropriate basis functions, which will lead to a similar decrease in the number of required multiplications. When boundary conditions are more complicated, e.g., a stress-free boundary, one can add additional functions in which only boundary conditions satisfied by the non-linear terms are implemented. Alternatively, regardless of boundary conditions, the non-linear terms can be represented as Chebyshev polynomial series, which will also decrease the number of multiplications. In this section we address two questions: does (i) orthogonalization of the basis or (ii) another polynomial basis change the final result? The answer is “no”, but it requires some additional explanations. After choosing the truncation number we seek the solution in the form of (3.3.8). The solution belongs to the linear space ℒ = 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠 �𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑥𝑥 , 𝑧𝑧 ) , 𝝋𝝋 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑦𝑦 , 𝑧𝑧 ) � as in (3.3.8). This space consists of divergence-free vectors that satisfy all the LHBC of the problem, and their components are polynomials of maximum order 𝑚𝑚𝑠𝑠𝑥𝑥�𝐿𝐿 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝐿𝐿 ( 𝑥𝑥 , 𝑦𝑦 ) , 𝐿𝐿 ( 𝑦𝑦 , 𝑧𝑧 ) � + 4 in the x- direction, with similar expressions in the two other directions. We denote the order of space ℒ as 𝑁𝑁 ℒ and store all the basis functions of (3.3.8) in a set of vectors 𝑸𝑸 = { 𝒒𝒒 𝒊𝒊 } 𝑖𝑖=1𝑁𝑁 ℒ . Clearly, this set forms a basis in ℒ , ℒ = 𝑠𝑠𝑝𝑝𝑠𝑠𝑠𝑠 { 𝑸𝑸 } . Assume now another basis in ℒ , denoted as 𝑸𝑸� = { 𝒒𝒒� 𝑖𝑖 } 𝑖𝑖=1𝑁𝑁 ℒ . The connection between the two bases is given by a matrix ℬ of order 𝑁𝑁 ℒ as 𝑸𝑸� = ℬ𝑸𝑸 , 𝑸𝑸 = ℬ −1 𝑸𝑸� . (6.4.1) Elements of the matrix ℬ are solutions of the following system of linear algebraic equations (summation over repeating indices is assumed) 〈𝒒𝒒� 𝑖𝑖 , 𝒒𝒒� 𝑖𝑖 〉 = ℬ 𝑖𝑖𝑘𝑘 〈𝒒𝒒 𝑘𝑘 , 𝒒𝒒� 𝑖𝑖 〉 (6.4.2) Equation (6.4.1) can be interpreted as a transformation to another polynomial basis, a particular case of which is an orthonormal polynomial basis. In this case the matrix ℬ is then the operator of the Gram-Schmidt or Householder orthogonalization procedure. Assume now that projection of the solution 𝒗𝒗 on each of the bases is described by coefficients 𝑋𝑋 𝑖𝑖 and 𝑋𝑋� 𝑖𝑖 . Since these coefficients describe the orthogonal projection of the vector 𝒗𝒗 onto the same space, the result of projection onto either basis must be identical, i.e., 𝑑𝑑 ≈ 𝑋𝑋 𝑖𝑖 𝒒𝒒 𝑖𝑖 = 𝑋𝑋� 𝑖𝑖 𝒒𝒒� 𝑖𝑖 = 𝑋𝑋� 𝑖𝑖 ℬ 𝑖𝑖𝑘𝑘 𝒒𝒒 𝑘𝑘 , (6.4.3) from which it follows that 𝑋𝑋 𝑖𝑖 = 𝑋𝑋� 𝑖𝑖 ℬ 𝑖𝑖𝑘𝑘 , 𝑿𝑿 = ℬ 𝑇𝑇 𝑿𝑿� , 𝑿𝑿� = ( ℬ 𝑇𝑇 ) −1 𝑿𝑿 = ( ℬ −1 ) 𝑇𝑇 𝑿𝑿 (6.4.4) The Galerkin procedure applied with either of the two bases will result in two different ODE systems similar to (6.1.1): 𝐺𝐺 𝑖𝑖𝑖𝑖 𝑋𝑋̇ 𝑖𝑖 = 𝐿𝐿 𝑖𝑖𝑖𝑖 𝑋𝑋 𝑖𝑖 + 𝑁𝑁 𝑖𝑖𝑖𝑖𝑘𝑘 𝑋𝑋 𝑖𝑖 𝑋𝑋 𝑘𝑘 + 𝐹𝐹 𝑖𝑖 , 𝐺𝐺� 𝑖𝑖𝑖𝑖
𝑋𝑋�̇ 𝑖𝑖 = 𝐿𝐿� 𝑖𝑖𝑖𝑖
𝑋𝑋� 𝑖𝑖 + 𝑁𝑁� 𝑖𝑖𝑖𝑖𝑘𝑘
𝑋𝑋� 𝑖𝑖 𝑋𝑋� 𝑘𝑘 + 𝐹𝐹� 𝑖𝑖 , (6.4.5) where 𝐺𝐺 𝑖𝑖𝑖𝑖 = 〈𝒒𝒒 𝑖𝑖 , 𝒒𝒒 𝑖𝑖 〉 , 𝐿𝐿 𝑖𝑖𝑖𝑖 = 〈∆𝒒𝒒 𝑖𝑖 , 𝒒𝒒 𝑖𝑖 〉 , 𝑁𝑁 𝑖𝑖𝑖𝑖𝑘𝑘 = 〈�𝒒𝒒 𝑖𝑖 ∙ ∇�𝒒𝒒 𝑘𝑘 , 𝒒𝒒 𝑖𝑖 〉 , 𝐹𝐹 𝑖𝑖 = 〈𝒇𝒇 , 𝒒𝒒 𝑖𝑖 〉 (6.4.6) 𝐺𝐺� 𝑖𝑖𝑖𝑖 = 〈𝒒𝒒� 𝑖𝑖 , 𝒒𝒒� 𝑖𝑖 〉 , 𝐿𝐿� 𝑖𝑖𝑖𝑖 = 〈∆𝒒𝒒� 𝑖𝑖 , 𝒒𝒒� 𝑖𝑖 〉 , 𝑁𝑁� 𝑖𝑖𝑖𝑖𝑘𝑘 = 〈�𝒒𝒒� 𝑖𝑖 ∙ ∇�𝒒𝒒� 𝑘𝑘 , 𝒒𝒒� 𝑖𝑖 〉 , 𝐹𝐹� 𝑖𝑖 = 〈𝒇𝒇 , 𝒒𝒒� 𝑖𝑖 〉 (6.4.7) Note that the two ODE systems in (6.4.5) describe the orthogonal projection of the residual on the same linear space ℒ , so that the result again must be identical. However, it is not clear yet whether the coefficients 𝑋𝑋 𝑖𝑖 and 𝑋𝑋� 𝑖𝑖 yielded by solution of the two systems will be connected via Eq. (6.4.4). Let us evaluate how the matrices in (6.4.6) and (6.4.7) are connected. 𝐹𝐹� 𝑖𝑖 = 〈𝒇𝒇 , 𝒒𝒒� 𝑖𝑖 〉 = 〈𝒇𝒇 , ℬ 𝑖𝑖𝑝𝑝 𝒒𝒒 𝑝𝑝 〉 = ℬ 𝑖𝑖𝑝𝑝 𝐹𝐹 𝑝𝑝 , (6.4.8) 𝐺𝐺� 𝑖𝑖𝑖𝑖 = 〈𝒒𝒒� 𝑖𝑖 , 𝒒𝒒� 𝑖𝑖 〉 = 〈ℬ 𝑖𝑖𝑘𝑘 𝒒𝒒 𝑘𝑘 , ℬ 𝑖𝑖𝑝𝑝 𝒒𝒒 𝑝𝑝 〉 = ℬ 𝑖𝑖𝑘𝑘 ℬ 𝑖𝑖𝑝𝑝 𝐺𝐺 𝑘𝑘𝑝𝑝 , (6.4.9) 𝐿𝐿� 𝑖𝑖𝑖𝑖 = 〈∆𝒒𝒒� 𝑖𝑖 , 𝒒𝒒� 𝑖𝑖 〉 = 〈ℬ 𝑖𝑖𝑘𝑘 ∆𝒒𝒒 𝑘𝑘 , ℬ 𝑖𝑖𝑝𝑝 𝒒𝒒 𝑝𝑝 〉 = ℬ 𝑖𝑖𝑘𝑘 ℬ 𝑖𝑖𝑝𝑝 𝐿𝐿 𝑘𝑘𝑝𝑝 , (6.4.10) 𝑁𝑁� 𝑖𝑖𝑖𝑖𝑘𝑘 = 〈�𝒒𝒒� 𝑖𝑖 ∙ ∇�𝒒𝒒� 𝑘𝑘 , 𝒒𝒒� 𝑖𝑖 〉 = 〈�ℬ 𝑖𝑖𝑚𝑚 𝒒𝒒 𝑚𝑚 ∙ ∇�ℬ 𝑘𝑘𝑚𝑚 𝒒𝒒 𝑚𝑚 , ℬ 𝑖𝑖𝑝𝑝 𝒒𝒒 𝑝𝑝 〉 = ℬ 𝑖𝑖𝑚𝑚 ℬ 𝑘𝑘𝑚𝑚 ℬ 𝑖𝑖𝑝𝑝 𝑁𝑁 𝑝𝑝𝑚𝑚𝑚𝑚 . (6.4.11) Consider now how all the terms of the second system of (6.4.5) are expressed via matrices and unknowns of the first system 𝐺𝐺� 𝑖𝑖𝑖𝑖
𝑋𝑋�̇ 𝑖𝑖 = ℬ 𝑖𝑖𝑘𝑘 ℬ 𝑖𝑖𝑝𝑝 𝐺𝐺 𝑘𝑘𝑝𝑝 ℬ 𝑗𝑗𝑞𝑞− 𝑋𝑋 ̇ 𝑞𝑞 = ℬ 𝑖𝑖𝑝𝑝 𝐺𝐺 𝑘𝑘𝑝𝑝 𝑋𝑋 ̇ 𝑘𝑘 , (6.4.12) 𝐿𝐿� 𝑖𝑖𝑖𝑖
𝑋𝑋� 𝑖𝑖 = ℬ 𝑖𝑖𝑘𝑘 ℬ 𝑖𝑖𝑝𝑝 𝐿𝐿 𝑘𝑘𝑝𝑝 ℬ 𝑗𝑗𝑞𝑞− 𝑋𝑋 𝑞𝑞 = ℬ 𝑖𝑖𝑝𝑝 𝐿𝐿 𝑘𝑘𝑝𝑝 𝑋𝑋 𝑘𝑘 , (6.4.13) 𝑁𝑁� 𝑖𝑖𝑖𝑖𝑘𝑘
𝑋𝑋� 𝑖𝑖 𝑋𝑋� 𝑘𝑘 = ℬ 𝑗𝑗𝑚𝑚 ℬ 𝑘𝑘𝑙𝑙 ℬ 𝑑𝑑𝑝𝑝 𝑁𝑁 𝑝𝑝𝑚𝑚𝑙𝑙 ℬ 𝑗𝑗𝑞𝑞− 𝑋𝑋 𝑞𝑞 ℬ 𝑘𝑘𝑟𝑟− 𝑋𝑋 𝑟𝑟 = ℬ 𝑑𝑑𝑝𝑝 𝑁𝑁 𝑝𝑝𝑗𝑗𝑘𝑘 𝑋𝑋 𝑗𝑗 𝑋𝑋 𝑘𝑘 . (6.4.14) Substituting (6.4.8), (6.4.12)-(6.4.14) into the second system of (6.4.5) we obtain ℬ 𝑖𝑖𝑝𝑝 𝐺𝐺 𝑘𝑘𝑝𝑝 𝑋𝑋 ̇ 𝑘𝑘 = ℬ 𝑖𝑖𝑝𝑝 𝐿𝐿 𝑘𝑘𝑝𝑝 𝑋𝑋 𝑘𝑘 + ℬ 𝑑𝑑𝑝𝑝 𝑁𝑁 𝑝𝑝𝑗𝑗𝑘𝑘 𝑋𝑋 𝑗𝑗 𝑋𝑋 𝑘𝑘 + ℬ 𝑖𝑖𝑝𝑝 𝐹𝐹 𝑝𝑝 , (6.4.15) and multiplying (6.4.15) by ℬ −1 we return to the first system of (6.4.5). This proves that both systems of (6.4.5) yield identical solutions if their initial conditions are connected via eq. (6.4.4). To conclude, seeking other polynomial basis functions is useless, since we’ll arrive to exactly the same approximate solution. On the other hand, the orthonormalization procedure can be meaningful, since it does not change the solution, but avoids inverting the Gram matrix.
7. Inner products with arbitrary weight
For a scalar problem, e.g. Orr-Sommerfeld or Burgers equations, the choice of the weight function in the inner product (2.18) is arbitrary. In the case of unit or Chebyshev weight, the inner products can be evaluated analytically using properties of the Chebyshev polynomials listed in Appendix A. In other cases the Gauss quadrature formulae can be efficiently applied. An appropriate choice of the weight function can drastically improve the convergence, as was demonstrated in Gelfgat (2006) for the Burgers equation. For the incompressible Navier-Stokes equation, use of the unit weight function allows one to eliminate the pressure by the Galerkin projection. The unity weight also yields important conservation properties of the resulting ODE system (6.1.3) and (6.1.4). At the same time, if the weight function can be optimized such that the convergence is noticeably improved, then the total number of degrees of freedom in the resulting dynamical system can be decreased. In this case it can be reasonable to give up the nice properties of the unit weight and proceed with the optimized one. Then it will be necessary to solve the pressure equation (2.11), so that the ODE system with the algebraic constraints (2.16), (2.17) will be considered. In the following we follow Gelfgat (2006) to show how the algebraic constraints can be removed in the framework of the Galerkin approach. First of all, the Laplacian operator in pressure equation (2.11) must be supplied by boundary conditions. It was shown in Gelfgat (2006) that the boundary conditions proposed by Gresho & Sani (1987), which are limits of the momentum equation at the boundaries, yield a correct pressure field. The boundary conditions on the boundary 𝛤𝛤 are � 𝜕𝜕𝑝𝑝𝜕𝜕𝑛𝑛 � 𝛤𝛤 = 𝒏𝒏 ∙ � ∆𝒗𝒗 − ( 𝒗𝒗 ∙ ∇ ) 𝒗𝒗 − 𝜕𝜕𝒗𝒗𝜕𝜕𝜕𝜕 + 𝒇𝒇� 𝛤𝛤 , (7.1) Where 𝒏𝒏 is the normal to 𝛤𝛤 , and the boundary conditions for velocity are assumed to be steady. Note that Rempfer (2006) argued that the numerical solution of the pressure problem (2.11), (7.1), together with the momentum equation (2.7), do not yield a divergence-free solution for velocity. The global Galerkin method with divergence-free velocity basis functions described here does not have this problem, because the continuity equation is satisfied analytically by the basis functions, before the numerical process starts. Thus, any approximation of a solution is analytically divergence free. For the following we represent the pressure as a truncated Chebyshev series 𝑝𝑝 ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 , 𝑡𝑡 ) = ∑ 𝜗𝜗 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) 𝑇𝑇 𝑖𝑖 ( 𝑥𝑥 ) 𝑇𝑇 𝑖𝑖 ( 𝑦𝑦 ) 𝑇𝑇 𝑘𝑘 ( 𝑧𝑧 ) 𝑖𝑖 , 𝑖𝑖 , 𝑘𝑘 . (7.2) Here we cannot introduce any boundary conditions into the basis functions, because we cannot propose any sufficiently general change of variables that will make the boundary conditions (7.1) homogeneous for the pressure, so that they will be incorporated in the pressure basis functions. To obtain a problem for the unknown coefficients 𝜗𝜗 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) we perform Galerkin projections of the residuals of (2.11) and (7.1) on the Chebyshev basis (7.2), and require that the projections vanish. In other words, we apply the Galerkin method separately in the domain and on the boundaries. Clearly the total number of unknowns 𝜗𝜗 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) must be equal to the total number of equations used. Recalling that the velocity coefficients are stored in the vector 𝑿𝑿 , we store the additional coefficients of pressure 𝜗𝜗 𝑖𝑖𝑖𝑖𝑘𝑘 ( 𝑡𝑡 ) in the vector 𝒀𝒀 . After the Galerkin process is completed, the system of linear algebraic equations for 𝒀𝒀 has the following form 𝑄𝑄 𝐼𝐼𝐽𝐽 ( 𝑝𝑝 ) 𝑌𝑌 𝐽𝐽 = 𝐵𝐵 𝐼𝐼𝐽𝐽 ( 𝑝𝑝 ) 𝑋𝑋̇ 𝐽𝐽 + 𝐿𝐿 𝐼𝐼𝐽𝐽 ( 𝑝𝑝 ) 𝑋𝑋 𝐽𝐽 + 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 ( 𝑝𝑝 ) 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹 𝐼𝐼 ( 𝑝𝑝 ) (7.3) These equations are formed from the projections of Eqs. (2.11) and (7.1). The superscript ( 𝑝𝑝 ) is used to underline that all the matrices belong to the pressure problem. The matrix 𝑄𝑄 𝑖𝑖𝑖𝑖 ( 𝑝𝑝 ) is singular because the pressure is defined to within an additive constant. This singularity can be easily removed by, e.g., an additional requirement 𝜗𝜗 = 0 , after which the matrix 𝑄𝑄 ( 𝑝𝑝 ) is regular and its inverse is denoted as 𝑄𝑄 −1 . The Galerkin coefficients of velocity must be calculated using the equations (2.16), which take the following form 𝐺𝐺 𝐼𝐼𝐽𝐽
𝑋𝑋̇ 𝐽𝐽 = 𝑃𝑃 𝐼𝐼𝐽𝐽 𝑌𝑌 𝐽𝐽 + 𝐿𝐿 𝐼𝐼𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁 𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹 𝐼𝐼 . (7.4) Here the first term of the r.h.s. is the projection of the pressure gradient onto the velocity basis. Substituting 𝑌𝑌 𝑖𝑖 from Eqs. (7.3) into Eqs. (7.4) we obtain 𝐺𝐺�
𝐼𝐼𝐽𝐽
𝑋𝑋̇ 𝐽𝐽 = 𝐿𝐿�
𝐼𝐼𝐽𝐽 𝑋𝑋 𝐽𝐽 + 𝑁𝑁�
𝐼𝐼𝐽𝐽𝐾𝐾 𝑋𝑋 𝐽𝐽 𝑋𝑋 𝐾𝐾 + 𝐹𝐹� 𝐼𝐼 , (7.5) where 𝐺𝐺� = 𝐺𝐺 − 𝑃𝑃𝑄𝑄 −1 𝐵𝐵 ( 𝑝𝑝 ) , 𝐿𝐿� = 𝐿𝐿 + 𝑃𝑃𝑄𝑄 −1 𝐿𝐿 ( 𝑝𝑝 ) , 𝑁𝑁� = 𝑁𝑁 + 𝑃𝑃𝑄𝑄 −1 𝑁𝑁 ( 𝑝𝑝 ) , 𝐹𝐹� = 𝐹𝐹 + 𝑃𝑃𝑄𝑄 −1 𝐹𝐹 ( 𝑝𝑝 ) (7.6) Thus, after some analytical and numerical evaluations we obtain the ODE system (7.5), whose structure is equivalent to that of (6.1.1). Generally, the matrices of (7.5) do not obey the properties (6.1.3) and (6.1.4), however everything else that we have written about the system (6.1.1) is applicable also to (7.5). Some numerical experiments comparing convergence of the Galerkin method with Chebyshev and unit weights are reported in Gelfgat (2006). There the lid-driven cavity flow and convection in a laterally cavity were taken as test problems. It was found that the Chebyshev weight allows for a better resolution of boundary layers in the convection flow, but slows down the convergence for the lid-driven flow. This was attributed to the problem of approximating the corner discontinuities. Fig. 3. Comparison of numerical results with the experimental phograph of Escudier (1984). From Gelfgat et al. (1996). The flow at
H/R =3.25, Re=2752. (a) calculation with the Galerkin method using 34×34 basis functions. (b), (c) calculation with the finite volume method using 100×100 and 200×200 grids, respectively.
The optimization of the weight function has never been considered for a realistic fluid dynamics problem. The only optimization example is given in Gelfgat (2006) for the Burgers equation. Considering weight functions of the form of ( 𝑥𝑥 − 𝑥𝑥 ) −𝛼𝛼 , 𝛼𝛼 ≥ , it was found that the convergence is fastest at 𝛼𝛼 = 1.3 .
8. Solved problems and other applications of the method
The effectiveness of the Galerkin approach follows from a decrease in the number of degrees of freedom of a numerical model. The decrease is usually of an order of magnitude or more. One of the earliest examples of this is shown in Fig. 3. There we consider flow in a cylinder with rotating lid. Steady states of this flow exhibit the vortex breakdown phenomenon, which was experimentally studied by Escudier (1984). At certain Reynolds number a weak reverse circulation attached to the cylinder axis appears. In tall cylinders up to three distinct recirculation zones are observed. The intensity of the reverse vortices is 3-5 orders of magnitude less than that of the main meridional circulation, which makes its numerical modeling quite challenging. We show in Gelfgat et al. 1996 that the size and position of the reverse circulations is well reproduced with 34 with 34×34 basis functions (Fig. 3a), as well as with 200×200 finite volume grid (Fig. 3c), but is resolved inaccurately with the 100×100 grid (Fig. 3b). The total number of degrees of freedom of the Galerkin method is the number of unknown Galerkin coefficients of the meridional velocity vector and that of the azimuthal velocity component, which separates as a scalar in the axisymmetric formulation. In the finite volume approach this is the number of unknown functions multiplied by the number of grid nodes. Clearly the total number of degrees of freedom consumed by the Galerkin method, ∙ =2312 , is smaller than that of the finite volume method, ∙ = 120,000 , by about 1.5 orders of magnitude. The next example presented in Fig. 4 is a thermocapillary convective flow in a laterally heated cavity with the aspect ratio length/height = 4. The velocity boundary condition at the upper surface is 𝑑𝑑 = 0, 𝜕𝜕𝜕𝜕𝜕𝜕𝑦𝑦 = −𝑀𝑀𝑠𝑠 𝜕𝜕𝜕𝜕𝜕𝜕𝑥𝑥 , (8.1) where 𝜃𝜃 is the dimensionless temperature, 𝑀𝑀𝑠𝑠 = 𝑀𝑀𝑠𝑠 𝑃𝑃𝑟𝑟⁄ , and
𝑀𝑀𝑠𝑠 and
𝑃𝑃𝑟𝑟 are the Marangoni and Prandtl numbers, respectively (other details are in Gelfgat, 2007a). The problem was treated with the truncation of 50×20 basis functions with 100 collocation points at the upper surface to satisfy the boundary condition (8.1) The critical Marangoni number, corresponding to transition from steady to oscillatory flow regime, was calculated to be
𝑀𝑀𝑠𝑠 𝑐𝑐𝑟𝑟 = 4781 and the dimensionless critical frequency (imaginary part of the leading eigenvalue) 𝑖𝑖 𝑐𝑐𝑟𝑟 = 6674 . This result was never published because the value of 𝑖𝑖 𝑐𝑐𝑟𝑟 seemed to be too high compared with the values already known for buoyancy convection (Gelfgat et al., 1999a). Owing to the computer restrictions of that time, the convergence could not be rigorously checked. Later, the same problem was solved using a 800×200 stretched finite volume grid (Gelfgat, 2007a) and the result was 𝑀𝑀𝑠𝑠 𝑐𝑐𝑟𝑟 = 4779 with the same value of the critical frequency. The convergence of the critical values with resolution obtained by the finite volume method was found to be very slow. The reason for that is seen in Fig. 4. The streamlines and the isotherms are smooth, so that the velocity and temperature fields can be easily calculated. At the same time, the perturbation patterns exhibit
Fig. 4. Patterns of flow and amplitude of the most unstable perturbation at the critical Marangoni number for thermocapillary convection of low-Prandtl-number fluid ( Pr =0.015) in a cavity of aspect ratio length/height=4. From Gelfgat (2007a). very steep maxima, which must be numerically resolved to obtain correct critical values. We observe here again that the Galerkin method yielded the correct result with a much smaller number of degrees of freedom. The smaller number of degrees of freedom, as well as analytic representation of the Jacobian matrix (6.1.6), allows one to effectively apply Newton iteration for calculation of the steady states. To follow different solution branches one also can apply arc-length or similar continuation technique. This is illustrated in Fig. 5 for the convection in a cavity with partially heated sidewall (Erenburg et al., 2003), the boundary conditions for which were discussed in section 4. Since the boundary conditions are symmetric, the flow at low Grashof number 𝐺𝐺𝑟𝑟 is symmetric. The symmetry is broken at
𝐺𝐺𝑟𝑟 = 180 . The diagram in Fig. 5 shows difference between the Nusselt numbers calculated at the left and right vertical boundaries, so that in the symmetric state the difference is zero. After the symmetry breaks, we observe several interconnected solution branches with qualitatively different flow patterns. Those depicted in red are stable, and those depicted in color are oscillatory unstable. Although most of the steady state branches are unstable, we speculate that there exist multiple oscillatory states with similar flow patterns.
Fig. 5. Bifurcation diagram for convection in a vertical cavity with partially heated sidewall. Red lines correspond to stable steady states, blue lines to unstable ones. Pr =10. From Erenburg et al. (2003). The most well-known results obtained with this method are stability diagrams of swirling rotating disk – cylinder flow and buoyancy convection flows in laterally heated rectangular cavities. The first results on the three-dimensional instability of rotating disk – cylinder flow were obtained in Gelfgat et al. (2001a) using the Galerkin method with 30×30 basis functions. Since then several research groups validated these results experimentally and numerically. These comparisons are very convincing and are shown in Figs. 6 and 7. According to our results, the instability is axisymmetric for the cylinder aspect ratio (height/radius) varying between 1.6 and 2.7, and is three-dimensional outside of this interval. The three-dimensionality sets in with the azimuthal wavenumber 𝑘𝑘 = 2 at small aspect ratios, and with 𝑘𝑘 = 3 or in taller cylinders. Several later studies tried to reproduce these results either by straightforward integration in time, or by means of stability analysis, and fully confirmed our conclusions. The quantitative comparison was carried out for the critical Reynolds numbers and critical frequencies, as well as for the azimuthal mode number. It was possible also to confirm values of the aspect ratio at which the modes replace each other. Fig. 6. Stability diagram for flow in a cylinder covered by rotating disk. The lines correspond to results obtained by the Galerkin method (Gelfgat et al., 2001) for different azimuthal wavenumbers k in (3.4.1). Symbols show results independent studies . A similar comparison, but experimental, was carried out by Sørensen et al (2006, 2009). Their result is shown in Fig. 7. The oscillations were measured by LDA, and the flow azimuthal periodicity by PIV. The symbols in Fig. 7 show experimentally measured points, and their color corresponds to the azimuthal wavenumber as shown in the figure. The pioneer results of Escudier (1984) are also shown. Taking into account experimental uncertainties, the agreement between experiment of Sørensen et al (2006) and the numerical predictions of the Galerkin method is quite impressive. The agreement with earlier results of Escudier (1984) is observed only for the axisymmetric mode, but his experiments were based only on visualizations by the aluminum powder, and the instability was observed only as oscillations of the vortex breakdown bubbles. Other stability
Fig. 7. Stability diagram for flow in a cylinder covered by rotating disk. The lines correspond to results obtained by the Galerkin method (Gelfgat et al., 2001) for different azimuthal wavenumbers k in (3.4.1). Symbols show experimental results. Fig. 8. Stability diagram for convective flow in laterally heated cavities. The curves color corresponds to the number of convective rolls in the flow pattern. The flows are stable below or between the curves of the same color. Dashed regions correspond to the stability regions of similar flows with broken rotational symmetry. Pr =0. From Gelfgat et al. (1999) A G r c r results obtained by the described Galerkin method for similar swirling flows can be found in Gelfgat at al. (1996b) for flow in a cylinder with independently rotating top and bottom, and in Marques et al. (2003) for independently rotating top and sidewall. The neutral curves shown in Fig. 6 are plotted through up to a hundred calculated critical points. The next example, relating to the oscillatory instability of buoyancy convection flows in laterally heated cavities and shown in Fig. 8, required several hundreds of critical points to complete the study. The calculations were performed with up to 60×20 basis functions. With increasing cavity aspect ratio A = length/height, and at large enough Grashof number, the single convective cell splits into several cells. The number of cells grows with the aspect ratio. At the same time several steady states with a different number of rolls can be stable at the same set of governing parameters, as shown in Fig. 9. The transition from one number of cells to another takes place at points where a neutral curve of a certain color continues with a different color. At these points the flows with, e.g., two and three rolls are indistinguishable. Other results on stability of convection in rectangular cavities can be found in Gelfgat et al. (1996, 1999a,b), Erenburg et al. (2003), and Gelfgat (2004). One particularly interesting result reported in Gelfgat (2004) showed that weakly non-linear approximation of limit cycles (6.9)-(6.11) yields results that are very close to those obtained by independent straightforward time integration. Several studies have been devoted to three-dimensional instabilities of axisymmetric buoyancy convection in vertical cylindrical containers. These studies were started in Gelfgat et al. (1999c), where we were able to reproduce a nice experimental result showing the breaking of axisymmetry leading to a spoke pattern with azimuthal wavenumber 𝑘𝑘 = 9 . Later we studied cylinders with non-uniformly heated sidewalls that mimicked conditions of Bridgman crystal growth (Gelfgat et al., 2000, 2001b). Later works were devoted to axisymmetric flows driven by rotating or traveling magnetic field (Gelfgat & Gelfgat, 2004; Gelfgat, 2005). Most of these results were reviewed in more detail in Gelfgat & Bar-Yoseph (2004). branch 1branch 2branch 3 Fig. 9. Three distinct stable steady states found at Pr =0, A=7, Gr=88,000. From Gelfgat & Bar-Yoseph (2004). Additional opportunities are provided by analytical representation of numerical solution via the Galerkin series. Clearly, one can differentiate or integrate the series without any noticeable loss of accuracy, contrary to low-order methods. An obvious example is the calculation of flow trajectories using a previously calculated steady or time-dependent flow. Since the velocity field is defined analytically in the whole domain, wherever the liquid particle arrives, its velocity is known without the need to interpolate between grid nodes. This fact was used in Gelfgat (2002), where trajectories were calculated over very long time to obtain a Poincare map in the midplane. Another application of the divergence-free bases (3.3.2), (3.3.3), and (3.3.6) is visualization of three-dimensional incompressible flows, as described in Gelfgat (2014). Without going into detail, we only mention that projection of flow on each divergence-free set can be interpreted as a divergence-free projection on coordinate planes 𝑥𝑥 = 𝑐𝑐𝑜𝑜𝑠𝑠𝑠𝑠𝑡𝑡 , 𝑦𝑦 = 𝑐𝑐𝑜𝑜𝑠𝑠𝑠𝑠𝑡𝑡 , 𝑧𝑧 = 𝑐𝑐𝑜𝑜𝑠𝑠𝑠𝑠𝑡𝑡 , which results in two-component divergence-free fields. These can be described by an analog of the streamfunction. Assembling all the planes, e.g., 𝑥𝑥 = 𝑐𝑐𝑜𝑜𝑠𝑠𝑠𝑠𝑡𝑡 , we obtain a scalar 3D function whose isosurfaces are tangent to the projected vectors. Three such projections of three sets of coordinate planes complete the visualization of a three-dimensional flow field. Details and illustrations can be found in Gelfgat (2014).
9. Similar approaches in studies of other authors
As mentioned, the idea to use linear superpositions of Chebyshev polynomials for definition of basis functions satisfying linear and homogeneous boundary conditions was introduced by Orszag (1971a,b) for the homogeneous two-point Dirichlet problem. Since then, similar linear-superposition-based basis functions were used to solve one-dimensional problems for, e.g., Orr-Sommerfeld and boundary layer equations, by Holte (1983), Pasquarelli (1991), Yueh & Weng (1996), Yang (1997), Borget et al. (2001), Yahata (2001), Bistrian et al. (2009), and Buffat & Le Penven (2011). In these works the basis functions were constructed from the Chebyshev polynomials. Recently, Wan & Yu (2017) applied the same idea to Legendre polynomials. Grants & Gerbeth (2001), Uhlmann & Nagata (2006), and Batina et al. (2009) used the same approach for a two-dimensional flow field, but their basis functions were not divergence-free. Picardo et al. used linear superpositions for a two-fluid problem, as was done in Gelfgat et al. (2001d). Moser et al. (1983) proposed to multiply linear superpositions of the Chebyshev polynomials by powers of the Chebyshev weight function in order to make better use of the polynomial orthogonality properties. For problems with two periodic directions, these authors constructed a divergence free basis, in which the non-periodic direction was treated by linear superpositions of the Chebyshev polynomials multiplied by additional weight-dependent functions. Such functions were used for either coordinate or projection bases in the weighted residuals method by Ganske et al. (1994), Goddeferd & Lollini (1999), Kerr (1996). It should be noted that multiplication by either function makes evaluation of derivatives and computations of their Galerkin projections more complicated. Yahata (1999) solved a problem of buoyancy convection in laterally heated cavities similar to those treated by Gelfgat & Tanasawa (1994) and later by Gelfgat et al. (1997 1999a,b). He used linear superpositions of Chebyshev polynomials to construct basis functions for the temperature and the streamfunction with subsequent orthogonalization of the basis. The inner product was formulated with an arbitrary weight function, however it is not clear which weight was used. By evaluating derivatives of the streamfunction basis of Yahata (1999) one would obtain the two-dimensional basis (3.2.3), so that both formulations are equivalent in the 2D case. An extension of Yahata’s approach to 3D formulation would require replacement of the stream function by a vector potential, which would complicate the formulation. Suslov & Paolucci (2002) also solved similar convection problem in a cavity with coordinate functions (3.2.3). For the projection system they used the same functions multiplied by the Chebyshev weight, which made the Gram matrix sparser and allowed to use the fast Fourier transform (FFT) to evaluate non-linear terms of the dynamical system. The whole approach was used for straightforward integration in time, but did not exhibit much advantage compared to other methods and, to the best of the author’s knowledge, had no further continuation.
10. To conclude: what else can be done?
To discuss further possible implementations of this Galerkin approach we mention that with the present growth of computational power and state-of-the-art methods of numerical linear algebra, the solution of two-dimensional and quasi-two-dimensional problems became feasible, and sometimes more efficient, with lower order methods. A very popular methodology of turning a time-stepping code into a steady state / stability solver can be found in Boronska & Tuckerman (2010a,b) and Tuckerman et al (2018). Another possible methodology together with several examples are given in Gelfgat (2007a,b). For these problems the Galerkin approach can be more suitable for weakly non-linear bifurcations analysis. It is not clear, however, whether results applicable only for small supercriticality will justify the effort. One of possible applications of the method is consideration of fully three-dimensional flows, steady and unsteady, in axisymmetric domains. In these problems the bases (3.4.5) and (3.4.6) can be combined with the Fourier decomposition in the circumferential direction, so that the Gram matrices will separate for each Fourier mode and thus not be too large, so that they will be easily inverted. Finally, one obtains an ODE system, where equations corresponding to different Fourier modes will be coupled via the non-linear terms. This system allows for computation of steady states, path-following, stability analysis and time-dependent calculations (see, e.g., Boronska & Tuckerman, 2010a,b). Possibly, the most challenging task would be to develop a fully three-dimensional solver for flow in a three-dimensional rectangular box. This would require orthogonalization of the whole set of bases (3.3.2), (3.3.3), and (3.3.6) with subsequent effective treatment of the non-linear terms. If successful, this approach would allow one to have steady state, stability, and time-dependent solvers within a single computational model, with which very complicated flows can be studied. It should be mentioned here that Krylov subspace based solvers are sometimes thought to be applicable only to sparse matrices, for which matrix-vector products can be quickly evaluated. However, matrix-vector products are also fast for the method that we have described, even though all the related matrices are densely filled. As mentioned in the very beginning of this paper, the Galerkin approach that we have described is limited to simple domains, which must be curvilinear rectangles, in other words, regions bounded by coordinate surfaces. This is a very stringent restriction since it excludes a very big set of important problems, in which the boundaries have a more complicated shape. Another restriction for implementation of this method is flows with deformable interfaces. One of the ways to solve such problems on fixed grids is the immersed boundary method and/or the diffuse interface approach (not described here). Implementation of these approaches for the spectral method would require a good approximation of the delta function, which can be difficult to do using smooth polynomials. Finally, we emphasize a technique for flow visualization made in Gelfgat (2014). This is applicable to flows calculated by any of numerical method, and can be very helpful for understanding the topology of complicated three-dimensional flows. References:
Batina J, Blancher S, Amrouche C, Batchi M, Creff R (2009) Convective heat transfer augmentation through vortex shedding in sinusoidal constricted tube, Int. J. Numer. Meths Heat Fluid Flow, 19: 374-395. Blackburn HM, Lopez JM (2000) Symmetry breaking of the flow in a cylinder driven by a rotating endwall Phys. Fluids 12: 2698–2701. Bistrian DA, Dragomirescu IA, Muntean S, Topor M (2009) Numerical methods for convective hydrodynamic stability of swirling flows, Proc. 13th WSEAS Int. Conf. on Systems, 283-288. Borget V, Bdéoui F, Soufiani A, Le Quéré P (2001) The transverse instability in a differentially heated vertical cavity filled with molecular radiating gases. I. Linear stability analysis, Phys. Fluids, 13:1492-1507.
Borońska K Tuckerman LS (2010a) Extreme multiplicity in cylindrical Rayleigh -Bénard convection. I. Time dependence and oscillations. Phys. Rev. E, 81: 036320. B orońska K, Tuckerman LS (2010b) Extreme multiplicity in cylindrical Rayleigh -Bénard convection. II. Bifurcation diagram and symmetry classification. Phys. Rev. E, 81: 036321. Boyd JP (2000) Chebyshev and Fourier Spectral Methods. Dover Publications, New York. Buffat M, Le Penven L (2011) A spectral fictitious domain method with internal forcing for solving elliptic PDEs, J. Comput. Phys., 230: 2433-2450. Canuto C, Quarteroni A, Hussaini MY, Zhang TA (2006) Spectral Methods. Fundamentals in Single Domains. Springer, Berlin Heidelberg New York. Dean WR (1928) ‘Fluid flow in a curved channel, Proc. R. Soc. London, Ser.A 121: 402-420. Dumas G, Leonard A (1994) A divergence free spectral expansion method for three-dimensional flows in spherical gap geometries, J. Comput. Phys., 11: 205-219. Erenburg V, Gelfgat AY, Kit E, Bar-Yoseph PZ, Solan A (2003) Multiple states, stability, bifurcations of natural convection in rectangular cavity with partially heated vertical walls, J. Fluid Mech., 492: 63-89. Escudier MP (1984) Observation of the flow produced in a cylindrical container by a rotating endwall, Exp. Fluids, 2:189-196. Fletcher CAJ (1984) Computational Galerkin Methods. Springer, New York. Ganske A, Gebhardt T, Grossmann S (1994) Modulation effects along stability border in Taylor–Couette flow, Phys. Fluids, 6: 3823-3832. Gelfgat AY (1988) Instability and oscillatory supercritical regimes of free convection in a laterally heated square cavity, PhD (Cand. Sci.) thesis. Latvian State University, Riga, Latvia. Gelfgat AY, Tanasawa I (1994) Numerical analysis of oscillatory instability of buoyancy convection with the Galerkin spectral method, Numer. Heat Transfer. Part A: Applications, 25: 627-648. Gelfgat AY, Bar-Yoseph PZ, Solan A (1996a) Stability of confined swirling flow with, without vortex breakdown., J. Fluid Mech., 311: 1-36. Gelfgat AY, Bar-Yoseph PZ, Solan A (1996b) Steady states, oscillatory instability of swirling flow in a cylinder with rotating top, bottom, Phys. Fluids, 8: 2614-2625. Gelfgat AY, Bar-Yoseph PZ, Yarin AL (1997) On oscillatory instability of convective flows at low Prandtl number, J. Fluids Eng., 119: 823-830. Gelfgat AY, Bar-Yoseph PZ, Yarin AL (1999a) Stability of multiple steady states of convection in laterally heated cavities, J. Fluid Mech., 388: 315-334. Gelfgat AY, Bar-Yoseph PZ, Yarin AL (1999b) Non-Symmetric convective flows in laterally heated rectangular cavities, Int. J. Computational Fluid Dynamics, 11: 261-273. Gelfgat AY, Bar-Yoseph PZ, Solan A, Kowalewski T (1999c) An axisymmetry- breaking instability in axially symmetric natural convection, Int. J. Transport Phenomena, 1: 173-190. Gelfgat AY (1999) Different modes of Rayleigh-Bénard instability in two-, three-dimensional rectangular enclosures, J. Comput. Phys. 156: 300-324. Gelfgat AY, Bar-Yoseph PZ, Solan A (2000) Axisymmetry breaking instabilities of natural convection in a vertical Bridgman growth configurations, J. Cryst. Growth, 220: 316-325. Gelfgat AY (2001) Two-, three-dimensional instabilities of confined flows: numerical study by a global Galerkin method, Computational Fluid Dynamics Journal, 9: 437-448. Gelfgat AY, Bar-Yoseph PZ, Solan A. (2001a) Three-dimensional instability of axisymmetric flow in a rotating lid - cylinder enclosure, J. Fluid Mech., 438: 363-377. Gelfgat AY, Bar-Yoseph PZ, Solan A (2001b) Effect of axial magnetic field on three-dimensional instability of natural convection in a vertical Bridgman growth configuration, J Cryst.Growth, 230: 63-72. Gelfgat AY, Bar-Yoseph PZ (2001c) The effect of an external magnetic field on oscillatory instability of convective flows in a rectangular cavity, Phys. Fluids, 13: 2269-2278. Gelfgat AY, Yarin AL, Bar-Yoseph PZ (2001d) Three-dimensional instability of a two-layer Dean flow, Phys. Fluids, 13: 3185-3195. Gelfgat AY (2002) Three-dimensionality of trajectories of experimental tracers in a steady axisymmetric swirling flow: effect of density mismatch, Theor. Comput. Fluid Dyn., 16: 29-41. Gelfgat AY (2004) Stability, slightly supercritical oscillatory regimes of natural convection in a 8:1 cavity: solution of benchmark problem by a global Galerkin method, Int. J. Numer. Meths. Fluids, 44: 135-146. Gelfgat AY, Bar-Yoseph PZ (2004) Multiple solutions, stability of confined convective, swirling flows – a continuing challenge, Int. J. Numer. Meth. Heat, Fluid Flow, 14: 213-241. Gelfgat YM, Gelfgat AY (2004) Experimental, numerical study of rotating magnetic field driven flow in cylindrical enclosures with different aspect ratios, Magnetohydrodynamics, 40: 147-160. Gelfgat AY (2005) On three-dimensional instability of a traveling magnetic field driven flow in a cylindrical container, J. Crystal Growth, 279: 276-288. Gelfgat AY (2006) Implementation of arbitrary inner product in global Galerkin method for incompressible Navier-Stokes equation, J. Comput. Phys., 211: 513-530. Gelfgat AY (2007a) Stability of convective flows in cavities: solution of benchmark problems by a low-order finite volume method, Int. J. Numer. Meths. Fluids, 53: 485-506. Gelfgat AY (2007b) Three-dimensional instability of axisymmetric flows: solution of benchmark problems by a low-order finite volume method, Int. J. Numer. Meths. Fluids, 54: 269-294. Gelfgat AY (2014) Visualization of three-dimensional incompressible flows by quasi-two-dimensional divergence-free projections, Computers & Fluids, 97 : Hassard BD, Kazarinoff ND, Wan YH (1981 ) Theory and Applications of Hopf Bifurcation . London, Math. Soc. Lecture Note Series, vol.41. Holte S (1983) Numerical experiments with a three-dimensional model of an enclosed basin, Continental Shelf Research, 2: 301-315. Iwatsu R (2005) Numerical study of flows in a cylindrical container with rotating bottom and top free surface, J. Phys. Soc. Japan, 74: 333:344. Kerr R (1996) Rayleigh number scaling in numerical convection, J. Fluid Mech., 310: 139-179. Marques F, Lopez JM (2001). Precessing vortex breakdown mode in an enclosed cylinder flow. Phys. Fluids, 13: 1679–1682. Marques F, Lopez JM, Shen J (2002). Mode interactions in an enclosed swirling flow: A double Hopf bifurcation between azimuthal wavenumbers 0 and 2. J. Fluid Mech., 455: 263–281. Marques F., Gelfgat AY, Lopez JM (2003) Tangent double Hopf bifurcation in a differentially rotating cylinder flow, Phys. Rev. E, 68: 06310-1 – 06310-13 . Mason JC, Handscomb DC (2003) Chebyshev Polynomials, CRC Press, London. Meseguer A , Mellibovsky F (2007) On a solenoidal Fourier–Chebyshev spectral method for stability analysis of the Hagen–Poiseuille flow, Appl. Numer. Math., 920-938. Moser RD, Moin P, Leonard A (1983) A spectral numerical method for the Navier-Stokes equations with applications to Taylor-Couette flow, J. Comput. Phys., 52: 524-544. Nore C, Tuckerman LS, Daube O, Xin S (2003) The 1:2 mode interaction in exactly counter-rotating von Karman swirling flow, J. Fluid Mech., 477: 51–88. Orszag SA (1971a) Galerkin approximations to flows within slabs, spheres, and cylinders, Phys. Rev. Lett., 26: 1100-1103. Orszag SA (1971b) Accurate solution of the Orr-Sommerfeld stability equation, J. Fluid Mech, 50: 689-703. Pasquali F (1991) Domain decomposition for spectral approximation to Stokes equations via divergence-free functions, Appl. Umr. Math., 8: 493-51. Paszkowski S (1975) Numerical Applications of Chebyshev Polynomials and Series, PWN, Warsaw (in Polish). Picardo JR, Garg P, Pushpavanam S (2015) Centrifugal instability of stratified two-phase flow in a curved channel, Phys. Fluids, 27: 054106. Rempfer D (2006) On boundary conditions for incompressible Navier-Stokes problems, Appl. Mech. Rev., 59: 107-125. Rubinov A, Erenburg V, Gelfgat AY, Kit E, Bar-Yoseph PZ, Solan A (2004) Three-dimensional instabilities of natural convection in a vertical cylinder with partially heated sidewalls, J. Heat Transfer, 126: 586-599. Sørensen JN, Naumov I, Mikkelsen R (2006) Experimental investigation of three-dimensional flow instabilities in a rotating lid-driven cavity, Exp. Fluids, 41: 425-440. Sørensen JN, Gelfgat AY, Naumov I, Mikkelsen R (2009) Experimental and numerical results on three-dimensional instabilities in a rotating disk–tall cylinder flow, Phys. Fluids, 21: 054102. Suslov SA, Paolucci S (2002) A Petrov-Galerkin method for flows in cavities: enclosure of aspect ratio 8, Int. J. Numer. Meths. Fluids, 40: 999-1007. Tuckerman LS, Langham J, Willis A., (2018) Stokes preconditioning in Channelflow and Openpipeflow for steady states and traveling waves. In:
Shifted Chebyshev polynomials of the first and the second type shifted onto the interval ≤ 𝑥𝑥 ≤ are defined as 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ) = 𝑐𝑐𝑜𝑜𝑠𝑠 [ 𝑠𝑠 𝑠𝑠𝑟𝑟𝑐𝑐𝑐𝑐𝑜𝑜𝑠𝑠 (2 𝑥𝑥 − 𝑈𝑈 𝑛𝑛 ( 𝑥𝑥 ) = 𝑠𝑠𝑖𝑖𝑛𝑛 [( 𝑛𝑛+1 ) 𝑎𝑎𝑟𝑟𝑐𝑐𝑐𝑐𝑎𝑎𝑠𝑠 ( )] 𝑠𝑠𝑖𝑖𝑛𝑛 [ 𝑎𝑎𝑟𝑟𝑐𝑐𝑐𝑐𝑎𝑎𝑠𝑠 ( )] , 0 ≤ 𝑥𝑥 ≤ (A1) The two systems { 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 )} 𝑛𝑛=1∞ and { 𝑈𝑈 𝑛𝑛 ( 𝑥𝑥 )} 𝑛𝑛=1∞ form bases in 𝐿𝐿 [0,1] and are connected via the relation 𝑇𝑇 𝑛𝑛′ ( 𝑥𝑥 ) = 2( 𝑠𝑠 + 1) 𝑈𝑈 𝑛𝑛−1 ( 𝑥𝑥 ) , (A2) which resembles connection between sine and cosine. Values of the polynomials and their derivatives in the points 𝑥𝑥 = 0 and 𝑥𝑥 = 1 that are needed to define basis functions for different boundary conditions are 𝑇𝑇 𝑛𝑛 (0) = ( − 𝑛𝑛 , 𝑇𝑇 𝑛𝑛 (1) = 1 (A3) 𝑈𝑈 𝑛𝑛 (0) = ( − 𝑛𝑛 ( 𝑠𝑠 + 1), 𝑈𝑈 𝑛𝑛 (1) = 𝑠𝑠 + 1 (A4) 𝑇𝑇 𝑛𝑛′ (0) = ( − 𝑛𝑛 𝑠𝑠 , 𝑇𝑇 𝑛𝑛′ (1) = 2 𝑠𝑠 (A5) 𝑈𝑈 𝑛𝑛′ (0) = ( − 𝑛𝑛 𝑠𝑠 ( 𝑠𝑠 + 1)( 𝑠𝑠 + 2) 3 ⁄ , 𝑈𝑈 𝑛𝑛′ (1) = 𝑠𝑠 ( 𝑠𝑠 + 1)( 𝑠𝑠 + 2) 3 ⁄ (A6) For the following we assume that for 𝑘𝑘 > 0 , 𝑇𝑇 −𝑘𝑘 ( 𝑥𝑥 ) = 𝑇𝑇 𝑘𝑘 ( 𝑥𝑥 ) and 𝑈𝑈 −𝑘𝑘 ( 𝑥𝑥 ) = 𝑈𝑈 −𝑘𝑘−1 ( 𝑥𝑥 ) = 𝑈𝑈 𝑘𝑘 ( 𝑥𝑥 ) . To evaluate inner products we need to decompose the polynomial derivatives and the polynomials products into polynomial sums. For the multiplication of a polynomial by a polynomial we have 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ) 𝑇𝑇 𝑚𝑚 ( 𝑥𝑥 ) = �𝑇𝑇 𝑛𝑛+𝑚𝑚 ( 𝑥𝑥 )+ 𝑇𝑇 𝑛𝑛−𝑚𝑚 ( 𝑥𝑥 ) � (A7) 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ) 𝑈𝑈 𝑚𝑚 ( 𝑥𝑥 ) = �𝑈𝑈 𝑛𝑛+𝑚𝑚 ( 𝑥𝑥 )+ 𝑈𝑈 𝑛𝑛−𝑚𝑚 ( 𝑥𝑥 ) � (A8) 𝑈𝑈 𝑛𝑛 ( 𝑥𝑥 ) 𝑈𝑈 𝑚𝑚 ( 𝑥𝑥 ) = ∑ 𝑈𝑈 𝑚𝑚−𝑛𝑛+2𝑘𝑘 ( 𝑥𝑥 ) 𝑛𝑛𝑘𝑘=0 (A9) The derivatives can be represented as Chebyshev series as 𝑇𝑇 𝑛𝑛+𝑚𝑚𝑚𝑚+1 ( 𝑥𝑥 ) = 2 𝑚𝑚+2 𝑙𝑙 ! ( 𝑠𝑠 + 𝑙𝑙 ) ∑ 𝑠𝑠 𝑛𝑛−1−2𝑖𝑖 �𝑗𝑗 + 𝑙𝑙𝑙𝑙 � �𝑠𝑠 − 𝑗𝑗 + 𝑙𝑙 − 𝑙𝑙 � 𝑇𝑇 𝑛𝑛−1−2𝑖𝑖 ( 𝑥𝑥 ) [( 𝑛𝑛−1 ) ] 𝑖𝑖=0 (A10) 𝑈𝑈 𝑛𝑛+𝑚𝑚𝑚𝑚+1 ( 𝑥𝑥 ) = 2 𝑚𝑚+3 ( 𝑙𝑙 + 1)! ∑ 𝑠𝑠 𝑛𝑛−1−2𝑖𝑖 �𝑗𝑗 + 𝑙𝑙 + 1 𝑙𝑙 + 1 � �𝑠𝑠 − 𝑗𝑗 + 𝑙𝑙𝑙𝑙 + 1 � 𝑇𝑇 𝑛𝑛−1−2𝑖𝑖 ( 𝑥𝑥 ) [( 𝑛𝑛−1 ) ] 𝑖𝑖=0 (A11) 𝑠𝑠 = 12 , 𝑠𝑠 𝑚𝑚>0 = 1 For example, 𝑇𝑇 𝑛𝑛′ ( 𝑥𝑥 ) = 4 𝑠𝑠 ∑ 𝑠𝑠 𝑛𝑛−1−2𝑖𝑖 𝑇𝑇 𝑛𝑛−1−2𝑖𝑖 ( 𝑥𝑥 ) [( 𝑛𝑛−1 ) ] 𝑖𝑖=0 (A12) Here �𝑚𝑚𝑠𝑠 � is the binomial coefficient. After the basis functions are built, the Galerkin projections can be computed with the unity weight 〈𝑓𝑓 , 𝑔𝑔〉 = ∫ 𝑓𝑓 ( 𝑥𝑥 ) 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 , (A13) or with the Chebyshev weight 〈𝑓𝑓 , 𝑔𝑔〉 𝐶𝐶ℎ = ∫ ( 𝑥𝑥 − 𝑥𝑥 ) −1 2⁄ 𝑓𝑓 ( 𝑥𝑥 ) 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 , (A14) or with an arbitrary weight. For example, in Gelfgat (2005) we used 〈𝑓𝑓 , 𝑔𝑔〉 𝑎𝑎 = ∫ ( 𝑥𝑥 − 𝑥𝑥 ) −𝛼𝛼 𝑓𝑓 ( 𝑥𝑥 ) 𝑔𝑔 ( 𝑥𝑥 ) 𝑑𝑑𝑥𝑥 , 𝛼𝛼 < 1 (A15) All the inner products needed to complete the Galerkin procedure can be evaluated analytically if either the unity or Chebyshev weight is implied. The following relations (base products) are needed for that 〈𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ), 𝑇𝑇 𝑚𝑚 ( 𝑥𝑥 ) 〉 = [1 + ( − 𝑛𝑛+𝑚𝑚−1 ] � − + − � (A16) 〈𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ), 𝑇𝑇 𝑚𝑚 ( 𝑥𝑥 ) 〉 𝐶𝐶ℎ = 𝑠𝑠 𝑛𝑛 𝜋𝜋𝛿𝛿 𝑛𝑛𝑚𝑚 (A17) And 𝛿𝛿 𝑛𝑛𝑚𝑚 is the Kronecker symbol. The relation (A2), (A7)-(A16) allow one to reduce all the inner products to the above ones. In the case of arbitrary inner product the base products 〈𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ), 𝑇𝑇 𝑚𝑚 ( 𝑥𝑥 ) 〉 must be evaluated numerically, which is usually done using the Gauss quadrature. Then all the other products can be expressed as sums using relations (A2) and (A7)-(A11). Additionally, for evaluation of inner products in orthogonal curvilinear coordinates, one may need the following relation 𝑥𝑥 𝑚𝑚 = 2 −2𝑚𝑚+1 ∑ 𝑠𝑠 𝑚𝑚−𝑖𝑖 � 𝑚𝑚𝑑𝑑 � 𝑇𝑇 𝑚𝑚−𝑖𝑖 ( 𝑥𝑥 ) 𝑚𝑚𝑖𝑖=0 (A18) Finally, to calculate the shifted Chebyshev polynomials in a point, the following recurrent formulae can be used 𝑇𝑇 𝑛𝑛 ( 𝑥𝑥 ) = (4 𝑥𝑥 − 𝑇𝑇 𝑛𝑛−1 ( 𝑥𝑥 ) − 𝑇𝑇 𝑛𝑛−2 ( 𝑥𝑥 ) (A19) 𝑈𝑈 𝑛𝑛 ( 𝑥𝑥 ) = (4 𝑥𝑥 − 𝑈𝑈 𝑛𝑛−1 ( 𝑥𝑥 ) − 𝑈𝑈 𝑛𝑛−2 ( 𝑥𝑥 ) (A20) Further details can be found in the books of Paszkowski (1975) and Mason & Handscomb (2003).(A20) Further details can be found in the books of Paszkowski (1975) and Mason & Handscomb (2003).