About Viiflow

Viiflow is an aerodynamic analysis method. It enables the fast evaluation of airfoil characteristics, such as lift and drag coefficients or transition locations. It is currently suitable for subsonic cases of single or multi-element airfoils without strong wake confluence or massive separation. It is especially useful for coupled problems, such as fluid-structure interaction or optimization. To cite viiflow, use [1].

Module Installation

Viiflow comes as a 64bit Cython wheel for Windows and Linux. Cython wheels are installation packages of compiled modules usable from Python. To use it, you therefore need Python (3.7+).

Now, if you are note comfortable already with Python you may want to start here. While no serious programming is necessary to use viiflow, it will be very beneficial if you know how statements and loops work to define your airfoil analysis.

If you have little experience with installing and running Python on Windows: I found it easiest to rely on scoop and use:

scoop install python

For the examples you will also need jupyter, so use pip install jupyter after installing python.

If you have little experience with Python on Linux: By default, a lot of Linux distributions see Python 2.7 as the default version, though this is changing. Check whether you have Python3.7 installed by using python3 in a shell. You may need to use python3 -m pip install jupyter, python3 -m jupyter notebook for all commands.

To install the wheels, use

pip install viiflow-[...].whl
pip install viiflowtools-[...].whl

from a shell.

While viiflow is compiled cython code, it calls numpys linear algebra solve routine at every iteration. For this to be efficient, you may need an efficient numpy installation. Here you can find wheel distributions of numpy for windows built with Intel MKL.

Module usage

The viiflow library has four methods:

  • setup: This function builds a setup structure that is used in the latter methods
  • init: This function uses the parameters to set up the data structures, calculate the panel operators, initialize the boundary layers and so on.
  • iter: This function uses the output of init to iterate the parameters until convergence.
  • set_forced_transition: Set current boundary layer to assume a forced transition


The setup routine can be called with a list of parameters and returns a parameter structure. The call

import viiflow as vf
s = vf.setup(Re=1e6,Alpha=10)
s.Ma = 0.1

defines the setup structure s with the parameters Reynolds number 1e6 and angle of attack 10°. After creation, the Mach number was changed from its default value (0) to 0.1. The following parameters can be used. The parameters not available in Basic are only available in the Pro version of viiflow.

Type Name Explanation Default Basic
double Re Reynolds Number 1E6
double Ma Mach Number (<1) 0
double Ncrit Critical amplification factor 9
int Itermax Maximum number of Newton iterations 100
bool IterateWakes Recalculate wake shape during iterations 0
int Substeps No of steps between two panel nodes (*1) 1
bool Silent Do not print info during iterations 0
bool Gradients Calculate gradients 0 ×
bool VirtualGradients Calculate gradients w.r.t. virtual displacement 0 ×
double Alpha Angle of attack 0
double PitchRate Pitching rate about (0,0) (*2) 0
double LocusA G-beta constant A (*3) 6.75
double LocusB G-beta constant B (*3) 0.83
double ViscPwrExp Viscosity-temperature model exponent 0.7
double WakeLength Length of all wakes w.r.t. chord. 1
double Tolerance Convergence criterion. 1e-4
int IncompressibleBL If>0, calculate boundary layer with Ma=0 0
double StepsizeLimit Additional limit on Newton step, between 0 and 1 1
double ShearLagLambdaWake Parameter \lambda in the lag equation (Wake) .9
double ShearLagLambdaFoil Parameter \lambda in the lag equation (Airfoil) 1
int HalfWakes If>0, use two distinct boundary layers for wake 0
int ShearLagType Defines choice of K_C function in shear lag equation (*4) 0
double TransitionStepLimit Restrict movement of transition inf

Substeps (*1)

During one iteration viiflow marches from the stagnation point along the pressure and suction side of the airfoil to solve the boundary layer equations. The boundary layer equations are essentially a set of ordinary differential equations (ODEs).

By default, the discretization of the ODE on the surface of the airfoil is the same as the panel nodes, i.e. a step along the surface is a single step from panel node to panel node.

By setting substeps to n, the path between to panel nodes is divided into n segments for the ODE solver. This does increase the numerical effort during the march, but not the effort for the much larger problem of the global Newton step.

Pitch rate (*2)

Pitching airfoils experience different flow conditions than non-pitching airfoils. A pitch rate is given as the physical pitch rate in rad per sec divided by the speed of the free flow. For example, setting pitch_rate to 40*pi/180/20 would be a pitch rate of 40°/s at a current speed of 20m/s. The pitching motion is assumed to be about (0,0), so arrange your geometry accordingly.

Usually the airfoils analyzed will have a chord length of 1m and the results are scaled to the application. To model a pitching airfoil (40°/s) for a model plane at 20 m/s that has a chord length of 7cm where viiflow is calculating the airfoil characteristics at a length of 1m, one would set the pitch rate to 40pi/180/200.07.

The quantity describes as well the chord to radius ratio for rotating airfoils, so set pitch_rate to c/R if given.

G-beta Locus (*3)

The parameters A and B modify the equations used in the turbulent boundary layer. They are part of an empirical correlation called the G-beta locus [2,5] and influence turbulent separation.

