A Note on the Solution of the Von Kármán Equations Using Series and Chebyshev Spectral Methods

The classical von Kármán equations governing the boundary layer flow induced by a rotating disk are solved using the spectral homotopy analysis method and a novel successive linearisation method. The methods combine nonperturbation techniques with the Chebyshev spectral collocation method, and this study seeks to show the accuracy and reliability of the two methods in finding solutions of nonlinear systems of equations. The rapid convergence of the methods is determined by comparing the current results with numerical results and previous results in the literature.

Most natural phenomena can be described by nonlinear equations that, in general, are not easy to solve in closed form. The search for computationally efficient, robust, and easy to use numerical and analytical techniques for solving nonlinear equations is therefore of great interest to researchers in engineering and science. The study of the steady, laminar, and axially symmetric viscous flow induced by an infinite disk rotating steadily with constant angular velocity was pioneered by von Kármán [1]. He showed that the Navier-Stokes equations could be reduced to a set of ordinary differential equations and solved using an approximate integral method. His solution, however, contained errors that were later corrected by Cochran [2] by patching together two series expansions.

Numerical and semianalytical methods including the cubic Hermite finite element, pseudospectral, Galerkin-B-Spline, and Chebyshev-collocation methods have been used previously to find solutions of the von Kármán equations [3–6]. These methods have their shortcomings, including instability, and hence the last few decades have seen the popularization of a number of new perturbation or nonperturbation techniques such as the Adomian decomposition method [7], the Lyapunov artificial small parameter method [8], the homotopy perturbation method [9, 10], and the homotopy analysis method [11].

The homotopy analysis method (HAM) was used recently by Yang and Liao [12] to find explicit, purely analytic solutions of the swirling von Kármán equations. Turkyilmazoglu [13] used the homotopy analysis method to solve the equations governing the flow of a steady, laminar, incompressible, viscous, and electrically conducting fluid due to a rotating disk subjected to a uniform suction and injection through the walls in the presence of a uniform transverse magnetic field. For this extended form of the von Kármám problem, the homotopy analysis method, however, produced secular terms in the series solution. Turkyilmazoglu [13] overcame this weakness by using initial guesses based on Ackroyd's (see the work of Ackroyd [14]) exponentially decaying functions, and a new linear operator which resulted in a method capable of tracking the shape of the exact solution. An alternative approach that serves to address these and other limitations of the HAM is the spectral homotopy analysis method; see the work of Motsa et al. [15, 16]. It is an efficient hybrid method that blends the HAM algorithm with Chebyshev spectral methods. The method retains the proven qualities of the HAM while effectively using Chebyshev polynomials as basis functions to ensure rapid convergence of the solution series. A novel quasilinearisation method—the successive linearisation method (see the work of Makukula et al. [17] and Motsa and Sibanda [18])—promises further improvement in accuracy and convergence rates compared to both the HAM and the SHAM.

In this study we apply the spectral homotopy analysis method (SHAM) and the successive linearisation method (SLM) to solve the von Kármán equations. The results are compared with those in the literature [11, 12] and against numerical approximations. Comparison of current results is further made with the recent results of Turkyilmazoglu [13] that include suction/injection and an applied magnetic field. We show, inter alia, that notwithstanding the fact that these two methods may involve more computations per step than the HAM, both the SHAM and SLM are efficient, robust, and converge much more rapidly compared to the standard homotopy analysis method.

Our focus in this section is on the original von Kármán equation for the steady, laminar, axially symmetric viscous flow induced by an infinite disk rotating steadily with angular velocity about the -axis with the fluid confined to the half-space above the disk. In cylindrical coordinates the equations of motion are given by

subject to the nonslip boundary conditions on the disk and boundary conditions at infinity

where is the fluid density, is the kinematic viscosity coefficient, is the pressure, , , and are the velocity components in the radial, azimuthal, and axial directions, respectively, and is the constant angular velocity. Using von Kármán's similarity transformations [1]

where is a nondimensional distance measured along the axis of rotation, the governing partial differential equations (2) reduce to a set of ordinary differential equations:

subject to the boundary conditions

Substituting (2.7) into (2.4) and (2.5) yields

subject to the boundary conditions

Equations (2.9) with the prescribed boundary conditions (2.10) are sufficient to give the three components of the flow velocity. The pressure distribution, if required, can be obtained from (2.6). This fully coupled and highly nonlinear system was solved using the spectral homotopy analysis method and the successive linearisation method. The results were validated using the Matlab bvp4c numerical routine and against results in the literature.

Following Boyd [19], we begin by transforming the domain of the problem from to using the domain truncation method. This approximates by the computational domain where is a fixed length that is taken to be larger than the thickness of the boundary layer. The interval is then transformed to the domain using the algebraic mapping

