INTRODUCTION


 

 

 

16 
 

 American Scientific Research Journal for Engineering, Technology,  and Sciences  (ASRJETS) 

ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 

© Global Society of Scientific Research and Researchers  

http://asrjetsjournal.org/  

 

Consistency, Stability and Convergence of the Viscous 

Garabedian and Korn (VGK) Method 

Dr. Salah.J.S. El-Ibrahim* 

Lecturer, College of Aviation Technology, P.O. Box 3357, Hawally 32034, Kuwait 

Email: s.alhilat@collegeofaviation.com, Salah13@gmail.com 

 

Abstract 

The Viscous Garabedian and Korn (VGK)
 
is a computational fluid dynamics (CFD) code, used for the 

prediction of viscous flows around two dimensional aerofoil sections (mapped onto a circular grid) in a subsonic 

free stream. The code is written in “Fortran 77” and solves the full potential flow equations using a finite 

difference algorithm. The code has been recently published commercially with a data preparation routine to 

facilitate its use
 
[1]. The present work investigates the code’s consistency, stability and convergence in the hope 

that it will shed light on the code’s performance and limitations under various flow conditions and code’s 

numerical parameters.  

Keywords: Artificial Viscosity; Backward Differencing Scheme; Central Differencing Scheme; Conformal 

Mapping; Consistency; Convergence; Grid Resolution; Numerical Oscillations; Relaxation Factor; Stability; 

Viscous Garabedian and Korn (VGK). 

1. Introduction 

Analytical solutions to the governing flow equations (the full Navier-Stokes equations) do not exist for the 

general case. Hence solutions must be obtained either experimentally or computationally. Computational fluid 

dynamics (CFD) involves the numerical solution of the flow equations governing fluid motion. The main uses of 

CFD in aeronautical applications fall under two main categories [2].
 
First, the provision of reliable aerodynamic 

predictions, enabling designers to produce “better aeroplanes”. Second, there is the possibility of using CFD for 

purely scientific investigations. It seems possible that numerical simulation of complex flows, not readily 

accessible to experimental measurements, can provide new insights into the underlying physical processes. In 

particular, computational methods offer a new tool for the study of turbulent structures and the mechanism of 

transition from laminar to turbulent flows. Cavity flows and surface irregularities may also fall under such 

complex flow categories, particularly in three-dimensional viscous flows.  

------------------------------------------------------------------------ 

* Corresponding author.  

 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

17 
 

Reflecting on the different physical behaviour of flow fields, most computational techniques and algorithms are 

tailored for a particular flow field [3], hence they tend to have their own (specific) applications and inherent 

limitations even within the intended application.Accordingly, any computational solution must be verified either 

experimentally or through comparison with other previously verified solutions to assess its validity and 

limitations. Once such criteria are established the use of the computational solution can be justified. 

Experimental work, therefore, -although time consuming and expensive on resources- is still a necessity, 

particularly if no other previously validated solution exists to compare with the one under investigation. In 

addition to the validation/verification process, the computational solution must also be assessed for the 

reliability and quality of its numerical solution. This can be done by analysing the Stability, Consistency, and 

Convergence of the algorithm(s) used in the numerical calculations. The present work attempts to investigate 

these aspects (Stability, Consistency and Convergence) of the VGK method and to describe the factors affecting 

its solutions under various flow conditions. 

2. The VGK Method 

The Viscous Garabedian and Korn (VGK)
 
is a computational fluid dynamics (CFD) code, used for the 

prediction of viscous flows around two dimensional aerofoil sections (mapped onto a circular grid) in a subsonic 

free stream. The code is written in “Fortran 77” and solves the full potential flow equations using a finite 

difference algorithm. The code was published commercially with a data preparation routine to facilitate its use 

[1].  Because of its relative simplicity in use and computational economy, VGK offers a convenient method for 

the assessment of the effects of small deviations from a specified aerofoil profile, such as may be caused by 

damage or manufacturing errors [4]. The reliable assessment of such imperfections enables a decision to be 

made as to whether to tolerate, repair, or reject the surface under investigation. In the aircraft industry the 

decision involves three important considerations: flight safety, money and time. However, to reach a decision 

the assessment process itself may be both costly and lengthy and in the end may not be fully reliable. This 

highlights the importance of relying on a code that will provide quick and reliable answers to the decision 

maker. VGK offers simplicity of use and had been proved [5] to provide valid solutions to high speed subsonic 

flows (up to the point where weak shock waves may exist) around aerofoils with surface imperfections provided 

that these imperfections do not sharply change the curvature of the profile, also under flow conditions where 

weak shock waves can present in the flow. Therefore, for a simple and direct assessment method it is hoped that 

VGK will provide quick and reliable answers.  

