Solving fluid dynamics problems with Matlab

Abstract

MATLAB (short for Matrix Laboratory) was created by Cleve Moler and Jack Little in the 1970’s. It is a programming language for technical computing. Its environment is easy to work with, the syntax is very simple and intuitive, it has powerful toolboxes to treat many different problems in engineering, and it allows us to produce fantastic graphics as the programme runs. It also allows us to create a graphical interface (via graphical user interfaces - GUIs) that gives our programme a look that is very close to professional software. Because of many of the mentioned features, a MATLAB code can be very compact, allowing anyone to have "the big picture" of any code without have to look at all its details. Another great advantage of Matlab is that, if the code is written in a vectorized form, the code can run much faster than if it was written in the traditional form (’a la C/fortran’). The fact that MATLAB allows us to use a powerful toolbox for sparse matrices, is also a great advantage since, many traditional linear algebra operations can be highly improved, allowing the codes to run much faster than it would run with the traditional linear algebra functions. In our work we have made extensive use of MATLAB to do ’proof of concept’ studies, especially when developing new algorithms and techniques for solving systems of coupled nonlinear partial differential equations, such as those which arise in fluid dynamics. This includes, for instance, codes for investigating instabilities in lid-driven cavities, Boppana and Gajjar (2010), instabilities in flow past circular cylinders, Boppana and Gajjar (2011), and transonic flow past aerofoils Pereira and Gajjar (2010). In some cases MATLAB is used in its own right for solving small problems, but the fact that MATLAB is an interpreted language means that for increasing problem sizes, the MATLAB version of the code can be much slower than equivalent versions in other languages especially when one is dealing with very large sparse matrices. On the other hand the beauty of MATLAB is that much of the hard work is buried in the simple syntax and hidden from the user. An example of this is the use of the backslash operator for solving linear systems. Whether the system is sparse or full, the manner in which the equations are solved is hidden from the user and this greatly facilitates code development. In the equivalent fortran versions of the code the replacement for the ’\’ operation requires considerable work and the code translation process is no longer a trivial exercise. In this chapter we will discuss the use of hybrid spectral methods to solve two and three-dimensional problems using MATLAB. There is an excellent book by Trefethen (2000) which discusses the application of spectral methods using MATLAB to solve ordinary and partial differential equations, and which provides the foundation for the techniques described below. To motivate the ideas we first consider the solution of a model elliptic equation of the form a(x, y)ψxx + b(x, y)ψyy + c(x, y)ψx + e(x, y)ψy + f (x, y)ψ = g(x, y), say with Dirichlet boundary conditions in a rectangular domain. To obtain a numerical solution to this problem the first step is to choose an appropriate method and discretization. In our work we have used a combination of spectral methods in one or two dimensions and high order finite difference methods in another dimension. The main reasons for this choice are that a hybrid approach combines the accuracy of spectral methods together with flexibility in comparison to using spectral methods on their own. One restriction with the use of spectral methods is that one needs to use somewhat simplified geometries. A hybrid approach gives more flexibility in this respect. Another important reason stems from consideration of the type of matrix patterns which arise. With finite differences in say the x direction and spectral collocation in the y direction, the coefficient matrix with the unknowns ordered in terms of increasing y values for a fixed x value, has a particular sparsity pattern dependent on the order of the finite-differencing used. Second order finite-differencing leads to a a block tridiagonal matrix whilst with fourth-order finite differences, the matrix is block-pentadiagonal of the form: AqΨq−2 + BqΨq−1 + CqΨq + DqΨq+1 + EqΨq+2 = Rq, q = 0, 1, ...,M. (1) Here Ψq is the vector of unknowns at location x = xq, M+1 is typically the number of points in the x direction and the block matrices are of size N +1 by N +1 where N +1 spectral points are used. UsingMATLAB it is not too difficult to generate a short code to solve the above discrete system and the book by Trefethen (2000) gives plenty of such examples. The problem becomes more challenging when N and M become large, as for example in some fluid flow applications where large N,M values are needed to resolve regions of the flow where the solution changes very rapidly. When using a large number of points the sparse matrix facilities of MATLAB come into their own. The whole coefficient matrix does not need to be stored and by declaring this as a sparse matrix, only the non-zero entries of the block matrices are calculated. This avoids having to store a very large and sparse matrix which can quickly lead to memory problems. Increasing the order of the scheme leads to increased bandwidths. This sparsity pattern can be exploited for 2nd, 4th or even 6th finite-differencing with a direct solver. However, with spectral methods in two directions, unless the differential operator involved has a special form, it is not immediately possible to utilize the sparse nature of the matrix. Whilst this does not pose any intrinsic difficulties if one is coding in MATLAB, with increased number of points the solution phase can become very memory intensive and requires a lot of processor time. The use of the hybrid approach in our work is motivated in part by the observation that the sparse matrix structure can be exploited to write efficient solvers, which not only work well with MATLAB, but can be coded directly in other languages. MATLAB provides for an excellent environment in which one can test and develop solvers of this type. The above techniques have been successfully applied to investigate a whole range of different flow problems governed by the Navier-Stokes and related equations. In the first example we consider the onset of instability in the lid-driven cavity flow. MATLAB was used to generate results on coarse grids and do preliminary eigenvalue computations. For very fine grids, the computations were performed in Fortran 95. The problem is described in detail in Boppana and Gajjar (2010). The second problem concerns the onset of instability in the flow past a row of circular cylinders. Again the same technqiues have been used but for a more complicated geometry. This problem is described in detail in Boppana and Gajjar (2011). The third problem we discuss concerns the inviscid transonic flows past thin airfoils. Here the governing equations are nonlinear and of mixed type and the flow can contain shock wave discontinuities for certain parameter values. The full details are given in Pereira and Gajjar (2010). The same methods as described above are used except now type differencing needs to be incorporated to allow for the different flow behaviours in regions of subsonic and supersonic flow. The method is fast and very robust and we are able to compute steady flows with strong shocks. The code was written in MATLAB, using vectorization when possible, and, in order to produce a good interface with the user, we used GUIs (graphical user interfaces) from MATLAB. The result was a fast and accurate code, with the extra bonus of a very good interface with the user, without a lot of effort in terms of programming. The fact that graphical results can be shown immediately, saves us a lot of work, both on the analysis of the results and on its presentation.