Solving Fluid Dynamics Problems with Matlab

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 12


Introduction
MATLAB (short for Matrix Laboratory) was created by Cleve Moler and Jack Little in the seventies.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 (2010a), instabilities in flow past circular cylinders, Boppana and Gajjar (2010b), 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 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: Here Ψ q is the vector of unknowns at location x = x q , 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.Using MATLAB 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 (2010a).
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 (2010b).
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.

Instabilities in lid-driven cavities.
To motivate the techniques used in our work we first consider a model elliptic equation of the form a(x, y) with say Dirichlet boundary conditions in a square domain 0 ≤ x, y ≤ 1.It is assumed that the functions a, b, c, e, f , g are smooth functions of x and y.
The discrete solution to the equations will be obtained at a set of grid points on an (M + 1) × (N + 1) grid say with the nodes, x = x j , j = 0, . . ., M with x 1 = 0, x M = 1, and y = y k , 0 ≤ k ≤ N with y 0 = 1, y N = 1, and at these points we set ψ j,k = ψ(x j , y k ).
We will assume that derivatives in the x− and y− directions may be approximated as Given the type and order of discretisation, the elements of the matrices D x , D xx , D y , D yy are known.For example with second-order central finite differences, at an interior point of a uniform grid in x with ∆ x = x j − x j−1 , we have and zero otherwise.If we take Chebychev collocation in the y−direction, the D y , D yy are the Chebychev differentiation matrices and are given in a number of places, see for example Weideman and Reddy (2003) or Trefethen (2000), where MATLAB code to generate the matrices is given.Discretization of the equation 2 leads to the set of equations At the boundaries we have ψ = 0.The discrete set (4) can be combined into a compact form as where Ψ denotes the vector of unknowns, A ⊕ B is the kronecker tensor product of two matrices A, B (which in MATLAB is represented by the kron(A,B) operator), and DIAG(v) is as in MATLAB the diagonal matrix with entries given by the vector v, and I M+1 is the identity matrix of order M + 1. Certain rows of the matrix operator in ( 5) are also modified when the boundary conditions are incorporated.
The linear sytem may be represented as Once the coefficient matrix S and the right hand sides have been computed, the solution just involves the use of the \ operator with Ψ = S\R.From the user perspective other than declaring that the matrices involved are sparse matrices, no additional special treatment is required to obtain the solution of the linear systems.
Given the discretization matrices, the above system is easy to code in MATLAB.For certain discretizations however, the linear systems outlined above can be huge and highly sparse.
The bandwidth of the coefficient matrix increases with increasing order of differences used, and with spectral methods in two directions, the coefficient matrix has the sparsity pattern similar to that in figure 1(c), obtained using the spy function in MATLAB.On the other hand by taking second-order finite differences in the x−direction and spectral collocation in the y−direction, with a particular ordering of grid points, the matrices can be written in block tridiagonal form as shown in figure 1(a),1(b), or with fourth order finite-differencing in x the linear system is block pentadiagonal.With 2nd, 4th or even 6th finite-differencing the linear systems can be solved with a direct solver but 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.In our work we have written our own direct block solvers as well as making direct use of the \ operator with a sparse matrix.
In the MATLAB codes when constructing the coefficient matrices, we use the functions speye or spalloc for initially creating the sparse matrices.After this only the non-zero elements of the matrices are assigned.

