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 suitable for subsonic cases of single or multi-element airfoils without strong wake confluence or massive separation and it is especially useful for coupled problems, such as fluid-structure interaction or optimization. In addtion to free flow simulations, it has the ability to predict flow in pitching, cavitating, free-surface and ground-effect cases. To cite viiflow, use [1].
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 need Python 3.12.x, numpy 2.4.x and scipy. The examples use jupyter as well as sometimes other packages depending on the notebook, but they are not necessary to use viiflow.
If you are not 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.
To install the wheels, use
pip install viiflow-[...].whl
pip install viiflowtools-[...].whl
from a shell (or the appropriate peotry/uv commands).
Consider installing numpy and scipy wheels from here: https://github.com/urob/numpy-mkl
Important:
CmRefX=.25 in the setupThe standard viiflow library has four methods:
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 |
|---|---|---|---|
| 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 |
| 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 |
| bool | HalfWakes | 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 (*5) | inf |
| bool | IsCascade | If true, the analysis assumes a cascade of airfoils(*6) | 0 |
| double | CascadeStaggerHeight | Stagger distance of airfoils (*6) | 1 |
| double | CascadeStaggerAngle | Stagger angle of airfoils (*6) | 0 |
| bool | BypassTransition | Use Bypass transition [14] | 1 |
| bool | BypassA | Parameter A [14] | .1 |
| bool | BypassB | Parameter B [14] | .3 |
| int | FreezeWake | If true, the wake stays as is when changing aoa (*7) | 0 |
| int | Mirror | 0: None. 1: ground effect (y=0). -1: free surface (*9) | 0 |
| double | IntermittencyGrowthRate | 10: Roughly transition within 1/10/c | 10 |
| double | CLTarget | If not inf, look for aoa to fit desired lift (*10) | np.inf |
| int | FreeSurface | If true, the | 0 |
| double | FroudeNumber | Froude Number for Free Surface calculations (*11) | 1.0 |
| double | FreeSurfMultLeft | Free surface length for visualization upwind | 4.0 |
| double | FreeSurfMultRight | Free surface length for visualization upwind | 6.0 |
| double | CmRefX | Reference position for moment coefficient | 0.0 |
| double | CmRefY | Reference position for moment coefficient | 0.0 |
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.
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.
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 but noted in [3]), Green et al.[6] (A=6.43, B=.80), RFOIL [3] and viiflow defaults (A=6.75, B=.83).
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.
This restricts the movement of the transition position within one Newton step. This does not change the converged solution, but may help convergence in some cases. Set to the (absolute) length the position along the surface may change.
The cascade flow model implementation used in the calculation mostly follows [11], which is similar to the original formulation of Hess and Smith [12]. The parameters are best explained with a graph, and therefore please take a look at the example in the notebooks section. The cascade is an infinite stack of airfoils, aligned on a common vertical axis.
This feature seems unhelpful. But e.g. if you are looking for the AOA for a specific lift and are not using the internal aoa search routine, you might want to use such a feature.
At a ray through zero in the direction of the free flow the geometry is mirrored. With Mirror=1 a ground effect is modelled and the flow is parallel to the ground, which is also the high Froude number limit of a free surface. With Mirror=-1 instead of parallel flow we enforce a pressure coefficient of 1 at y=0, which is the low Froude number limit of a free surface. Mirror=0 disables it.
This treats the angle of attack as an additional unknown inside the call to "iter", and x is increased by one item.
However, res and grad are not increased, the still are the residual and gradient assuming constant aoa.
Internally iter solves the problem with the additional gradients with respect to aoa, which are available in the gradients output.
Following [15], though with a closed form solution of the integral, a free surface computation is conducted at the ray through 0 in the direction of the free flow. THe Froude number is a setup parameter. The velocity at the free surface can be inspected in p.SurfaceVel which are meant at p.SurfaceNodes. The position of the nodes are defined using FreeSurfMultLeft and Right, the length is defined by these multipliers times the total geometry span in x.
The init routine can be called to initialize the variables that are used during the iterations. These are
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.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.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
print(p.gamma_inviscid)
print(bl.bl_fl.nodes.theta)
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 maximum 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!)
break
This allows for a fine-grained control of the iterations. Here, the iterations are stopped when the residual falls below 1e-5 and we take only 0.01 times the Newton direction as a step.
For the usage of gradient information I suggest looking into the Fluid-Structure Interaction example.
The methods initand iter work with arrays of boundary layer objects bl.
For every airfoil there is a boundary layer object one can access by index, i.e. the boundary layer for the first surface is bl[0].
Each object has the following structure.
bl[0].S # Airfoil surface coordinate, CCW from trailing edge
bl[0].Sw # A wake coordinate, increasing along distance
bl[0].Re # Reynolds number
bl[0].Ma # Mach number
bl[0].ncrit # Critical amplification factor
bl[0].hte # Trailing edge thickness
bl[0].ST # Stagnation point along surface coordinate
bl[0].sti # Stagnation occurs between sti-1 and sti
bl[0].N # Number of surface elements
bl[0].NW # Number of wake elements
bl[0].substeps # Number of substeps, usually 1
bl[0].CD # Drag coefficient
bl[0].S_trans_forced_up # Surface coordinate at which transition is enforced
bl[0].S_trans_forced_lo # Surface coordinate at which transition is enforced
bl[0].bl_fl.nodes # Boundary layer nodes on the surface
bl[0].bl_fl.node_tr_up # Extra node at transition
bl[0].bl_fl.node_tr_lo # Extra node at transition
bl[0].bl_wk.nodes # Boundary layer nodes on the wake
All boundary layer variables for surfaces and wakes are in numpy recarray substructs in that object. Each element of these recarrays I call node. Each node has the following structure, following the nomenclature of XFOIL [2,13].
node = bl[0].bl_fl.node_tr_up # Look into transition node
node.typ # 1: Laminar. 2: turbulent. 3: wake
node.delta # Displacement thickness
node.theta # Momentum thickness
node.ue # Edge velocity
node.nct # Laminar: Amplification factor. Turbulent/Wake: Shear stress coefficient
node.xi # Surface coordinate, starting from 0 at stagnation on suction and pressure side
node.UE # Mach corrected edge velocity (Eq 35 in [13])
node.H # Shape factor delta/theta
node.HK # Kinematic shape parameter
node.HS # Kinetic energy shape parameter, H* in [2]
node.HC # Density shape parameter, H** in [2]
node.RT # Momentum thickness Reynolds number
node.MSQ # Local mach number
node.CF # Skin-friction coefficient
node.DI # Dissipation coefficient
node.NX # Local rise in amplification factor, d/dxi n
node.DE # Boundary layer thickness (Eq 23 in [2])
node.US # Slip velocity (Eq 21 in [2])
node.CQ # Equilibrium shear stress coefficient (Eq 24 in [2])
node.xDt # xi/theta
node.delta_wake # Additional displacement thickness on wake near blunt trailing edge
# Plot a quantity
ax.plot(bl[0].bl_fl.nodes.xi,bl[0].bl_fl.nodes.H)
The methods initand iter work with the object p containing everything related to the inviscid part of the solver.
It contains
p.NS # Number of surface coordinates
p.NW # Number of wake coordinates
p.alpha # Angle of attack (in deg)
p.pitch_rate # Pitch rate divided by airspeed (or flow curvature chord/flow radius c/R)
p.IsCascade # Flag, indicating if a cascade is calculated
p.CascadeStaggerHeight # Quantity describing cascade geometry. See cascade example.
p.CascadeStaggerAngle # Quantity describing cascade geometry. See cascade example.
p.mach # Mach number
p.CL # Lift coefficient
p.CLi # Inviscid lift coefficient (no boundary layer)
p.CM # Moment coefficient (around 0 with respect to the airfoil coordinates).
# Note: If you have only one airfoil that has a length of 1, you may want to shift the X coordinate of your airfoil by -.25 to get the usual moment coefficient around .25c.
p.CascadeDeflectionAngle # Result of cascade calculation, deflection angle of flow downstream
p.CascadeDownstreamVelocity # Result of cascade calculation, downstream velocity
p.CascadeDeflectionAngleInviscid # Same as above, but without boundary layer
p.CascadeDownstreamVelocityInviscid # Same as above, but without boundary layer
p.gamma_inviscid_0 # Signed surface velocity at AOA = 0°, no boundary layer
p.gamma_inviscid_90 # At AOA = 90°
p.gamma_inviscid_pitch # At AOA = 0°, pitch_rate = 1
p.gamma_inviscid # # Signed surface velocity at current AOA, no boundary layer.
p.cp # Pressure coefficient
p.cp_inviscid # Pressure coefficient, no boundary layer
p.gamma_viscid # Signed surface velocity including boundary layer
p.max_thickness # For multi-element configurations, the boundary layer or virtual displacement should not cross into another geometry. For every surface coordinate, this is the maximum of the local displacement.
p.min_thickness # See above, but in the other direction (inwards).
In addition, the airfoil and wake coordinates and some helpful quantities are also store in p. I will assume in the following the first airfoil is to be accessed. Otherwise, use another index.
p.foils[0].X # Airfoil coordinates, 2xN
p.foils[0].S # Airfoil surface coordinate, counting CCW from TE
p.foils[0].normals # Airfoil surface normals, 2xN
p.foils[0].N # Number of coordinates
p.foils[0].sharp # Flag indicating blunt or sharp TE
p.foils[0].VD # Virtual displacements used
p.foils[0].delta # Boundary layer displacement used
p.foils[0].hte # Thickness of blunt TE
p.wakes[0].X # Wake surface coordinates, 2xN
p.wakes[0].S # Wake surface coordinate, counting CCW from TE
p.wakes[0].normals # Wake surface normals, 2xN
p.wakes[0].N # Number of coordinates
p.wakes[0].delta # Boundary layer displacement used
The gradients are split into partial and total gradients.
Viiflow solves the inviscid equations together with the viscous boundary layer equations by using a Newton method to set res to 0 (or small enough) using the gradient of the residual with respect to the variables x.
At each iteration, a step is given by dx=inv(grad)*res.
At convergence, res should be sufficiently close to 0.
At every iteration, we can also get the gradient of the residual or any other quantity with respect to other variables, such as the angle of attack or the virtual displacements.
These gradients are within the partial structure.
These are used e.g. in the Fluid-Structure Interaction notebook, where a common solver is used to drive the viiflow res and a structural equation to 0 at the same time.
However, usually one lets viiflow converge to a solution of the fluid dynamics equations and is only interested in the outcome of a converged solution. In this case, the total gradients need to be used, which return the sensitivity of a converged solution.
Let NF be the number of airfoils, N be the number of airfoil surface coordinates, NVD the amount of virtual displacements and NX the number of viiflow variables.
Then these are the available gradients and their dimensions. _vd always refers to the gradient with respect to the virtual displacements, _aoa with respect to the angle of attack (in rad) and x with respect to the viiflow variables.
| Name | Dimensions | Explanation (w/o repetitions) |
|---|---|---|
| partial.res_vd | (NX,NVD) | Equations residual |
| partial.res_aoa | (NX,1) | Assuming AOA in rad |
| partial.gam_x | (N,NX) | Signed surface velocity |
| partial.gam_aoa | (N,1) | |
| partial.gam_vd | (N,NVD) | |
| partial.cl_x | (1,NX) | Total lift |
| partial.cd_x | (NF,NX) | Drag for every airfoil |
| partial.cm_x | (1,NX) | Total moment |
| partial.cm_aoa | (1,1) | |
| partial.cl_vd | (1,NVD) | |
| partial.cm_vd | (1,NVD) | |
| total.gam_aoa | (N,1) | |
| total.cl_aoa | (1,1) | |
| total.cm_aoa | (1,1) | |
| total.cd_aoa | (NF,1) | |
| total.delta_vd | (NX,NVD) | Boundary-layer displacement thickness |
| total.gam_vd | (NX,NVD) | |
| total.cl_vd | (1,NVD) | |
| total.cm_vd | (1,NVD) | |
| total.cd_vd | (NF,NVD) |
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):
vf.inviscid_velocity(p,x[:,k-1:k],v[:,0:1])
# Normalize for streamline integration
v/=np.sqrt(v[0]**2+v[1]**2)
# Tiny step forward
x[:,k] = x[:,k-1] + v[:,0]*1e-3
# Plot with airfoil
fig,ax = plt.subplots(1,1)
ax.plot(x[0,:],x[1,:])
ax.plot(p.foils[0].X[0,:],p.foils[0].X[1,:])
ax.axis('equal')
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]
vf.viscid_profile(node,Y,U)
# Plot profile
fig,ax = plt.subplots(1,1)
ax.plot(U,Y)
# Plot function in wake
vf.viscid_profile(bl[0].bl_wk.nodes[3],Y,U)
# Plot profile
ax.plot(U,Y)
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
vf.set_forced_transition(bl,p,[trans_upper],[trans_lower])
[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.
The cavitation solver has another set of methods:
problem = vf.cav_init(AF,s): Initialize a cavitation problem using the Airfoil(s) AF and the setup sflag = vf.cav_solve(problem,problem.s,sig,nsubiter,callback,tol,restart,cavWeight): Solve a cavitation problem with Sigma sig.You cannot combine this at this point with a target-lift-coefficient, the angle of attack stays fixed.
sigthat was used to solve the problemCavResidual,Step,NormalResidual which can be used to observe what the solver does along the iterations.problem as an input. See example notebook. Use None to ignore.True (default), the solver starts by solving the non-cavitating case and then goes on from there. Otherwise it starts from whatever currently resides in problem. False is experimental.cp>cp_cav against cavity only where cp=cp_cav. If 1 does not result in a converged case because the already found separation after a cavity end does not want to "wander" to the true cavity end, another factor may. But only use values between .1 and 20 or so. See example notebook.There exists a Basic version and a Pro version of viiflow. The Basic version does not allow
[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.
[11] Mokry, M. The Complex Variable Boundary Element Method for potential flow problems. Boundary Elements and Other Mesh Reduction Methods Twenty-eight 42 (2006): 211.
[12] Hess, J. L., & Smith, A. O.. Calculation of potential flow about arbitrary bodies. Progress in Aerospace Sciences, 8 (1967), 1-138.
[13] Drela, Mark. XFOIL: An analysis and design system for low Reynolds number airfoils. Low Reynolds number aerodynamics. Springer, Berlin, Heidelberg, 1989. 1-12.
[14] Drela, Mark. MISES implementation of modified Abu-Ghannam/Shaw transition criterion. Mises User’s Guide, MIT (1995).
[15] Chen, Zhi-min, and W. G. Price. Dissipative free-surface solver for potential flow around hydrofoil distributed with doublets. Applied Mathematics and Mechanics 33 (2012): 1467-1480.