2.1 VGK Capabilities and Limitations 

The Viscous Garabedian and Korn (VGK) code is a computational fluid dynamics code written in Fortran 77, 

and based on the coupled solution of the inviscid and viscous flow regions of two-dimensional transonic 

attached flow past an aerofoil in a subsonic free stream. The full potential equations governing the inviscid flow 

region, formulated on the assumptions of steady and compressible flow of air (with specific heats ratio  =1.4), 

are solved by a finite difference technique. The viscous flow region is solved using the integral equations 

governing the laminar and turbulent boundary layer components, and solutions involve the displacement and 

momentum thickness distributions of the boundary layer and wake, which are calculated allowing for the effects 

of the pressure distribution obtained from the inviscid flow solution. The boundary conditions for the inviscid 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

18 
 

flow solution specify a non-zero normal velocity at the aerofoil surface to allow for the growth of the boundary 

layer displacement surface, and a jump in velocity across the dividing streamline to allow for the wake thickness 

and curvature effects. This introduces the concept of “Equivalent Inviscid Flow-EIF”. This means that the 

“actual inviscid” flow region is replaced by an “equivalent inviscid flow region, that, by definition, extends right 

onto the wall (the aerofoil surface), and a convenient line in the wake dividing the flow of the upper and lower 

surfaces of the aerofoil. VGK employs a central-difference and a backward difference schemes (for fully 

subsonic and supersonic flow regions respectively) of the Taylor series expansion of the full potential flow 

equation (equation (1) below).  VGK uses computational grid in the solution which is an O-type grid. It is 

generated by conformally mapping the infinite region outside the aerofoil into the inside of a unit circle. The 

grid points are defined by the intersection of radial lines emanating from the aerofoil surface and concentric 

circles circumferentially placed around the aerofoil, so the grid is allowed to touch the aerofoil surface. The 

computational grid is calculated by the VGK code utilising the aerofoil surface co-ordinates. VGK requires that 

for the grid to be satisfactory, it is necessary (but not sufficient) that the input aerofoil co-ordinates must be 

closely pitched, specially near the leading and trailing edges, and also be sufficiently smooth to yield smooth 

first and second order derivatives. VGK employs an iterative scheme to calculate the computational grid using a 

spline interpolation of the input aerofoil co-ordinates for each iteration step, and with about 30 iteration steps 

convergence is rapid for most of the aerofoil geometries of interest. Grid sequencing is employed in VGK to 

reduce the truncation error (residuals) at all the grid nodes and only two grids are used in VGK for this purpose: 

a coarse grid constituting 80 radial lines x 15 circumferential circles, and a fine grid that constitutes 160 lines x 

30 circles. The VGK solution process starts with the evaluation of the computational grid. Once this is achieved, 

the iteration solution of the finite difference equations for the inviscid flow potential proceeds by a relaxation 

method in which the values of potential obtained at the end of each iteration step are multiplied by a factor 

(called a relaxation factor). When this factor is greater than unity the process is termed over relaxation and when 

under unity the term under-relaxation is used. The associated factors are called over relaxation or under-

relaxation factors respectively. Over- or under-relaxation is used to increase the rate of convergence or avoid 

divergence depending on the form of the finite difference equations. In VGK different relaxation factors are 

used for the finite difference equations corresponding to subsonic and supersonic flow regions and the default 

values of these factors are set at 1.9 and 1.0 respectively. A relaxation process is also used in introducing the 

boundary conditions corresponding to the viscous flow element into the inviscid flow. A number of inviscid 

flow iterations (typically 5) are carried out between each occasion the boundary conditions are updated, and the 

calculated changes to these boundary conditions are decreased by an under-relaxation factor (typically 0.15 for 

the fine grid and 0.075 for the coarse grid). The final coupled (viscous-inviscid) solution is considered to be 

converged when the fine grid residuals in both the viscous and inviscid equations have reduced to a level 

beyond which any changes, as a result of further iterations, will have zero or negligible effects on the calculated 

flow parameters. In addition to the basic pressure distribution over the aerofoil surface, VGK provides a number 

of local surface parameters, such as Mach number (M), pressure ratio (local static pressure/total stagnation), and 

pressure coefficient, and overall parameters such as lift, drag and pitching moment coefficients. Data related to 

boundary layers and wake are also output, these include: shape parameter, displacement and momentum 

thicknesses, skin friction coefficient and boundary layer thickness for the upper and lower surfaces boundary 

layer, and pressure coefficient, pressure ratio, and Mach number for the lower and upper parts of the wake. 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

19 
 

Certain “numerical control parameters” are used in VGK to control the grid density, rate of convergence and the 

form of the finite difference scheme employed in any supersonic flow region. These parameters influence the 