Lid-driven cavity flow.
In this section we discuss the steps needed to go from the continuous model described by the set of coupled nonlinear partial differential equations, the Navier-Stokes equations, to solution using MATLAB.To focus attention we will consider the classic problem of flow in a lid-driven cavity which has been extensively studied in the literature.The governing equations for the problem are ∇ 2 ψ = ω, and Here Re is the Reynolds number defined as Uw ν , where U is the velocity of the lid, ν is the kinematic viscosity of the fluid, and w is the width of the cavity that are used to non-dimensionalize the velocity and length-scale variables respectively.The boundary conditions (see figure 2(a)) are given by where A is the aspect ratio of the cavity.The discrete solution to the equations will be obtained at a set of grid points on an (M + 1) × (N + 1) grid say with the nodes, x = x j , j = 0, . . ., M with x 1 = 0, x M = 1, and y = y k , 0 ≤ k ≤ N with y 0 = 1, y N = 1, and at these points we set We will assume that derivatives in the x− and y− directions may be approximated as in (3).
The discretization leads to a set of nonlinear algebraic equations for the unknowns To solve these we use Newton linearization and this will lead to linear system of equations which are solved using MATLAB.Linearization of the nonlinear equations is performed by setting ψ j,k = ψj,k + G j,k , ω j,k = ωj,k + H j,k , the barred quantities representing a current iterate and G j,k , H j,k begin small corrections, and substituting into the nonlinear equations and neglecting second order small terms.This leads to Here G, H are the vector of unknown corrections.Depending on the discretization, the above can be coded directly in MATLAB by constructing the coefficient matrix multiplying the vector of unknowns (G, H) T .Note that the size of the coefficient matrix is 2(N + 1)(M + 1) × 2(N + 1)(M + 1) and even for modest N, M the above procedure leads to very large matrices and is not efficient.The approach we have adopted is to make use of the sparsity patterns for particular types of discretizations.

Use of Matlab for the solution of the discrete system
For the case when we use second order finite-differeces in the x−direction and chebychev collocation in the y−direction, the linear system may be written as SΦ = (S 0 , S 1 , . . ., S M ) T = r, where This represents a particular ordering of unknowns and gives rise to the block tridiagonal system in (9).With 4th order finite-differencing in x the linear system is block pentadiagonal.The coefficient matrices A p , B P , C p can be extracted from the discrete equations above and with The above excludes the boundary conditions, but these just alter certain rows of the matrices.
In MATLAB the individual entries of the block matrices are easily computed and the S matrix is updated via

Results for lid-driven cavity
The method described above was used to compute the flow in a lid-driven cavity.The same techniques used were adapted to firstly compute the steady flow, and then investigate the instability of the flow via simulations as well as solving the linear eigenvalue problem assuming normal mode disturbances proportional to e λt .Further details of the techniques may be found in Boppana and Gajjar (2010a).MATLAB was also used in the eigenvalue analysis.In fact the eigenvalue problem required to be solved takes the form which is a generalised eigenvalue problem.The matrix S is the same as the Jacobian matrix of the linear system after Newton linearization, and T is a singular diagonal matrix.The eigenvalue problem is solved in MATLAB using the eig function.In other related problems the routine sptarn available in the PDE toolbox was used for the solution of the generalised eigenvalue problem.
In figure 3 we have shown the real and imaginary parts of eigenfunctions for the disturbance streamfunction ψ for aspect ratios of A = 1 and A = 2 at the onset of instability, obtained from a solution of the eigenvalue problem.The advantage of working with MATLAB is that such plots can be generated at the same time as the computation is in progress.More extensive results and details are documented in Boppana and Gajjar (2010a).

Flow past circular cylinders
The techniques described above have also been used to solve for the uniform flow past a row of circular cylinders, see figure 2(b).Here the Navier-Stokes equations are solved with boundary conditions of no slip on the cylinder surfaces and uniform flow far upstream.The two important parameters are the gap-width between the cylinder centres W (normalised with respect to cylinder radius) and the Reynolds number Re.Of particular interest is to ascertain when the flow first becomes unstable and the mode of instability.Experiments, see for example Mizushima & Ino (2008), for the flow past two such cylinders show a complicated dynamics as the parameters are varied.For large Reynolds numbers, and for large gap widths the flow is like that past an isolated cylinder and the observed shedding frequencies also similar.For Re >> 1 and intermediate gap widths (of 0.5 to 1 cylinder diameters), the flow can become deflected to one side and becomes asymmetric.For small gap widths, the flow is similar to that past a single bluff body.At low Reynolds numbers the flow dynamics is also complicated and with conflicting results.
In our work, MATLAB was used to first compute the base flows, and then for studies of the onset of instability, by solving the generalised eigenvalue problem.Full details of the numerical methods and results can be found in Boppana and Gajjar (2010b).In figures 5, 6 we have shown the results for the streamfunction eigenfunctions.Modes which are symmetric and asymmetric with respect to the cylinder centreline, can be observed and it found that the critical Reynolds numbers for the onset of instabilties are very similar for both modes.However the symmetric modes as shown in 5(b), 6(b) for instance, have much lower critical Strohal frequencies as compared to asymmetric mode.The latter is the one which tends to that for the flow past an isolated cylinder as the gap width increases.