For convenience we make the boundary conditions homogeneous by applying the transformations

where and are chosen so as to satisfy boundary conditions (2.10). The chain rule gives

Substituting (3.2) and (3.3)-(3.4) in the governing equations gives

where prime denotes derivative with respect to and

As initial guesses we employ the exponentially decaying functions used by Yang and Liao [12], namely,

The initial solution is obtained by solving the linear parts of (3.5), namely,

subject to

The system (3.8)-(3.9) is solved using the Chebyshev pseudospectral method where the unknown functions and are approximated as truncated series of Chebyshev polynomials of the form

where and are the th Chebyshev polynomials with coefficients and , respectively, are Gauss-Lobatto collocation points defined by

and is the number of collocation points. Derivatives of the functions and at the collocation points are represented as

where is the order of differentiation and is the Chebyshev spectral differentiation matrix (see, e.g., [20, 21]). Substituting (3.10)–(3.12) in (3.8)-(3.9) yields

subject to the boundary conditions

where

The superscript denotes the transpose, is a diagonal matrix, and is an identity matrix of size . We implement boundary conditions (3.14) in rows 1, , and of in columns 1 through to by setting all entries in the remaining columns to be zero. The second set (3.15) is implemented in rows and , respectively, by setting , and setting all other columns to be zero. We further set entries of in rows , , , , and to zero.

The values of are determined from the equation

which provides the initial approximation for the solution of (3.5).

We now seek the approximate solutions of (3.5) by first defining the following linear operators:

where is the embedding parameter and and are unknown functions. The zeroth-order deformation equations are given by

where is the nonzero convergence controlling auxiliary parameter and and are nonlinear operators given by

The th-order deformation equations are given by

subject to the boundary conditions

where

Applying the Chebyshev pseudospectral transformation to (3.21)–(3.23) gives

subject to the boundary conditions

where and are as defined in (3.16) and

Boundary conditions (3.26) are implemented in matrix on the left-hand side of (3.25) in rows , , , , and , respectively, as before with the initial solution above. The corresponding rows, all columns, of on the right-hand side of (3.25), and are all set to be zero. This results in the following recursive formula for :

The matrix is the matrix on the right-hand side of (3.25) but with the boundary conditions incorporated by setting the first, , , , and , rows and columns to zero. Thus, starting from the initial approximation, which is obtained from (3.17), higher-order approximations for can be obtained through recursive formula (3.28).

The spectral homotopy analysis method, just like the original HAM, depends for its convergence rate on the careful selection of an embedded arbitrary parameter . Turkyilmazoglu [13] showed that the solution of the von Kármán problem by the homotopy analysis method is prone to wild oscillations when suction/injection is present. In this section we apply the successive linearisation method that requires no artificial parameters to control convergence to solve the governing equations (2.9)-(2.10). The method assumes that the unknown functions and can be expanded as

where , are unknown functions and and () are approximations that are obtained by recursively solving the linear part of the equation system that results from substituting (4.1) in the governing equations (2.9)-(2.10). Substituting (4.1) in the governing equations gives

where the coefficient parameters , (), , and are defined as

To facilitate direct comparison of the methods, we use the same initial approximations as in the case of the spectral homotopy analysis method of Yang and Liao [12]:

The solutions for , , , are obtained by successively solving the linearized form of (4.2), namely,

subject to the boundary conditions

Once each , () has been found, the approximate solutions for and are obtained as

where is the order of the SLM approximation. In coming up with (4.7), we have assumed that

Equations (4.5)-(4.6) can be solved using analytical techniques (whenever possible) or any numerical method. In this work the equations were solved using the Chebyshev spectral collocation method in the manner described in the previous section. This leads to the matrix equation

where is a square matrix and and are column vectors defined by

with

In the above definitions, , () are diagonal matrices of size and with being the Chebyshev spectral differentiation matrix. After modifying the matrix system (4.9) to incorporate boundary conditions, the solution is obtained as

The study of the magnetohydrodynamic swirling boundary layer flow over a rotating disk with suction or injection through the porous surface of the disk has recently been undertaken by Turkyilmazoglu [13]. In this case the Navier-Stokes equations reduce to a set of ordinary differential equations

subject to the boundary conditions

where is the magnetic interaction parameter due to the externally applied magnetic field and denotes uniform suction () or blowing () through the surface of the disk.

Turkyilmazoglu [13] utilized a twin strategy, using Ackroyd's series expansion and the homotopy analysis method to find purely analytic solutions to (5.1)–(5.5). In this study we use the SLM to obtain solutions to this system of equations.

Eliminating in (5.1) and (5.2) using (5.4) gives the following system of equations:

subject to the boundary conditions