results obtained by any VGK run, hence they must be carefully chosen to generate the most reliable data. 

However, VGK uses default values for these parameters (known as “default parameters”) which are generally 

suited for most flow conditions provided no strong shock waves or separated boundary layers are present. For 

flows with moderate shock waves or boundary layer being close to separation these parameters may be changed 

to enable convergence and “reliable” data to be obtained. Table (1) gives a summary of these parameters and 

their related influence on the calculated flow. The VGK code has some limitations; for example, there is no 

method in VGK for predicting boundary layer transition location, & the user has to input the location for the 

upper and lower surfaces. If the calculated local skin friction coefficient becomes less than a limiting value 

(=2x10
-6

) separation of the laminar boundary layer is assumed to occur at this location and the transition point is 

moved to this location instead of the specified one. This may be thought of as a simplified model for a short 

bubble laminar separation/ turbulent re-attachment which would normally occur in practice [6].  There is also no 

check applied within VGK that the boundary layer could actually become turbulent at the specified transition 

points and remain turbulent thereafter. Solutions in VGK will not converge for difficult flow conditions 

associated with strong shock waves (Mshock >1.4) or regions of extensive boundary layer separation.  

Table 1: VGK Numerical parameters 

PARAMETER & 

DEFINITION 

VALUE EFFECTS OF CHANGES ON 

THE FLOW 
DEFAULT MAX. MIN. 

a). NSG, Number of the 

radial mesh lines in the 

fine grid
*
 

160 160 96 Very slight effects for flows 

without shocks, 

Lower values give less steep 

shock pressure rises, and milder 

effects on the boundary layer 

b). NSC, Number of 

iterations of the inviscid 

flow element using the 

coarse mesh 

100 May be 

double

d 

(see f 

below) 

0 Default value generally 

adequate to reduce the level of 

residuals for the commencing 

fine grid convergence 

c). NSF, Number of 

iterations of the inviscid 

flow elements using the 

fine mesh 

200 May be 

double

d 

(see f 

below) 

---- Default value generally provides 

adequate convergence level; a 

further few hundred iterations 

may be needed if the residuals 

are declining, or if the CL value 

has not stabilised after 200 

iterations 

d).XS, The subsonic flow 

relaxation parameter 

1.9 <2 >0 For difficult flows lower values 

may be needed to avoid 

divergence, values down to 1.6 

may confer advantages 

e).XM, The supersonic 

flow relaxation parameter 

1 1 > 0 Same as XS above, lower values 

down to 0.8 may confer 

advantages 

f).EP( ),Artificial 

viscosity parameter
**

, and 

QCP( ),Partially-

conservative parameter
*** 

 =0.8 

 =0.25 

  These two parameters jointly 

influence the computed 

supersonic flow. 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

20 
 

NOTES on table (1):  

* There is no facility in VGK for changing the number of circumferential lines in the fine grid. 

** The artificial viscosity parameter is  is added into the finite differencing equations to allow for the inclusion 

of higher order terms and reduce the discretisation error. 

*** The partially conservative parameter  is a shock point operator [13] , it is assigned any value between  0.0 

and 1.0 and is meant to simulate conservative  1 or non-conservative schemes  0 . Thus improving 

the pressure jump through a shock wave.  

However, for flows with moderate difficulties (1.2<Mshock<1.4), convergence may be achieved in steps, that is to 

arrive at the run under question as a continuation run to previous runs achieved with lower Mach number, lower 

incidence or higher Reynolds’ number. It has also been found that VGK solutions can exhibit numerical 

oscillations when large changes in profile curvature are encountered (aerofoils with indentations). The present 

work attempts to investigate the numerical behaviour of VGK under various flow conditions for the purposes of 

assessing its overall reliability and identifying the factors (in the numerical schemes) that may cause calculations 

issues in the hope of finding ways to improve the code’s performance. 

2.2 Limitations/Constraints of the present work 

For brevity, the present work is limited to the theoretical analysis of the numerical behaviour of the VGK code 

solution to the 2-D potential flow partial differential equations. VGK runs with changing numerical parameters 

(table 1 above) and flow conditions to simulate difficult flow conditions were not performed here as these will 

have to cover various free stream conditions (different free stream Mach numbers, M ,and angles of attack), 

and in the present author’s opinion these are to be left to a separate future work. 

3. Definitions 

In any numerical scheme three main features must be addressed in order to provide sound judgement on the 

behaviour of the scheme in terms of the overall accuracy and representation of the original or discretized partial 

differential equation. These features are the concepts of consistency, stability and convergence of the scheme 

[7], which are briefly defined as follows: 

3.1 Consistency 

A requirement that imposes a condition on the structure of the numerical formulation (discretized equation) and 

how well it represents the original partial differential equation. This provides an insight into the order of the 