Exemplary values from [3] are : XFOIL default [2] (A=6.7, B=.75), Boeing (A=6.935, B=.70, reference not found by the author), Green et al.[6] (A=6.43, B=.80), RFOIL [3] and viiflow defaults (A=6.75, B=.83).

Shear-Lag type (*4)

Depending on this integer value, the factor K_c in the shear-lag equation is calculated differently. For a value of 0 the factor is a constant 5.6. This is the default case. For a value of 1, the RFOIL K_C function is used[3], that is K_C = 4.65-.95*tanh(.275*H-3.5). The shear stress follows the equilibrium value more slowly for large H. For a value of 2, the K_C function is calculated as described by Ye[10] following Thomas[9], that is K_C = 3.75*(3*H)/(H+2). This is similar to the behavior in XFOIL, for high H the shear stress follows the equilibrium faster. This has a significant effect near and post stall.


The init routine can be called to initialize the variables that are used during the iterations. These are

  • the panel structure p, which contains the airfoil geometries, wake geometries, viscid and inviscid solutions and the lift and moment coefficient. E.g. p.CL returns the calculated lift coefficient and p.foils[0].X returns the airfoil geometry as a 2xN ndarray.
  • the boundary layer structure bl, which is a list of structures for every airfoil. Every structure, among other things, contains the substructures bl_fl, the airfoil surface boundary layer, and bl_wk, the wake boundary layer.
  • The array x0, which is the variable used in the Newton iteration in iter.
import viiflowtools.vf_tools as vft

# Read Airfoil coordinates into numpy 2x220 array using a function from viiflowtools
RAE = vft.repanel(vft.read_selig("RAE2822.dat"),220)

# Init takes a list of airfoil geometries, here this list is a single airfoil.
# A single 2xN array is fine as well.
[p,bl,x] = vf.init([RAE],s) # [p,bl,x] = vf.init(RAE,s) works, too 


The above code read an airfoil into a numpy array, and lets the init function initialize p, bl and x. the ' print statement display first the vector of the inviscid solution (the edge velocity) and secondly the initial momentum thickness.


This routine is called to drive the problem towards a solution using s.itermax Newton iterations. A very simple call would be

for AOA in range(0,10):
    s.Alpha = AOA
    [x,flag] = vf.iter(x,bl,p,s)
    if flag: # if flag = 1: converged
       print('AOA %u CL %f CD %f'%(AOA,p.CL,bl[0].CD))

Above, the call to iter lets it run until convergence or the defined maxiumum number of iterations, overwriting p,bl and x in the process. If the iterations were successful, the lift and drag coefficients are printed.

A more advanced call is