The SLM is applied to (5.6) to (5.8) in the manner described in Section 4, and for brevity we omit the finer details. The intrinsic parameters of the SLM are essentially the same as those defined in Section 4 except for the following:

An appropriate initial approximation for finding in this case is

In this section we present the results for the velocity distributions and . To check the accuracy of the successive linearisation method and the spectral homotopy analysis method, comparison is made with numerical solutions obtained using the Matlab routine, which is an adaptive Lobatto quadrature scheme (see [22]). The current results are compared with previously published results by Liao [11], Yang and Liao [12], and Turkyilmazoglu [13]. The results presented in this work were generated using mostly collocation points and .

Table 1 gives a comparison of the values of obtained at different orders of the SLM and the SHAM approximations against the homotopy analysis method results, the homotopy-Padé results, and the numerical results. Our finding is that the SLM results converge most rapidly to the numerical result of 0.884474. Full convergence is achieved at the very low third order. Comparatively, convergence (to 6 decimal places) was achieved at the twentieth order using the homotopy analysis method and at the fifteenth order in the case of the homotopy-Padé method. When the same value is used, convergence of the spectral homotopy analysis method is achieved at the eighth order compared to the twentieth order for the homotopy analysis method approximations. This suggests that the SLM is a very useful computational tool that converges much more rapidly than the homotopy analysis method, the homotopy-Padé method, and the spectral homotopy analysis method, although, the SLM may, in fact, require more computations per step than the other methods.

Table 2 gives a comparison of the pressure difference at different orders of the homotopy analysis method, SHAM, and SLM against the numerical results. A similar pattern as in Table 1 emerges where the SLM results converge rapidly to the numerical result of with full convergence achieved at the third order. In the case of the HAM, convergence up to four decimal places was achieved at the tenth order. For the same values, the SHAM converges at the sixth order.

Tables 3–6 give a comparison between the SLM and the results reported by Turkyilmazoglu [13] for several suction/injection velocities and magnetic parameter values. Comparison of the results of Turkyilmazoglu [13] with the SLM seems most appropriate since the former study also partly utilizes a linearizing technique, the Newton-Raphson method to compute elements of the solutions. Turkyilmazoglu [13] showed that for large injection velocities, the number of terms required to attain convergence of the series solution increases dramatically, for instance, for injection velocities , up to 2000 terms are required to achieve convergence of the series solution method, and hence the study resorts to the Chebyshev collocation method to solve the governing equations. Nonetheless, our findings indicate that with only a few terms of the SLM series good levels of accuracy are achieved for all suction and injection velocities. For the suction and injection velocities in the range and in Tables 3-4 there is an excellent agreement between the fourth-order SLM, the numerical, and the results reported by Turkyilmazoglu [13].

Table 5 gives a comparison between the numerical and the SLM results for larger values of , up to when . Moderate increases in the suction/injection velocities appear to have no effect on the precision of the method with convergence again achieved at the fourth order of the SLM series. In Table 6, is fixed and the magnetic parameter varied. We compare the convergence rate of the SLM to the numerical computations and show that increasing this parameter has no effect either on the convergence rate of the successive linearisation method.

Figure 1 gives a comparison between the fourth-order SHAM, second-order SLM, and numerical results for the dimensionless velocity distributions and , respectively. There is an excellent agreement among the three sets of results. For purposes of comparison, it is worth noting that in case of the HAM in the work of Yang and Liao [12], agreement between the numerical and the HAM results was only observed at the 30th order of approximation for and at the 10th order for . As with most iterative methods, it is worth noting that the convergence rate may depend on the initial approximation used. However, since we have used the same initial approximations as Yang and Liao [12], the graphical results suggest that the SLM converges much more rapidly than both the HAM and SHAM. This may, however, be offset by the fact that the SLM may require more computations per step than the other two methods.

In this work two relatively new methods, the spectral homotopy analysis method and the successive linearisation method, have been successfully used to solve the von Kármán nonlinear equations for swirling flow with and without suction/injection across the disk walls and an applied magnetic field. The velocity components were compared with numerical results generated using the built-in Matlab bvp4c solver and against the homotopy analysis method and homotopy-Padé results in the literature. The results indicate that both the spectral homotopy analysis method and the successive linearisation method may give accurate and convergent results using only few solution terms compared with the homotopy analysis method and the Homotopy-Padé methods. Comparison has also been made with the recent findings by Turkyilmazoglu [13]. The successive linearisation method gives better accuracy at lower orders than the spectral homotopy analysis method. The tradeoff, however, is that both the spectral homotopy analysis method and the successive linearisation method may involve more computations per step compared to the methods in the literature.