Transonic flows past airfoils.
The study of transonic flows is motivated in part by the observation that many modern airplane carriers operate most efficiently when cruising at speeds which fall in the transonic range, that is close to the speed of sound.Mathematically the study of transonic flows is fascinating because the governing equations are nonlinear and of mixed type.The methods

The mathematical model
In this subsection we describe briefly a numerical method to deal with transonic flows past thin airfoils, based on the Kárman-Guderley equations.This method was discussed in detail in Pereira and Gajjar (2010).However the aim of this section in not to repeat once again the same work, but rather to focus on the use of the graphics interface which is extremely easy to program in MATLAB.First we will give an overview of how the mathematical model was built, and then we will explain how to use MATLAB graphics user interfaces (GUIs) in the context of our problem.
The governing equation may be obtained using asymptotic methods as described in Cole and Cook (1986).The starting point to our model was the full potential equation: where Φ represents the velocity potential, a is the local speed of sound, U ∞ is the velocity in the far field, a ∞ is the speed of sound in the far field, and M ∞ = U ∞ /a ∞ is the free-stream Mach number.The velocity components (U, V) are defined as follows, The density ρ and pressure p can be determined via the relationships, where γ is the ratio of specific heats.
In order to define the boundary conditions we assume that the flow is uniform in the far field and that the flow is tangent to the airfoil on its surface.To construct this theory we also assume that we have a thin aerofoil with width (δ → 0) and that the air flow speed is close to sonic so M 2 ∞ = 1 − kµ(δ), µ(δ) → 0 where k is the transonic similarity parameter and µ is a function of the airfoil width (δ).The oncoming flow is assumed to be aligned with the x-direction.In order to define the boundary conditions we assume that the flow is uniform in the far field and that the flow is tangent to the airfoil on its surface.To construct this theory we also assume that we have a thin aerofoil with width (δ → 0) and that the air flow speed is close to sonic so M 2 ∞ = 1 − kµ(δ), µ(δ) → 0 where k is the transonic similarity parameter and µ is a function of the airfoil width (δ).The oncoming flow is assumed to be aligned with the x-direction.The airfoil is defined by, y = δF(x).
Introducing non dimensional variables, then the full potential equation becomes, The expansion for Φ is described in Cole and Cook (1986) and is given by, It is well known that as M ∞ → 1, the perturbations extend in the y direction significantly.Because of this, stretched coordinates were used and the governing equations and boundary conditions reduce to, where ( 12) is the so called Kárman-Guderley equation.Equation ( 12) in conservative form is written as follows, ∂ ∂x and, if we denote, ψ = kx − (γ + 1)φ ( 16) then, we may rewrite (15) as, The boundary conditions become, When considering the non symmetric case, one has to introduce a new boundary conditionthe Kutta condition.The Kutta condition is used at the trailing edge to ensure that the jump obtained in the integration of ψ x along the lower and upper surfaces at the tail is zero.
Next we give a brief overview of the numerical method used to solve the above problem.Ee used finite differences for the derivatives in the x direction and a Chebyshev collocation method to describe the derivatives in the Y direction.As described in Cole and Cook (1986) each point of the domain may be either subsonic, sonic, supersonic or a shock point.Let, then equation ( 17) becomes, Let (ψ x ) i+1/2,j represent the derivative with respect to x of ψ at the point (x i+1/2 , Y j ), and let h i = x i − x i−1 .Using central differences for the x derivatives we may write (19) as, or, As discussed in Cole and Cook (1986) if we have a subsonic point, the equation is elliptic and central differences should be used to calculate both (ψ x ) i+1/2,j and (ψ x ) i−1/2,j .If we have a supersonic point, then equation ( 20) becomes hyperbolic and backwards differences should be used to calculate both (ψ x ) i+1/2,j and (ψ x ) i−1/2,j .In the first case we can rewrite (20) as, In the second case, we can rewrite (20) as, where, and, Two other cases are considered, the case when the flow accelerates from subsonic to supersonic in which case we define a sonic point, and the opposite, when the flow decelerates form supersonic to subsonic in which case we call it a shock point.To deal with these cases we use artificial viscosity (µ i,j ) Murman and Cole (1971), and equations ( 21) and ( 22) can be condensed as, where the values of µ i,j and µ i−1,j are taken from the following table: Note that a shock point can be seen as an addition of both elliptic and hyperbolic x difference operators.
In the Y direction, the physical domain was first truncated to y ∞ and mapped into the Chebyshev space, as in Canuto et. all (1998), where, and, First and second derivatives in the Y directions were calculated as described earlier.
After applying the above discretizations we obtain a set of coupled nonlinear algebraic equations.These are linearized using Newton-Raphson linearization by setting where ψ i,j represents the value of ψ i,j in a previous iteration and G i,j represents the update for ψ i,j .This results in a linear system of equations for the G i,j of the form Details on how to calculate the coefficients A i,j , B i,j ..., F i,j can be found in Pereira and Gajjar (2010).The block pentadiagonal system of equations was solved directly using routines described in Korolev et all. (2002).