res = None
grad = None
s.Itermax = 0 #No internal iteration
for iter in range(100):
    [x,flag,res,grad] = vf.iter(x,bl,p,s,res,grad)
    x -= 0.01*np.linalg.solve(grad,res)
    if np.sqrt(np.dot(res.T,res))<1e-5
        print('Now close enough for me!)

This allows for a fine-grained control of the iterations. Here, the iterations are stopped when the residual falls below 1e-5. For the usage of gradient information I suggest looking into the Fluid-Structure Interaction example.

grad contains

  • res_vd: gradient of residual w.r.t. virtual displacement
  • gam_vd: gradient of surface speed w.r.t. virtual displacement
  • cl_vd: gradient of lift coefficient w.r.t. virtual displacement
  • delta_vd: gradient of boundary layer displacement thickness w.r.t. virtual displacement
  • cd_vd: gradient of drag coefficient w.r.t. virtual displacement
  • gam_x = gradient of surface speed w.r.t. x
  • cl_x = gradient of lift coefficient w.r.t. x
  • cd_x = gradient of drag coefficient w.r.t. x.


This function calculates the velocity at points X given as numpy 2xN arrays. The calculation is based on the inviscid (panel) method, but includes the displacement body due to the boundary layer. Therefore, it is suitable to estimate the velocities sufficiently far away from the boundary layer. The function needs a memory buffer where the velocity is written to. The function needs 2xN buffers and positions and for single point use the 2D vectors need to be given to the function as shown in the example.

# Assuming we have a solution (p,bl) from viiflow we generate a single streamline
N = 1000
x = np.zeros((2,N))
x[0,0] = 0.1
x[1,0] = 0.1
v = np.zeros((2,1))
for k in range(1,N):
    # Normalize for streamline integration
    # Tiny step forward
    x[:,k] = x[:,k-1] + v[:,0]*1e-3

# Plot with airfoil
fig,ax = plt.subplots(1,1)


This function estimates the velocity profile of a single boundary layer node. For laminar boundary layer nodes, a Pohlhausen function is assumed. For turbulent surface boundary layer nodes the explicit Musker function and the Finley wake function as described in [7] is used. For turbulent wake nodes with low shape factors, a simple power law is assumed. The function takes to 1D buffers for the Y coordinate and the velocity U.

# Assuming we have a solution (p,bl) from viiflow we generate a single streamline
# Buffers
Y = np.zeros((50))
U = np.zeros((50))
# Calculate profile at trailing edge, pressure side
node = bl[0].bl_fl.nodes[-1]

# Plot profile
fig,ax = plt.subplots(1,1)

# Plot function in wake
# Plot profile


This function uses the function viscid_profile to calculate all velocity profiles in a boundary layer and calculates the normal flow velocities from the (incompressible) equation d/dx u + d/dy v = 0. The returned matrices Y,U and V are centered on the panels, between the nodes.

# Assuming we have a solution (p,bl) 
N = 50 # Num. of normal elements
[Y,U,V] = vf.viscous_flowfield(bl[0],N)


This routine is called to set a boundary layer to force transition at a given foil-coordinate x location

[p,bl,x] = vf.init([RAE],s)
trans_upper = 0.025 # Transition location on suction side
trans_lower = 0.4 # Transition location on pressure side
for AOA in range(0,10):
    s.Alpha = AOA
    if AOA>5: # Only use forced transition if AOA>5 for some reason
    [x,flag] = vf.iter(x,bl,p,s)
    if flag: # if flag = 1: converged
       print('AOA %u CL %f CD %f'%(AOA,p.CL,bl[0].CD))

The lists at the third and fourth argument can have as many entries as there are airfoils, or can be empty. If the given transition location is empty or some coordinate not within the airfoil coordinate range, no forced transition occurs. If only forced transition is allowed, set in addition s.ncrit = np.inf.


There exists a Basic version and a Pro version of viiflow. The Basic version does not allow

  • Multiple airfoil geometries.
  • Gradient calculations.



  • Added ShearLagType parameter to switch between different K_C functions in the shear-lag equation. The new default behavior is a constant value, to switch to the previous default (RFOIL) set ShearLagType to 1.
  • Added TransitionStepLimit parameter to restrict movement of transition point within iterations. This may improve convergence for certain scenarios, but do not use by default as it usually does not.
  • Renamed all setup parameters to use CamelCase, so you need to adapt your scripts. The confusing parameter equal_wakes is now called WakeLength.
  • Default Tolerance is now 1e-4, default wake length is now 1


  • Added functions to aid visualization
    • inviscid_velocity: Calculates the velocity at arbitrary points away from the boundary layer
    • viscid_profile: Calculates the local velocity profile of the boundary layer
    • viscous_flowfield: Used viscid_profile to calculate the complete boundary layer velocities
  • Added two-wake method, using two individual boundary layers for the top and bottom wake half (still work in progress)
  • Fixed error with (very) large blunt trailing edges, resulting in very high lift coefficients


  • Speed and stability improvements
  • Added weight parameter for the lag equation, set for wake and airfoil separately
  • If there are virtual displacements at the TE, iterate_wakes=True will move the wake with the displacement


  • Speed and stability improvements
    • Added setup parameter stepsize_limit
    • Increased internal tolerance of boundary layer solve to 1e-6
  • Fixed wrong assertions in airfoil check from 1.2.2


  • Check airfoil geometry on init


  • Added incompressible_bl parameter
  • Added inviscid CL
  • Added set_forced_transition function
  • Added tolerance as parameter


  • Adapted pressure calculation for pitching airfoils [8]
  • Added equal_wakes parameter


[1] Ranneberg, Maximilian. Viiflow—A New Inverse Viscous-Inviscid Interaction Method. AIAA Journal 57.6 (2019): 2248-2253.

[2] Drela, Mark, and Michael B. Giles. Viscous-inviscid analysis of transonic and low Reynolds number airfoils. AIAA journal 25.10 (1987): 1347-1355.

[3] Van Rooij, R. P. J. O. M. Modification of the boundary layer calculation in RFOIL for improved airfoil stall prediction. Report IW-96087R TU-Delft, the Netherlands (1996).

[4] Drela, Mark. Integral boundary layer formulation for blunt trailing edges. 7th Applied Aerodynamics Conference. 1989.

[5] Francis H. Clauser. Turbulent boundary layers in adverse pressure gradients. Journal of the Aeronautical Sciences 21.2 (1954): 91-108.

[6] Green, J. E., D. J. Weeks, and J. W. F. Brooman. Prediction of turbulent boundary layers and wakes in compressible flow by a lag-entrainment method. ARC R&M 3791 (1973).

[7] Rona, Aldo, Manuele Monti, and Christophe Airiau. On the generation of the mean velocity profile for turbulent boundary layers with pressure gradient under equilibrium conditions. The Aeronautical Journal 116.1180 (2012): 569-598.

[8] van der Horst, Sander, et al. Flow Curvature Effects for VAWT: a Review of Virtual Airfoil Transformations and Implementation in XFOIL. 34th Wind Energy Symposium. 2016.

[9] Ye, Boyi. The Modeling of Laminar-to-turbulent Transition for Unsteady Integral Boundary Layer Equations with High-order Discontinuous Galerkin Method. Thesis TU-Delft (2015).

[10] Thomas, J. Integral boundary-layer models for turbulent separated flows. 17th Fluid Dynamics, Plasma Dynamics, and Lasers Conference. 1984.