Nonetheless, the successive linearisation method has been shown to be very efficient in that it rapidly converges to the numerical results. The study by Turkyilmazoglu [13] shows that whenever suction/blowing through the disk walls is present, the homotopy analysis method is prone to give wildly oscillating solutions. These oscillations have to be controlled by a careful choice of the embedded parameter . The advantage of the successive linearisation method is that such a parameter does not exist and no such oscillations are observed in the solution of the von Kármán equations for swirling flow.

The authors wish to acknowledge financial support from the National Research Foundation (NRF).

1. von Kármán, T: Uberlaminare und turbulence reibung. Zeitschrift für Angewandte Mathematik und Mechanik. 52(1), 233–252 (1921). PubMed Abstract

2. Cochran, WG: The flow due to a rotating disc. Proceedings of the Cambridge Philisophical Society. 30, 365–375 (1934). Publisher Full Text

3. Chien, C-S, Shih, Y-T: A cubic Hermite finite element-continuation method for numerical solutions of the von Kármán equations. Applied Mathematics and Computation. 209(2), 356–368 (2009). Publisher Full Text

4. Kirby, RM, Yosibash, Z: Solution of von-Kármán dynamic non-linear plate equations using a pseudo-spectral method. Computer Methods in Applied Mechanics and Engineering. 193(6–8), 575–599 (2004)

5. Yosibash, Z, Kirby, RM, Gottlieb, D: Collocation methods for the solution of von-Kármán dynamic non-linear plate systems. Journal of Computational Physics. 200(2), 432–461 (2004). Publisher Full Text

6. Zerarka, A, Nine, B: Solutions of the von Kàrmàn equations via the non-variational Galerkin–spline approach. Communications in Nonlinear Science and Numerical Simulation. 13(10), 2320–2327 (2008). Publisher Full Text

7. Adomian, G: A review of the decomposition method and some recent results for nonlinear equations. Computers & Mathematics with Applications. 21(5), 101–127 (1991). PubMed Abstract | Publisher Full Text

8. Lyapunov, AM: The General Problem of the Stability of Motion,p. x+270. Taylor & Francis, London, UK (1992)

9. He, J-H: Homotopy perturbation method for solving boundary value problems. Physics Letters A. 350(1-2), 87–88 (2006). Publisher Full Text

10. He, J-H: A coupling method of a homotopy technique and a perturbation technique for non-linear problems. International Journal of Non-Linear Mechanics. 35(1), 37–43 (2000). Publisher Full Text

11. Liao, S: Beyond Perturbation. Introduction to the Homotopy Analysis Method, CRC Series: Modern Mechanics and Mathematics,p. xii+322. Chapman & Hall/CRC, Boca Raton, Fla, USA (2004)

12. Yang, C, Liao, S: On the explicit, purely analytic solution of von Kármán swirling viscous flow. Communications in Nonlinear Science and Numerical Simulation. 11(1), 83–93 (2006). Publisher Full Text

13. Turkyilmazoglu, M: Purely analytic solutions of magnetohydrodynamic swirling boundary layer flow over a porous rotating disk. Computers and Fluids. 39(5), 793–799 (2010). Publisher Full Text

14. Ackroyd, JAD: On the steady flow produced by a rotating disc with either surface suction or injection. Journal of Engineering Physics. 12, 207–220 (1978)

15. Motsa, SS, Sibanda, P, Shateyi, S: A new spectral-homotopy analysis method for solving a nonlinear second order BVP. Communications in Nonlinear Science and Numerical Simulation. 15(9), 2293–2302 (2010). Publisher Full Text

16. Motsa, SS, Sibanda, P, Awad, FG, Shateyi, S: A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem. Computers and Fluids. 39(7), 1219–1225 (2010). Publisher Full Text

17. Makukula, Z, Motsa, SS, Sibanda, P: On a new solution for the viscoelastic squeezing flow between two parallel plates. Journal of Advanced Research in Applied Mathematics. 2(4), 31–38 (2010). Publisher Full Text

18. Motsa, SS, Sibanda, P: A new algorithm for solving singular IVPs of Lane-Emden type. In: Proceedings of the 4th International Conference on Applied Mathematics, Simulation, Modelling, July 2010, Corfu Island, Greece, NAUN International Conferences. 176–180

19. Boyd, JP: Chebyshev and Fourier Spectral Methods,p. xvi+668. Dover, Mineola, NY, USA (2001)

20. Canuto, C, Hussaini, MY, Quarteroni, A, Zang, TA: Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics,p. xiv+557. Springer, New York, NY, USA (1988)

21. Don, WS, Solomonoff, A: Accuracy and speed in computing the Chebyshev Collocation Derivative. SIAM Journal on Scientific Computing. 16(6), 1253–1268 (1995). Publisher Full Text

22. Kierzenka, J, Shampine, LF: A BVP solver based on residual control and the MATLAB PSE. ACM Transactions on Mathematical Software. 27(3), 299–316 (2001). Publisher Full Text