The use of GUIs
In this subsection we will focus on how to generate a good graphics interface using MATLAB GUIs in the context of the problem described in the previous subsection.In order to obtain a graphics interface to our programme, the first thing we have to do is to type guide in MATLAB's command window.The result is that MATLAB opens a window that has a graphics interface working environment(GIWE).The next thing to do, is to save it as name.fig.
This action has as a result not only of saving the work done so far, but also to generate a file name.m, that has the MATLAB corresponding instructions.Suppose next that we want to introduce a title on top of the graphics interface window.We select the "Static text" button in the GIWE.With the mouse one selects the location and the size of the text to input.Then after double clicking on it a new window appears, where one can choose options for the static text we want to introduce.The results of both examples agree with the literature, see for example Cole and Cook (1986) and Camilo (2003).As we can see, by using GUIs, we obtained a user friendly professional looking graphics interface, which allowed any other user to perform simulations and input their parameters for the run.The other users of the software did not need to be familiar with the underlying code and equations.

Conclusion
-The environment of MATLAB is easy to work, 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 program runs.
-A MATLAB code can be very compact, allowing anyone to have "the big picture" of any code without having 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.-The main drawback with the use of MATLAB is that when the number of points used in the discretization grows, which is necessary for finely resolved computations, the execution time and memory requirements can grow substantially.On multicore machines, the use of the parallel toolbox is recommended and this is supposed to provide superior performance making use of parallel linear algebra routines with minimal code changes.However we have not been able to test this with our codes.

Fig. 1 .
Fig. 1.Sparsity patterns for model problem with (a) second order finite-differences in both x and y with N = 10, M = 11, (b) second-order finite-difference in x and chebychev collocation in y with N = 10, M = 11, (b) chebychev collocation in both x and y with N = 10, M = 11.

Fig. 2 .
Fig. 2. Sketch of (a) the lid-driven cavity with boundary conditions, and (b) flow past a row of circular cylinders.

Fig. 3 .
Fig. 3. Eigenfunctions of the streamfunction for lid-driven cavity at a critical Reynolds number of Re = 8026.7 and aspect ratio A = 1.

Fig. 4 .
Fig. 4. Eigenfunctions of the streamfunction for lid-driven cavity at a critical Reynolds number of Re = 5861 and aspect ratio A = 2.

Fig. 5 .
Fig. 5. Perturbation eigenfunctions the flow past a circular cylinder for a gap width W = 5 with (a) anti-phase oscillatory mode, (b) in phase oscillatory mode.

Fig. 6 .
Fig. 6.Perturbation eigenfunctions the flow past a circular cylinder for a gap width W = 50 with (a) anti-phase oscillatory mode, (b) in phase oscillatory mode.described above work have been applied to partial differential equations of mixed type containing shocks.