discretization error. This is analysed by writing out the Taylor series expansion of the discretized elements and 

comparing the result to the original partial differential equation.  The discretization error will then be the 

difference between the two and the consistency requirement states that “The discretization error will tend to 

zero when the mesh spacing in the discretization direction (or a time step) tends to zero” [7]. 

3.2 Stability 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

21 
 

This condition establishes the relationship between the numerical and exact solutions of the discretized equation. 

Within the stability requirements it is emphasized that this difference (error) should diminish (or be bound) for 

n , where n represents the total number of iterations (or time steps). 

3.3 Convergence  

This represents the relationship between the numerical solution and the exact solution of the differential 

equation. The difference between these two solutions is known as the solution error and the convergence 

condition establishes that this error should reduce to zero as the grid is refined. It is usually very difficult to 

establish the convergence criterion theoretically particularly for flow cases where the exact solution is not 

known. The convergence of a finite difference scheme is related to the stability and consistency by the Lax 

equivalence theorem [9], which states that: “Given a well-posed linear initial value problem and a finite 

difference approximation to it that is consistent, stability is the necessary and sufficient condition for 

convergence” The Lax equivalence theorem is of great importance since it is relatively easy to show the 

stability of an algorithm and its consistency with the original differential equation. However as most real flow 

problems are non-linear and are boundary or mixed initial/boundary value problems the Lax theorem cannot 

always be applied hence it should be interpreted as providing the necessary but not always sufficient condition 

for convergence. 

4. Investigating the Consistency, Stability & Convergence of the VGK Scheme 

The investigations here are based on analysing the central and backward difference schemes of the Taylor series 

expansions of the full potential flow equation applied to the circular grid points (i. j).   

The potential flow equation used in VGK is given by (see definitions of mathematical symbols): 

 
 
   0

2

22)(

221

222

22222









r

r

rrr

fvrfuvur

vuar

vuvarvurua





                 (1) 

Where: a is the local speed of sound , vu,  the tangential and normal velocities components respectively in the 

computational plane, ,r   the grid radial and angular coordinates, the velocity potential in the computational 

plane. (With subscripts indicating first, second and mixed derivatives with respect to ( r ,  ), f  the modulus of 

the mapping function with subscripts in  ,r representing derivatives of  f with respect to  ,r  The code 

solves for the flow potential ( ) using central differencing scheme for the subsonic flow regions and 

forward/backward differencing scheme for the supersonic flow regions [10]. Central differencing, applied 

throughout the subsonic region along the tangential and radial lines: 











2

,1,1 jiji
 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

22 
 

r

jiji

r







2

1,1,
 

2

1,,1,

)(

2

r

jijiji

rr






 

2

1,,,1

)(

2









 jijiji
 











r

jijijiji

r
4

1,11,11,11,1
 

Where the (i, j) subscripts representing the  i
th

 and j
th  

steps and  ,,, rr  are the step sizes in the 

,r directions respectively in the computational plane: 

In the supersonic (hyperbolic) region backward differences are used on the upper surface )2(    and 

forward differences on the lower surface ( )0   in the following manner: 

)33()2()( ,3,2,1,,2,1,

2

jijijijijijiji      

and 

)22(

)(24

1,1,11,21,21,11,

1,1,11,11,





jijijijijiji

jijijijirr











 
 

where for the indices  ji, , the upper sign is used for the backward difference scheme and the lower sign for 

the forward scheme. 

The governing equation (1) is of the form: 

DCBA rrr                                                                                                                               (2) 

where: 

)( 22 uaA  , rvuB 2 , )( 222 varC                                                                                                                                      

    rfvfuvurvuarvuD rr  


2212222                                                               

(3) 

and applies outside the viscous layers. To solve for  a finite difference representation of (2) is formed by a tri-

diagonal system of equations: 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

23 
 

ijjiijjiijjiij dcba   1,,1,                                                                                                                   (4) 

Here jijijiji dcba ,,,, ,,,  are the coefficients of the potential equation the finite difference form.  

4.1 Analysing the Consistency of VGK 

it can be shown that for the centre-differencing scheme, the total truncation error is given by: 

(Error) total =         rrrrr

B
r

CA


222

61212
                                                               (5) 

Equation (5) is second order in   ,r , and shows that the central difference scheme is consistent. This, also, 

gives the coefficients of the fourth derivatives of  in  and r (or the third derivatives of the velocities  and 

r ) i.e. an odd-order derivative as the dominant term in the truncation error. This represents a diffusive term 

(dispersive) implicitly imbedded in the central difference scheme [7]. The exact amount and sign (whether 

lagging or leading) of this dispersion will, of course, be dependent on the actual values of the coefficients 

,A ,B and C at every point in the grid. It is, however, difficult to assess analytically the value of (5) but it will 

not be of major concern here as it will be very small and tending to zero as both  and 0r . 

In the backward difference scheme, the truncation error takes the following forms: 

(Error) total =         rrrrr r
CB

A 
2

12
1

2
1                                                        (6) 

with   110   , equation (6) becomes: 

(Error) total =     rrrrr r
CB

A 









22

122
                                                                             (7) 

 With ( ) being an artificial viscosity parameter. When  0 equation (6) is first order accurate in   but 

second order in  r , however with the addition of   this becomes second order accurate in 

both  r & as shown in (7). Equation (7) represents the coefficients of the second derivatives of the 

velocities  & r and thus an artificial viscosity term (dissipative) that is either implicitly present in the 

original scheme (i.e. when 0 ) or explicitly modified by the introduction of the parameter  . The 

introduction of this parameter reduces the value of the implicit dissipation by a factor of  (or by  1 ) 

so that not too much dissipation is introduced in the original scheme. It is therefore immediately concluded that 

reducing  (by increasing the number of mesh radial lines) in the hope of increasing the accuracy, may in fact 

cause an oscillatory behaviour (because of the reduced dissipation). A qualitative assessment of this can be 

arrived at by explicitly writing the coefficients A and B in (7) using equation (3), hence (7) becomes: 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

24 
 

(Error) total =      rrvuua   )(222
                                                                            (8) 

In hyperbolic regions, for which this equation holds, and if flow conditions are leading to the formation of a 

shock wave, it is expected that u will be of the same order as a . Hence the coefficient of  will be very 

small making (8) tend to a negative value unless r is negative (maximum r ). With the presence of a small 

factor  2   and these conditions, the dissipation term is reduced and oscillatory behaviour may ensue. 

 Two other important aspects of the dissipation behaviour can be revealed through (8): 

 As  and r increase towards their maximum values then it is expected that their second derivatives 

(  and r ) will become negative. This will help make (8) positive thus aiding a stable behaviour of the 

scheme when au  . This situation happens when the flow is completely supersonic, for example upstream of 

a shock wave. 

 The velocity components in the computational plane ( u and v ) are related to the potential velocities 

(  and r ) and the mapping modulus (f) by:  

f
u

r









    , v f

r





 

Substituting these into (8) gives: 

(Error) total =     rrrfrfa     )( 4222222
                                                                (9) 

 It is apparent that for small values of f equation (9) may become negative, this is most likely the case when 

  is positive corresponding to a reducing tangential velocity; that is after passing the maximum suction 

point on the upper surface. This is related to the aerofoil geometry since 
d

ds
f  . It is thus to be expected that 

aerofoils with a sharp change in profile, (a sudden or a sharp decrease in the first derivative) occurring after the 

maximum suction point may exhibit an oscillatory behaviour particularly at reducing speeds (e.g. downstream 

of a shock wave). In summary, the above merely states that in a hyperbolic region dissipation is reduced as 

0 , and oscillatory behaviour will be promoted by a sharp reduction (or a sudden change) in profile first 

derivative after the maximum suction point and in the presence of  pressure gradient downstream of a shock 

wave.  Figure (1) below shows the variation of the dissipation factor  2   with the artificial viscosity 

parameter  at various values of  where it is shown that to keep a satisfactory level of dissipation at a 

reduced  the value of  must also be reduced. 

 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

25 
 

 

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

0.05

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1

 2 
Decreasing 

(3,2.25,2,1.5,1) deg.





 

Figure 1: Dissipation factor variation with artificial viscosity at different values of   

4.2 Analysis of the Stability of the VGK 

The Von Neumann method [7,8]
 
for stability analysis was used to check the stability of the VGK method. The 

analysis for the central and backward difference schemes reveals that the central difference scheme is 

unconditionally stable, thus any oscillatory behaviour at totally subsonic speeds (non-existing supersonic 

regions) will be geometry-related, and caused by the application of the cubic spline in the conformal mapping 

procedure to regions of highly varying curvature.  

For the hyperbolic regions the stability criterion is given by the following formula: 

         0.02222

4

2

3

222

2

222

1  vauaKuvKvaKuaK                                        (10) 

where ( vu, ) are the local stream-wise and normal velocity components in the physical plane, and 1K to 

4K are positive constants that depend on the mesh spacing (  & r ) and the artificial viscosity parameter  , 

defined by the following: 

  















4

2

1

441
4




K , 

  













42

2

r
K , 

  

2

3

2

8

1














r
K




, 

    













224

1
4

r
K




 

Figures (2,3) show the variations of 1K to 4K with the total number of mesh radial lines for 0.0 and 

0.1 respectively, (in VGK r is fixed at 
30

1
). These show that for 0 , 1K dominates the other 

coefficients for 220m where m is the total number of grid radial lines. This corresponds 

to rad03.0 approximately. On the other hand for 1 , 1K starts to dominate the other coefficients for 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

26 
 

160m , that is rad039.0 . With  0 this domination is also more apparent when the 

coefficients 1K to 4K are normalized by 1K as shown in figure (4). Here the region of domination of  1K is 

shown (inset of figure 4). 

0.00E+00

2.00E+06

4.00E+06

6.00E+06

8.00E+06

1.00E+07

1.20E+07

1.40E+07

0 50 100 150 200 250 300 350 400

Total number of radial grid lines

K1

K2

K3

K4

m

 

Figure 2: Variations of the stability equation coefficients with number of grid lines at 0  

0.00E+00

2.00E+07

4.00E+07

6.00E+07

8.00E+07

1.00E+08

1.20E+08

1.40E+08

0 50 100 150 200 250 300 350 400

 Total number of radial grid lines

K1

K2

K3

K4

m

 

Figure 3: Variations of the stability equation coefficients with number of grid lines at 1  

-500

0

500

1000

1500

2000

2500

3000

3500

4000

4500

0 50 100 150 200 250 300 350 400

Total number of radial grid lines (m)

K2/K1

K3/K1

K4/K1

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

40 80 120 160 200 240 280 320 360

Region of K1 domination

 

Figure 4: Stability equation coefficients normalized by K1 )0(   



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

27 
 

 Examining figures (2-4) in conjunction with equation (16) reveals that when: 

 au   

 The term associated with 4K is negative, and together with the 2K  term dominate the 1K term when 

18080 m approximately, so in this case equation (16) is satisfied and the stability of the scheme is 

ensured. For 80m the 3K term dominates the other terms and unless this is offset by the value of uv , 

instabilities may occur. 

 u  aO and < a  

 Domination of the 2K term over the 3K  term will govern the stability condition, this is the case for 220m  

approximately, (see inset of figure 7.8), unless  22 vauv  . The above only forms a general guideline on 

the effects of the local flow conditions (when coupled with m ) on the overall stability behaviour of the 

hyperbolic scheme. The actual behaviour will of course depend on the particular values in equation (10). 

4.3 Analysing the Convergence of the VGK method 

The consistency and stability are necessary but not sufficient conditions for convergence. However, since VGK 

solves a tri-diagonal system (formed by equation 4), to satisfy convergence the diagonal dominance of the 

coefficient matrix in (4) must also be ensured [11]. The diagonal dominance of a matrix is ensured if the 

following condition is satisfied: 

iii bca                                                              (11) 

where ii ba , and ic  constitute the row elements of the tri-diagonal coefficient matrix. It can be shown that the 

values for ii ba , and ic in terms of the coefficients ,A ,B ,C and D of the difference equation (2) can be 

written as: 

 For the central difference (subsonic region)  

2

, Ca ji                                

 2

, 2 CAb ji   

2

, Cc ji   

 For the backward Difference (Supersonic region): 

 







 
 BCa ji 




4

22

,                                          



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

28 
 

  ACb ji   12 2

,       

 







 
 BCc ji 




4

22

,
                                          

where 
22 uaA  ,  rvuB 2 ,  

22 vaC   and 
r





 . 

Applying these to the dominance equation (11) the following conditions are arrived at: 

 For the central difference (subsonic) scheme 

0.042 4222   CACA                                                                                                   (12) 

 

 For the backward difference (supersonic) scheme 

 
 

 
 

2

2

22
4

2

1

2

2

1

1
241 B

A
C

A
CA

















                                                                (13) 

Equations (12 and 13) apply to all internal mesh points, but exclude the unit circle representing the boundary 

where in the tri-diagonal matrix the boundary is represented by the elements of the first row  1j . The 

central difference scheme (equation 12) shows unconditional diagonal dominance and since in this flow regime 

the scheme was shown to be unconditionally stable it is immediately concluded that convergence is assured 

provided no geometry-related instabilities take place. In the hyperbolic regions, however, equation (13) places 

certain restrictions on the mesh spacing  and r , the local speed conditions ( u , v ) and the amount of 

artificial viscosity  . Accordingly it is difficult to assess the conditions for which diagonal dominance can be 

ensured in these regions of the flow except to state that it appears from equation (19) that for all 1 (that is 

30

1
 r , or 188m approx.), increasing  (reducing m) helps convergence. Also since 

rvuB 2 and r decreases away from the boundary towards the centre of the mesh (i.e. towards physical 

infinity), then it can also be deduced from (13) that, given the same flow conditions diagonal dominance is more 

attainable away from the boundary. On the boundary  1r the coefficients in equations (12 and 13) take the 

following forms (see appendix B of reference 12). 

 Central difference (subsonic) scheme: 

01, ia  

 1

2

11, 2 CAbi                                                                                        (14) 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

29 
 

1

2

1, 2 Cci                                                                                                                                                               

 Backward difference (supersonic) scheme: 

01, ia                   

  1

2

11, 21 CAbi                                                                                                      (15) 

1

2

1, 2 Cci                                                                                                                                                                      

where the subscript 1 indicates conditions on the boundary. In the central difference scheme, diagonal 

dominance of the first row of the coefficient matrix is then given by: 

011

2
1

2  CAA   

which again gives unconditional diagonal dominance. 

 In the backward difference scheme the equivalent condition is: 

     0411 1

2

11  CAA   

this gives 

 

1

12

4

1

C

A



                                                                                                   (16) 

where the upper sign refers to positive 1A  (that is au  , but with 
222 avu  to stay supersonic) and the 

lower sign refers to negative 1A . This shows that the favourable increase in  for convergence as stipulated by 

equation (13) must be restricted by an upper limit given by equation (16).  This analysis depends on an assumed 

prior knowledge of 1A and 1C , which is difficult since these form part of the solution. Nevertheless, for each 

row of the diagonal matrix representing a new circle in the mesh an iterative procedure is used whereby each 

iteration cycle moves along the present circle calculating the potential values at each radial point. For each circle 

a number of iterations is carried out for the coarse and fine grids (default values are 100 and 200 respectively).  

A relaxation process is used to speed up and help convergence. This relaxation process is given by the following 

formula: 

  K

jiji

K

ji xxR ,,

1

, 1  
 

where K is an iteration number, jiR ,  is the initial calculated value of the potential obtained by solving the tri-

diagonal system at the point ( ji, ) and consequently forms the value of ji,  from the previous iteration. The  

relaxation factor ( x ) is defined by: 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

30 
 

}},{,{  ms xMaxxMinx                                                                                         (17) 

where, again, sx is an over-relaxation factor specified by the user,( 21  sx ,default value=1.9), mx  is an 

under-relaxation factor also specified by the user ( 10  mx ,default value=1), and  a variable which 

depends on local and free stream conditions and is defined by: 

 2

2
1

1
q

q

xx
x

crit

ms
s 












                                                                                                   (18) 

with q and critq being the total local and critical speeds respectively and q is given by: 

2222 avuq  . 

Thus the relaxation parameter is constant ( sx or mx ) unless   is used. Equation (17) determines the condition 

for which the use of  is invoked, and that is when: 

sm xx   

substituting for  from (18) gives: 

  s

crit

ms
sm xq

q

xx
xx 












 2

2
1

1
 

that is  

0
1

1
1

2

2















critq

q
 

or finally the relaxation factor x is used when  

122  qqcrit                                                                                        (19) 

which indicates that in the supersonic range given by (19) a variable relaxation factor is used and is dependent 

on the local speed conditions. The shaded area in figure (5) shows this. 

To investigate the effects of  it can be written from (18) as: 

2

2

1 LqL   

where from (18) 1L , 2L  are constants and are given by: 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

31 
 

21
1 crit

ms

q

xx
L




    ,    

2

2

2
1 crit

critsm

q

qxx
L




  

1<Xs<2

0<Xm<P

X=1

122  qqcrit
X=P

 

Figure 5: Region where a variable relaxation factor is used (shaded) 

Since q represents the gradient of the potential, , and forms part of the iterative solution, any oscillations 

imply that q itself will be oscillatory. This oscillatory behaviour may be further amplified through  via the 

2

1qL term, leading to an unstable (diverging) solution. To help convergence in such regions it would seem 

preferable to set 1mx so that  is never used, thus it is hoped, that a final non-oscillatory solution will be 

reached through increasing the number of iterations, if necessary. A final word of warning regarding 

convergence in VGK should be introduced here.  Iterations proceed either if the total number of iterations is 

below a pre-specified number or the tolerance is exceeded, whichever occurs first. If the specified maximum 

number of iterations is reached then VGK outputs the various flow parameters (surface and wake pressures, 

velocities...etc.) even if the tolerance is exceeded. Though these may not represent a final converged (or 

accurate) solution. To investigate the accuracy of the solution the convergence history against the number of 

iterations (which can be increased successively for a number of runs) must be checked. This can be done by 

inspecting the residuals (output in the file with extension “.IEH”) for each of the successive runs. The most 

reliable solution will then be the one after which any increase in the number of iterations does not produce a 

significant change in the residuals.  

5. Results and Conclusions 

The preceding analysis showed that the numerical calculations performed by the Viscous Garabedian & Korn 

(VGK) code are: 

A. Consistency: 

I. The central differencing scheme is unconditionally stable. (the truncation error is second order accurate 



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

32 
 

in (  ) and ( r ) 

II. The backward differencing scheme is first order accurate in (  ) and second order accurate in ( r ). 

An artificial viscosity parameter ( ) is added to make it second order in (  ). However, reducing 

(  ) may lead to oscillatory behaviour particularly in regions of the flow where strong shock wave(s) 

may occur and a large value of ( ) is used. Default values for ( ) are shown in table 1. 

III. In the backward differencing scheme, oscillatory behaviour may also occur due to the mapping 

procedure specifically if a large reduction in profile (aerofoil) first derivative occurs downstream of the 

maximum suction point (lowest pressure point on the upper surface of the aerofoil). 

B. Stability 

I. The central differencing scheme is unconditionally stable. However, any oscillatory behaviour in fully 

subsonic flow will only be due to the profile input data (if not a smooth function) and the resulting 

cubic spline used in the conformal mapping procedure. 

II. The backward differencing scheme is conditionally stable based on equation (10) and the situations 

described therein. 

C. Convergence 

I. The central difference scheme is unconditionally convergent. 

II. For the backward differencing scheme, relaxation factors are used to speed up the iterative process and 

help conversion. These are stipulated in equation (18, 19) and figure (5). 

6. Recommendations 

Based on the results of this study it is recommended that: 

A. The VGK method can be used and will provide reliable data (very close to experimental results) over 

profiles defined by smooth functions up to free stream speeds and angles of attack that will not cause 

the flow to reach a critical Mach number (fully subsonic flow). 

B. For flows where local weak shock waves (Mshock<1.4) VGK can still be used but a careful examination 

of the iterative procedure must be considered. (an output file of VGK with an (.IEH) extension 

facilitates this check)  

D. For difficult flows (strong shock waves) and/or non-smooth profiles (for example aerofoils with 

indentations), the grid resolution effects  (  , ( r ) on the consistency, stability and convergence of 

the VGK must be investigated further. The present author modified the code by increasing the grid 

resolution to (360) radial lines instead of the original maximum (160) lines. The present author intends 

to investigate the effects of this modification on the VGK performance. The results of the intended 

study will be in a separate future paper.  



 
American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2020) Volume 70, No  1, pp 16-33 

 

33 
 

References 

[1]. Freestone, M.M. "ESDU Transonic Aerodynamics Vol.3, 1997" 

[2]. Jameson,A. "Full potential,Euler and Navier-Stokes Schemes",chapter 3 in "Applied computational 

aerodynamics",progress in aeronautics and astronautics,Henne,P.A, Editor,Vol.125,1990 

[3]. Anderson, J.D.  "Computational fluid dynamics, the basics with applications". McGraw Hill Int., 1996. 

[4]. Hutchison M J. "Investigation into the ability of two CFD codes to model aerodynamic effects of 

airframe surface imperfections", ESDU Paper 365,ESDU. July 1994. 

[5]. Salah.J.S EL-IBRAHIM. “Prediction of the effects of aerofoil surface irregularities at high subsonic 

speeds using the viscous Garabedian & Korn method”. PhD thesis. Department of Civil Aeronautical 

and Mechanical Engineering, University of Hertfordshire, U.K., unpublished work,2000. 

[6]. Freestone, M.M. "VGK Method for two-dimensional aerofoil sections- part1: Principles and results". 

Draft Data Item, ESDU international, Oct.1996. 

[7]. Hirsch C. "Numerical computation of internal and external flows". Volume 1, "Fundamentals of 

numerical discretization". Wiley series in numerical methods in engineering ,John Wiley & Sons. 1989. 

[8]. Hirsch C. "Numerical computation of internal and external flows". Volume 2, "Computational methods 

for inviscid and viscous flows". Wiley series in numerical methods in engineering ,John Wiley & Sons. 

1990. 

[9]. Hoffman J.D. "Numerical methods for engineers and scientists". McGraw-Hill International 

editions.1993. 

[10]. Murman, E.M. Cole, J.D. "Calculation of plane steady flows". AIAA Journal.Vol.9, No.1. January 

1971. 

[11]. Fletcher C.A. "Computational techniques for fluid dynamics, Volume I, -Fundamentals and general 

techniques". Springer-Verlag ,2nd edition,1991. 

[12]. Collyer, M.R. "An extension to the method of Garabedian & Korn for the calculation of transonic flow 

past an aerofoil to include the effects of a boundary layer and wake". Reports &Memoranda 

No.3828.RAE, July 1977.  

[13]. Lock, R. C. , Collyer,M.R " A modification to the method of Garabedian & Korn". Numerical methods 

for the calculation of inviscid transonic flows. FREDR VIEWE GUND S. 1981