# RCAIDE/Library/Methods/Aerodynamics/Vortex_Lattice_Method/VLM.py
#
# Created: Aug 2025, M. Clarke
# Modified:Apr 2026, S. Shekar, A. Molloy, M. Clarke
# May 2026 M. Clarke
# ----------------------------------------------------------------------
# Imports
# ----------------------------------------------------------------------
# package imports
from RCAIDE.Framework.Core import Data
from .compute_wing_induced_velocity import compute_wing_induced_velocity
from .generate_vortex_distribution import generate_vortex_distribution
from .compute_RHS_matrix import compute_RHS_matrix
from scipy.integrate import trapezoid
import numpy as np
# ----------------------------------------------------------------------
# Vortex Lattice
# ----------------------------------------------------------------------
[docs]
def VLM(conditions,settings,geometry):
"""Uses the vortex lattice method to compute the lift, induced drag and moment coefficients.
The user should be forwarned that this will cause very slight differences in results for 0 deflection due to
the slightly different discretization.
The user has the option to use the boundary conditions and induced velocities from either RCAIDE
or VORLAX. See build_RHS in compute_RHS_matrix.py for more details.
By default in Vortex_Lattice, VLM performs calculations based on panel coordinates with float32 precision.
The user may also choose to use float16 or float64, but be warned that the latter can be memory intensive.
The user should note that fully capitalized variables correspond to a VORLAX variable of the same name
Assumptions:
The user provides either global discretezation (number_spanwise/chordwise_vortices) or
separate discretization (wing/fuselage_spanwise/chordwise_vortices) in settings, not both.
The set of settings not being used should be set to None.
The VLM requires that the user provide a non-zero velocity that matches mach number. For
surrogate training cases at mach 0, VLM uses a velocity of 1e-6 m/s
Source:
1. Miranda, Luis R., Robert D. Elliot, and William M. Baker. "A generalized vortex
lattice method for subsonic and supersonic flow applications." (1977). (NASA CR)
2. VORLAX Source Code
Inputs:
geometry.
reference_area [m^2]
wing.
spans.projected [m]
chords.root [m]
chords.tip [m]
sweeps.quarter_chord [radians]
taper [Unitless]
twists.root [radians]
twists.tip [radians]
xz_plane_symmetric [Boolean]
aspect_ratio [Unitless]
areas.reference [m^2]
vertical [Boolean]
origin [m]
fuselage.
origin [m]
width [m]
heights.maximum [m]
lengths.nose [m]
lengths.tail [m]
lengths.total [m]
lengths.cabin [m]
fineness.nose [Unitless]
fineness.tail [Unitless]
settings.number_of_spanwise_vortices [Unitless]
settings.number_of_chordwise_vortices [Unitless]
settings.use_surrogate [Unitless]
settings.propeller_wake_model [Unitless]
settings.use_VORLAX_matrix_calculation [boolean]
settings.floating_point_precision [float16/32/64]
conditions.aerodynamics.angles.alpha [radians]
conditions.aerodynamics.angles.beta [radians]
conditions.freestream.mach_number [Unitless]
conditions.freestream.velocity [m/s]
conditions.static_stability.pitch_rate [radians/s]
conditions.static_stability.roll_rate [radians/s]
conditions.static_stability.yaw_rate [radians/s]
Outputs:
results.
CL [Unitless], CLTOT in VORLAX
CDi [Unitless], CDTOT in VORLAX
CM [Unitless], CMTOT in VORLAX
CY [Unitless], Total y force coeff
CRTOT [Unitless], Rolling moment coeff (unscaled)
CL_mom [Unitless], Rolling moment coeff (scaled by b_ref)
CNTOT [Unitless], Yawing moment coeff (unscaled)
CN [Unitless], Yawing moment coeff (scaled by b_ref)
CL_wing [Unitless], CL of each wing
CDi_wing [Unitless], CDi of each wing
cl_y [Unitless], CL of each strip
cdi_y [Unitless], CDi of each strip
alpha_i [radians] , Induced angle of each strip in each wing (array of numpy arrays)
CP [Unitless], Pressure coefficient of each panel
gamma [Unitless], Vortex strengths of each panel
Properties Used:
N/A
"""
S_ref = geometry.reference_area
c_ref = geometry.reference_chord
b_ref = geometry.reference_span
x_m = geometry.mass_properties.center_of_gravity[0][0]
z_m = geometry.mass_properties.center_of_gravity[0][2]
# ---------------------------------------------------------------------------------------
# Generate Panelization and Vortex Distribution
# ------------------ --------------------------------------------------------------------
VD = generate_vortex_distribution(conditions,settings,geometry)
settings.vortex_distribution.chord_lengths = VD.chord_lengths[VD.leading_edge_indices].reshape(len(VD.n_sw),np.sum(VD.n_sw[0]))
settings.vortex_distribution.n_sw = VD.n_sw
settings.vortex_distribution.n_cw = VD.n_cw
settings.vortex_distribution.n_w = VD.n_w
settings.vortex_distribution.chord_widths = VD.chord_widths
settings.vortex_distribution.leading_edge_sweeps = VD.leading_edge_sweeps
settings.vortex_distribution.XA1 = VD.XA1
settings.vortex_distribution.XA2 = VD.XA2
settings.vortex_distribution.XB1 = VD.XB1
settings.vortex_distribution.XB2 = VD.XB2
settings.vortex_distribution.YA1 = VD.YA1
settings.vortex_distribution.YA2 = VD.YA2
settings.vortex_distribution.YB1 = VD.YB1
settings.vortex_distribution.YB2 = VD.YB2
settings.vortex_distribution.ZA1 = VD.ZA1
settings.vortex_distribution.ZA2 = VD.ZA2
settings.vortex_distribution.ZB1 = VD.ZB1
settings.vortex_distribution.ZB2 = VD.ZB2
settings.vortex_distribution.X = VD.X
settings.vortex_distribution.Y = VD.Y
settings.vortex_distribution.Z = VD.Z
settings.vortex_distribution.XC = VD.XC
settings.vortex_distribution.YC = VD.YC
settings.vortex_distribution.ZC = VD.ZC
settings.vortex_distribution.Y_SW = VD.Y_SW
# unpack conditions--------------------------------------------------------------
pwm = settings.propeller_wake_model
K_SPC = settings.leading_edge_suction_multiplier
aoa = conditions.aerodynamics.angles.alpha
mach = conditions.freestream.mach_number
len_mach = len(mach)
#For angular values, VORLAX uses degrees by default to radians via DTR (degrees to rads).
#RCAIDE uses radians and its Units system. All algular variables will be in radians or var*Units.degrees
PSI = conditions.aerodynamics.angles.beta
PITCHQ = conditions.static_stability.pitch_rate
ROLLQ = conditions.static_stability.roll_rate
YAWQ = conditions.static_stability.yaw_rate
VINF = conditions.freestream.velocity
#freestream 0 velocity safeguard
if not conditions.freestream.velocity.all():
if settings.use_surrogate:
velocity = conditions.freestream.velocity
velocity[velocity==0] = np.ones(len(velocity[velocity==0])) * 1e-6
conditions.freestream.velocity = velocity
else:
raise AssertionError("VLM requires that conditions.freestream.velocity be specified and non-zero")
# Unpack vortex distribution
CHORD = VD.chord_lengths
chord_breaks = VD.chordwise_breaks
span_breaks = VD.spanwise_breaks
RNMAX = VD.panels_per_strip
LE_ind = VD.leading_edge_indices
ZETA = VD.tangent_incidence_angle
RK = VD.chordwise_panel_number
exposed_leading_edge_flag = VD.exposed_leading_edge_flag
YAH = VD.YAH*1.
YBH = VD.YBH*1.
XA1 = VD.XA1*1.
XB1 = VD.XB1*1.
# Compute X and Z BAR ouside of generate_vortex_distribution to avoid requiring x_m and z_m as inputs
VD.XBAR = np.ones(( len_mach,sum(LE_ind[0]))) * x_m
VD.ZBAR = np.ones(( len_mach,sum(LE_ind[0]))) * z_m
# ---------------------------------------------------------------------------------------
# STEP 10: Generate A and RHS matrices from VD and geometry
# ------------------ --------------------------------------------------------------------
# Compute flow tangency conditions
phi = np.arctan((VD.ZBC - VD.ZAC)/(VD.YBC - VD.YAC)) # dihedral angle
delta = np.arctan((VD.ZC - VD.ZCH)/((VD.XC - VD.XCH))) # mean camber surface angle
# Build the RHS vector
rhs = compute_RHS_matrix(VD,delta,phi,conditions,settings,geometry,pwm)
RHS = rhs.RHS*1 # this matches numpy=1.26 in terms of dimension
ONSET = rhs.ONSET*1
# Build induced velocity matrix, C_mn
C_mn, s, RFLAG, EW = compute_wing_induced_velocity(VD,mach,compute_EW=True)
# Turn off sonic vortices when Mach>1
RHS = RHS*RFLAG
# To ensure compatibility for np.linalg.solve across numpy1.0 and numpy2.0
RHS = np.atleast_3d(RHS)
# Build Aerodynamic Influence Coefficient Matrix
use_VORLAX_induced_velocity = settings.use_VORLAX_matrix_calculation
if not use_VORLAX_induced_velocity:
A = np.multiply(C_mn[:,:,:,0],np.atleast_3d(np.sin(delta)*np.cos(phi))) \
+ np.multiply(C_mn[:,:,:,1],np.atleast_3d(np.cos(delta)*np.sin(phi))) \
- np.multiply(C_mn[:,:,:,2],np.atleast_3d(np.cos(phi)*np.cos(delta))) # validated from book eqn 7.42
else:
A = EW
# Compute vortex strength
GAMMA = np.linalg.solve(A,RHS)
# To ensure compatibility for np.linalg.solve across numpy1.0 and numpy2.0
RHS = RHS.squeeze(axis=2)
GAMMA = GAMMA.squeeze(axis=2)
# ---------------------------------------------------------------------------------------
# STEP 11: Compute Pressure Coefficient
# ------------------ --------------------------------------------------------------------
#VORLAX subroutine = PRESS
# spanwise strip exposure flag, always 0 for RCAIDE's infinitely thin airfoils. Needs to change if thick airfoils added
RJTS = 0
# COMPUTE FREE-STREAM AND ONSET FLOW PARAMETERS. Used throughout the remainder of VLM
B2 = np.tile((mach**2 - 1),VD.n_cp[0])
SINALF = np.sin(aoa)
COSALF = np.cos(aoa)
TANALF = np.tan(aoa)
SINPSI = np.sin(PSI)
COPSI = np.cos(PSI)
COSIN = COSALF *SINPSI *2.0
COSINP = COSALF *SINPSI
COSCOS = COSALF *COPSI
PITCH = PITCHQ /VINF
ROLL = ROLLQ /VINF
YAW = YAWQ /VINF
# reshape CHORD
dim_1 = len(np.sum(LE_ind, axis=1))
dim_2 = np.sum(LE_ind, axis=1)[0]
# COMPUTE EFFECT OF SIDESLIP on DCP intermediate variables. needs change if cosine chorwise spacing added
FORAXL = COSCOS
FORLAT = COSIN
TAN_LEi= (VD.XB1[:,LE_ind[0]]-VD.XA1[:,LE_ind[0]])/ np.sqrt((VD.ZB1[:,LE_ind[0]]-VD.ZA1[:,LE_ind[0]])**2 + (VD.YB1[:,LE_ind[0]]-VD.YA1[:,LE_ind[0]])**2)
TAN_TE = (VD.XB_TE - VD.XA_TE)/ np.sqrt((VD.ZB_TE-VD.ZA_TE)**2 + (VD.YB_TE-VD.YA_TE)**2)
TAN_LE = np.repeat( TAN_LEi, RNMAX[LE_ind].reshape(dim_1,dim_2)[0] , axis=1)
TAN_LE = TAN_LE
TNL = TAN_LE * 1 # VORLAX's SIGN variable not needed, as these are taken directly from geometry
TNT = TAN_TE * 1
XIA = np.broadcast_to((RK-1)/RNMAX, np.shape(B2))
XIB = np.broadcast_to((RK )/RNMAX, np.shape(B2))
TANA = TNL *(1. - XIA) + TNT *XIA
TANB = TNL *(1. - XIB) + TNT *XIB
# cumsum GANT loop if KTOP > 0 (don't actually need KTOP with vectorized arrays and np.roll)
GFX = VD.chord_lengths
GANT = strip_cumsum(GFX*GAMMA, chord_breaks[0], RNMAX[LE_ind].reshape(dim_1,dim_2)[0] )
GANT = np.roll(GANT,1)
GANT[LE_ind] = 0
GLAT = GANT *(TANA - TANB) - GFX *GAMMA *TANB
cos_DL = (YBH-YAH)[LE_ind].reshape(dim_1,dim_2)/VD.D
COS_DL = np.repeat( cos_DL, RNMAX[LE_ind].reshape(dim_1,dim_2)[0] , axis=1)
DCPSID = FORLAT * COS_DL *GLAT /(XIB - XIA)
FACTOR = FORAXL + ONSET
# COMPUTE LOAD COEFFICIENT
GNET = GAMMA*FACTOR
GNET = GNET *RNMAX /CHORD
DCP = 2*GNET + DCPSID
CP = DCP
# ---------------------------------------------------------------------------------------
# STEP 12: Compute aerodynamic coefficients
# ------------------ --------------------------------------------------------------------
# Flip coordinates on the other side of the wing
boolean = YBH<0.
XA1[boolean], XB1[boolean] = XB1[boolean], XA1[boolean]
YAH[boolean], YBH[boolean] = YBH[boolean], YAH[boolean]
# Leading edge sweep. VORLAX does it panel by panel. This will be spanwise.
TLE = TAN_LE[LE_ind].reshape(dim_1,dim_2)
B2_LE = B2[LE_ind].reshape(dim_1,dim_2)
T2 = TLE*TLE
STB = np.zeros_like(B2_LE)
STB[B2_LE<T2] = np.sqrt(T2[B2_LE<T2]-B2_LE[B2_LE<T2])
# DL IS THE DIHEDRAL ANGLE (WITH RESPECT TO THE X-Y PLANE) OF
# THE IR STREAMWISE STRIP OF HORSESHOE VORTICES.
COD = np.cos(phi[LE_ind]).reshape(dim_1,dim_2) # Just the LE values
SID = np.sin(phi[LE_ind]).reshape(dim_1,dim_2) # Just the LE values
# Now on to each strip
PION = 2.0 /RNMAX
ADC = 0.5*PION
# XLE = LOCATION OF FIRST VORTEX MIDPOINT IN FRACTION OF CHORD.
XLE = 0.125 *PION
GAF = 0.5 + 0.5 *RJTS**2
# CORMED IS LENGTH OF STRIP CENTERLINE BETWEEN LOAD POINT
# AND TRAILING EDGE THIS PARAMETER IS USED IN THE COMPUTATION
# OF THE STRIP ROLLING COUPLE CONTRIBUTION DUE TO SIDESLIP.
X = VD.XCH #x-coord of load point (horseshoe centroid)
XTE = (VD.XA_TE + VD.XB_TE)/2 #Trailing edge x-coord behind the control point
CORMED = XTE - X
# SINF REFERENCES THE LOAD CONTRIBUTION OF IRT-VORTEX TO THE
# STRIP NOMINAL AREA, I.E., AREA OF STRIP ASSUMING CONSTANT
# (CHORDWISE) HORSESHOE SPAN.
SINF = ADC * DCP # The horshoe span lengths have been removed since VST/VSS == 1 always
# Split into chordwise strengths and sum into strips
# SICPLE = COUPLE (ABOUT STRIP CENTERLINE) DUE TO SIDESLIP.
CNC = np.add.reduceat(SINF ,chord_breaks[0],axis=1)
SICPLE = np.add.reduceat(SINF*CORMED,chord_breaks[0],axis=1)
# COMPUTE SLOPE (TX) WITH RESPECT TO X-AXIS AT LOAD POINTS BY INTER
# POLATING BETWEEN CONTROL POINTS AND TAKING INTO ACCOUNT THE LOCAL
# INCIDENCE.
XX = (RK - .75) *PION /2.0
TX = VD.SLOPE - ZETA
CAXL = -SINF*TX/(1.0+TX**2) # These are the axial forces on each panel
BMLE = (XLE-XX)*SINF # These are moment on each panel
# Sum onto the panel
CAXL = np.add.reduceat(CAXL,chord_breaks[0],axis=1)
BMLE = np.add.reduceat(BMLE,chord_breaks[0],axis=1)
SICPLE *= (-1) * COSIN * COD * GAF
DCP_LE = DCP[LE_ind].reshape(dim_1,dim_2)
# COMPUTE LEADING EDGE THRUST COEFF. (CSUC) BY CALCULATING
# THE TOTAL INDUCED FLOW AT THE LEADING EDGE. THIS COMPUTATION
# ONLY PERFORMED FOR COSINE CHORDWISE SPACING (LAX = 0).
# ** TO DO ** Add cosine spacing (earlier in VLM) to properly capture the magnitude of these earlier.
# Right now, this computation still happens with linear spacing, though its effects are underestimated.
CLE = compute_rotation_effects(VD, settings, EW, GAMMA, X, CHORD, XLE, VD.XBAR, rhs, COSINP, SINALF,COSCOS, PITCH, ROLL, YAW, STB, RNMAX)
# Leading edge suction multiplier. See documentation. This is a negative integer if used
# Default to 1 unless specified otherwise
SPC = K_SPC*np.ones_like(DCP_LE)
# If the vehicle is subsonic and there is vortex lift enabled then SPC changes to -1
VL = np.repeat(VD.vortex_lift,VD.n_sw[0], axis=1)
m_b = np.atleast_2d(mach[:,0]<1.)
SPC_cond = VL*m_b.T
SPC[SPC_cond] = -1.
SPC = SPC * exposed_leading_edge_flag
CLE = CLE + 0.5* DCP_LE *np.sqrt(XLE[LE_ind].reshape(dim_1,dim_2))
CSUC = 0.5*np.pi*np.abs(SPC)*(CLE**2)*STB
# TFX AND TFZ ARE THE COMPONENTS OF LEADING EDGE FORCE VECTOR ALONG
# ALONG THE X AND Z BODY AXES.
SLE = VD.SLOPE[LE_ind].reshape(dim_1,dim_2)
ZETA = ZETA[LE_ind].reshape(dim_1,dim_2)
XCOS = np.cos(SLE-ZETA)
XSIN = np.sin(SLE-ZETA)
TFX = 1.*XCOS
TFZ = -1.*XSIN
# If a negative number is used for SPC a different correction is used. See VORLAX documentation for Lan reference
TFX[SPC<0] = XSIN[SPC<0]*np.sign(DCP_LE)[SPC<0]
TFZ[SPC<0] = np.abs(XCOS)[SPC<0]*np.sign(DCP_LE)[SPC<0]
CAXL = CAXL - TFX*CSUC
# Add a dimension into the suction to be chordwise
CNC = CNC + CSUC*np.sqrt(1+T2)*TFZ
# FCOS AND FSIN ARE THE COSINE AND SINE OF THE ANGLE BETWEEN
# THE CHORDLINE OF THE IR-STRIP AND THE X-AXIS
FCOS = np.cos(ZETA)
FSIN = np.sin(ZETA)
# BFX, BFY, AND BFZ ARE THE COMPONENTS ALONG THE BODY AXES
# OF THE STRIP FORCE CONTRIBUTION.
BFX = - CNC *FSIN + CAXL *FCOS
BFY = - (CNC *FCOS + CAXL *FSIN) *SID
BFZ = (CNC *FCOS + CAXL *FSIN) *COD
# CONVERT CNC FROM CN INTO CNC (COEFF. *CHORD).
CHORD_strip = CHORD[LE_ind].reshape(dim_1,dim_2)
CNC = CNC * CHORD_strip
BMLE = BMLE * CHORD_strip
# BMX, BMY, AND BMZ ARE THE COMPONENTS ALONG THE BODY AXES
# OF THE STRIP MOMENT (ABOUT MOM. REF. POINT) CONTRIBUTION.
X = VD.XCH[LE_ind].reshape(dim_1,dim_2) # These are all LE values
Y = VD.YCH[LE_ind].reshape(dim_1,dim_2) # These are all LE values
Z = VD.ZCH[LE_ind].reshape(dim_1,dim_2) # These are all LE values
BMX = BFZ * Y - BFY * (Z - VD.ZBAR)
BMX = BMX + SICPLE
BMY = BMLE * COD + BFX * (Z - VD.ZBAR) - BFZ * (X - VD.XBAR)
BMZ = BMLE * SID - BFX * Y + BFY * (X - VD.XBAR)
CDC = BFZ * SINALF + (BFX *COPSI + BFY *SINPSI) * COSALF
CDC = CDC * CHORD_strip
ES = 2*s[:,0,:][LE_ind].reshape(dim_1,dim_2)
STRIP = ES *CHORD_strip
LIFT = (BFZ *COSALF - (BFX *COPSI + BFY *SINPSI) *SINALF)*STRIP
MOMENT = STRIP * (BMY *COPSI - BMX *SINPSI)
FY = (BFY *COPSI - BFX *SINPSI) *STRIP
RM = STRIP *(BMX *COSALF *COPSI + BMY *COSALF *SINPSI + BMZ *SINALF)
YM = STRIP *(BMZ *COSALF - (BMX *COPSI + BMY *SINPSI) *SINALF)
# Lift coefficient
Clift_y = LIFT/CHORD_strip/ES
CL_wing = np.add.reduceat(LIFT,span_breaks[0],axis=1)/VD.wing_areas
CLift = np.atleast_2d(np.sum(LIFT,axis=1)/S_ref).T
# Drag coefficient
results = compute_trefftz_plane_induced_drag(conditions, VD,Clift_y, X, Y, Z, CHORD_strip,S_ref,b_ref)
# force coefficeints
CX_for = (TANALF * CLift - results.CDrag_induced)/(COSALF - SINALF*TANALF)
CZ_for = (results.CDrag_induced+ CX_for*COSALF)/SINALF
CY_for = np.atleast_2d(np.sum(FY,axis=1)/S_ref).T
# moment coefficients
CM_mom = np.atleast_2d(np.sum(MOMENT,axis=1)/S_ref).T/c_ref
CL_mom = np.atleast_2d(np.sum(RM,axis=1)/S_ref).T /b_ref
CN_mom = np.atleast_2d(np.sum(YM,axis=1)/S_ref).T /b_ref
# ---------------------------------------------------------------------------------------
# STEP 13: Pack outputs
# ------------------ --------------------------------------------------------------------
results.CLift = CLift
results.CX = CX_for
results.CY = CY_for
results.CZ = -CZ_for
results.CL = CL_mom
results.CM = CM_mom
results.CN = -CN_mom
results.spanwise_stations = Y
results.CLift_wing = CL_wing
results.sectional_CLift = Clift_y
results.CP = np.array(CP , dtype=settings.floating_point_precision )
results.gamma = np.array(GAMMA , dtype=settings.floating_point_precision )
results.V_distribution = rhs.V_distribution
results.V_x = rhs.Vx_ind_total
results.V_z = rhs.Vz_ind_total
i = 0
dim_wing_lifts = results.CLift_wing * VD.wing_areas
dim_wing_drags = results.CDrag_induced_wing * VD.wing_areas
Clift_wings = Data()
Cdrag_wings = Data()
# Assign the lift and drag and non-dimensionalize
for wing in geometry.wings.values():
ref = wing.areas.reference
if wing.xz_plane_symmetric:
Clift_wings[wing.tag] = np.atleast_2d(np.sum(dim_wing_lifts[:,i:(i+2)],axis=1)).T/ref
Cdrag_wings[wing.tag] = np.atleast_2d(np.sum(dim_wing_drags[:,i:(i+2)],axis=1)).T/ref
i+=1
else:
Clift_wings[wing.tag] = np.atleast_2d(dim_wing_lifts[:,i]).T/ref
Cdrag_wings[wing.tag] = np.atleast_2d(dim_wing_drags[:,i]).T/ref
i+=1
results.CLift_wings = Clift_wings
results.CDrag_induced_wings = Cdrag_wings
results.VD = VD
return results
# ----------------------------------------------------------------------
# CLE rotation effects helper function
# ----------------------------------------------------------------------
[docs]
def compute_rotation_effects(VD, settings, EW_large, GAMMA, X, CHORD, XLE, XBAR,
rhs, COSINP, SINALF,COSCOS, PITCH, ROLL, YAW, STB, RNMAX):
""" This computes the effects of the freestream and aircraft rotation rate on
CLE, the induced flow at the leading edge
Assumptions:
Several of the values needed in this calculation have been computed earlier and stored in VD
Normally, VORLAX skips the calculation implemented in this function for linear
chordwise spacing (the if statement below). However, since the trends are correct,
albeit underestimated, this calculation is being forced here.
"""
LE_ind = VD.leading_edge_indices
RNMAX = VD.panels_per_strip
dim_1 = len(np.sum(LE_ind, axis=1))
dim_2 = np.sum(LE_ind, axis=1)[0]
dim_3 = len(LE_ind[0])
# Computate rotational effects (pitch, roll, yaw rates) on LE suction
# pick leading edge strip values for EW and reshape GAMMA -> gamma accordingly
EW = EW_large[LE_ind, :].reshape(dim_1, dim_2, dim_3)
gamma = np.array(np.split(np.repeat(GAMMA, dim_2, axis=0), dim_1))
CLE = (EW*gamma).sum(axis=2)
# Up till EFFINC, some of the following values were computed in compute_RHS_matrix().
# EFFINC and ALOC are calculated the exact same way, except for the XGIRO term.
# LOCATE VORTEX LATTICE CONTROL POINT WITH RESPECT TO THE
# ROTATION CENTER (XBAR, 0, ZBAR). THE RELATIVE COORDINATES
# ARE XGIRO, YGIRO, AND ZGIRO.
XGIRO = X - CHORD*XLE - np.repeat( XBAR, RNMAX[LE_ind].reshape(dim_1,dim_2)[0] , axis=1)
YGIRO = rhs.YGIRO
ZGIRO = rhs.ZGIRO
# VX, VY, VZ ARE THE FLOW ONSET VELOCITY COMPONENTS AT THE LEADING
# EDGE (STRIP MIDPOINT). VX, VY, VZ AND THE ROTATION RATES ARE
# REFERENCED TO THE FREE STREAM VELOCITY.
VX = (COSCOS - PITCH*ZGIRO + YAW *YGIRO)
VY = (COSINP - YAW *XGIRO + ROLL *ZGIRO)
VZ = (SINALF - ROLL *YGIRO + PITCH*XGIRO)
# CCNTL, SCNTL, SID, and COD were computed in compute_RHS_matrix()
# EFFINC = COMPONENT OF ONSET FLOW ALONG NORMAL TO CAMBERLINE AT
# LEADING EDGE.
EFFINC = VX *rhs.SCNTL + VY *rhs.CCNTL *rhs.SID - VZ *rhs.CCNTL *rhs.COD
CLE = CLE - EFFINC[LE_ind].reshape(dim_1,dim_2)
CLE = np.where(STB > 0, CLE /RNMAX[LE_ind].reshape(dim_1,dim_2) /STB, CLE)
return CLE
# ----------------------------------------------------------------------
# Vectorized cumsum from indices
# ----------------------------------------------------------------------
[docs]
def strip_cumsum(arr, chord_breaks, strip_lengths):
""" Uses numpy to to compute a cumsum that resets along
the leading edge of every strip.
Assumptions:
chordwise_breaks always starts at 0
"""
cumsum = np.cumsum(arr, axis=1)
offsets = cumsum[:,chord_breaks-1]
offsets[:,0] = 0
offsets = np.repeat(offsets, strip_lengths, axis=1)
return cumsum - offsets
[docs]
def compute_trefftz_plane_induced_drag(conditions, VD, cl, x_dist, y_dist, z_dist, chord_dist, SREF, b_ref, v_inf=1):
"""Compute induced drag using a Trefftz-plane (far-field wake) analysis.
The Trefftz plane is a control surface placed infinitely far downstream,
where the trailing vortex wake has fully rolled up. In this plane the
induced drag can be evaluated as a line integral of the downwash times
the local lift, which is computationally cleaner and more accurate than
integrating surface pressures in the near field.
The procedure follows the classical vortex-lattice / lifting-line approach
used in codes such as AVL (Drela & Youngren) and VORLAX:
1. **Circulation distribution** – The bound-vortex circulation on each
spanwise strip is reconstructed from the sectional lift coefficient and
the Kutta–Joukowski theorem:
Γ = 0.5 · c · V∞ · cl · (DS / |ΔY|)
where DS = √(ΔY² + ΔZ²) is the true arc-length of each strip, and the
DS/|ΔY| factor corrects for dihedral (for a flat wing DS = |ΔY| and the
factor is unity).
2. **Shed vortices** – Horseshoe vortices are shed at the spanwise edges of
each strip. The strength of the shed (trailing) vortex at each node is
the *difference* in bound circulation between adjacent strips (Helmholtz
vortex theorem). Tip vortices close the horseshoe with the full
circulation at the tip. Symmetric wings are handled by zeroing the
inboard root vortex; asymmetric wings shed a root vortex of full strength.
3. **Induced velocity in the Trefftz plane** – Each shed vortex induces a
velocity at the centerpoint of every strip. The 2-D Biot–Savart kernel
for a semi-infinite straight vortex filament reduces to:
w = Γ / (4π r)
where r is the perpendicular distance from the filament to the field
point. The velocity vector is projected onto the outward normal of the
receiving strip so that dihedral is handled correctly.
4. **Induced angle of attack** – The local induced downwash angle is:
α_i = arctan(w_i / V∞) ≈ w_i / V∞
5. **Induced drag** – Integrating along the span (using the composite
trapezoid rule over arc-length) gives:
CDi = ∫ (α_i · cl · c / S_ref) ds
where s is the spanwise arc-length coordinate. This integral is
evaluated separately for each wing surface and then summed.
Assumptions:
- Incompressible, inviscid flow (Trefftz-plane analysis is exact only in
this limit; compressibility corrections are applied elsewhere).
- The freestream speed is normalised to v_inf = 1 unless explicitly
overridden; lift coefficients and circulation are consistent with this
normalisation.
- Strip widths are taken from the VLM panel geometry stored in VD; the
chordwise discretisation does not affect the Trefftz-plane result.
- Each wing is treated as an independent lifting surface when assembling
the shed-vortex strengths; mutual induction *between* wings is
accounted for in the induced-velocity loop over all panels.
- Symmetric wings (xz-plane symmetry) have their inboard root vortex
strength set to zero (the image on the other side is implicit).
Parameters
----------
conditions : RCAIDE.Framework.Mission.Common.Results
Flight conditions data structure. The following fields are accessed:
- ``conditions.aerodynamics.angles.alpha`` (n_cases, 1) [rad]
Angle-of-attack for each flight condition.
VD : RCAIDE.Framework.Analyses.Aerodynamics.Vortex_Lattice_Method.Data
Vortex-lattice geometry descriptor populated by the VLM geometry
build. Key arrays (all indexed by case ``k``):
- ``VD.Y[k]`` – global Y-coordinates of all VLM nodes
- ``VD.Z[k]`` – global Z-coordinates of all VLM nodes
- ``VD.n_sw[k]`` – number of spanwise strips per wing (n_wings,)
- ``VD.n_cw[k]`` – number of chordwise panels per strip (n_wings,)
- ``VD.normals[k]`` – unit panel normal vectors (n_panels, 3)
- ``VD.symmetric_wings[k]`` – boolean list; True if wing has xz symmetry
- ``VD.wing_areas`` – reference area of each wing (n_wings,)
cl : ndarray, shape (n_cases, n_strips)
Sectional lift coefficient at the chordwise leading-edge panel of each
spanwise strip, as computed by the near-field VLM pressure integration.
x_dist : ndarray, shape (n_cases, n_strips)
Chordwise (X) coordinate of the leading-edge control points [m].
Currently not used directly inside this function (kept for API
consistency and potential future extensions).
y_dist : ndarray, shape (n_cases, n_strips)
Spanwise (Y) coordinate of the leading-edge control points [m].
Currently not used directly inside this function.
z_dist : ndarray, shape (n_cases, n_strips)
Vertical (Z) coordinate of the leading-edge control points [m].
Currently not used directly inside this function.
chord_dist : ndarray, shape (n_cases, n_strips)
Local chord length at each spanwise strip [m].
SREF : float
Vehicle reference area [m²] used for non-dimensionalising forces.
b_ref : float
Vehicle reference span [m]. Reserved for future use (aspect-ratio
corrections, Oswald efficiency output, etc.).
v_inf : float, optional
Freestream speed [m/s]. Default is 1 (non-dimensional). All
circulation and induced-velocity calculations are proportional to
v_inf; the resulting drag coefficient is independent of this value
in the incompressible, linear limit.
Returns
-------
results : RCAIDE.Framework.Core.Data
Container with the following fields:
- ``CDrag_induced`` ndarray (n_cases, 1)
Total vehicle induced-drag coefficient.
- ``sectional_CDrag_induced`` ndarray (n_cases, n_strips)
Spanwise distribution of the induced-drag coefficient integrand
``α_i · cl`` at each strip (before chord-weighting and integration).
- ``CDrag_induced_wing`` ndarray (n_cases, n_wings)
Induced-drag coefficient contribution from each individual wing
surface (referenced to SREF).
- ``alpha_induced`` ndarray (n_cases, n_strips)
Local induced angle of attack at each strip [rad].
Notes
-----
**Arc-length vs. projected-span integration**
RCAIDE normalises the sectional lift coefficient by the arc-length area
``c × DS`` rather than the projected area ``c × ΔY``. The DS/|ΔY|
correction in the circulation formula compensates for this convention so
that Γ is always the true bound circulation.
**Integration direction**
The sign of the tip shed vortex determines whether the strip ordering runs
from root-to-tip or tip-to-root. The integrand is flipped as needed so
that the trapezoid integration always proceeds from smaller to larger
arc-length values, which is required for ``scipy.integrate.trapezoid``.
**Mutual induction**
The double loop over wing panels (i, j outer; l, m inner) computes the
full mutual induction matrix, including the self-induction of each wing
on itself as well as between-wing terms. This is equivalent to the
Trefftz-plane integration in AVL but performed in the body (not wind)
frame.
**Numerical singularity guard**
Spanwise strips with zero projected width (ΔY = 0) are given a small
floor value of 1e-6 m to avoid division-by-zero in the arc-length
correction. This situation can arise for winglet panels that are
oriented in the XZ plane.
References
----------
Drela, M. and Youngren, H., "AVL 3.36 User Primer," MIT, 2017.
Katz, J. and Plotkin, A., "Low-Speed Aerodynamics," 2nd ed.,
Cambridge University Press, 2001, Chap. 8.
Lan, C. E., "A Quasi-Vortex-Lattice Method in Thin Wing Theory,"
Journal of Aircraft, Vol. 11, No. 9, 1974, pp. 518–527.
"""
# Initialize results storage
alpha = conditions.aerodynamics.angles.alpha
n_cases = len(alpha)
CDi_total = np.zeros(n_cases)
CDi_wing = np.zeros((n_cases, len(VD.n_sw[0])))
Cd_i_distribution = np.zeros_like(cl)
alpha_i = np.zeros_like(cl)
for k in range(n_cases):
alpha = conditions.aerodynamics.angles.alpha [k]
n_wings = len(VD.n_sw[k])
divisions = np.cumsum(VD.n_sw[k])[:-1]
# Split into per-wing lists — no np.stack, wings may have different n_sw
cl_split = np.split(cl[k], divisions)
chord_split = np.split(chord_dist[k], divisions)
# VD.Y[k]/VD.Z[k] are flat node arrays: each wing w contributes (n_sw[w]+1)*(n_cw[w]+1)
# corner nodes stored row-major. Striding by (n_cw[w]+1) picks one point per spanwise
# station. Using the first wing's n_cw as a global stride is wrong when wings differ.
node_sizes = (VD.n_sw[k] + 1) * (VD.n_cw[k] + 1)
node_splits = np.cumsum(node_sizes)[:-1]
y_nodes_per_wing = np.split(VD.Y[k], node_splits)
z_nodes_per_wing = np.split(VD.Z[k], node_splits)
y_control_points = [seg[::(VD.n_cw[k][w] + 1)] for w, seg in enumerate(y_nodes_per_wing)]
z_control_points = [seg[::(VD.n_cw[k][w] + 1)] for w, seg in enumerate(z_nodes_per_wing)]
is_symmetric = np.array(VD.symmetric_wings[0], dtype=bool)
is_vertical = np.array(VD.vertical_wing[0], dtype=bool)
symmetric_wing_flags = np.concatenate([np.repeat(is_symmetric & ~is_vertical, 2), np.zeros(np.count_nonzero(~is_symmetric), dtype=bool)])[:n_wings]
wing_areas = (symmetric_wing_flags+1)*VD.wing_areas
# ------------------------------------------------------------------------------------------
# Trefftz Plane Drag — per-wing quantities built into lists
# ------------------------------------------------------------------------------------------
y_centerpoints = []
z_centerpoints = []
shed_vortices = []
for w in range(n_wings):
ycp = y_control_points[w]
zcp = z_control_points[w]
cl_w = cl_split[w]
ch_w = chord_split[w]
y_centerpoints.append((ycp[:-1] + ycp[1:]) / 2)
# The Trefftz plane is evaluated in body-axis coordinates, meaning no angle-of-attack
# rotation is applied to the Z positions. This is consistent with how AVL's TPFORC
# routine handles the Prandtl-Glauert transform (ALFAT = 0 in that context): the wake
# lies in the Y-Z body plane and the bound-vortex geometry is already expressed in
# body axes, so rotating Z by alpha would double-count the incidence angle.
z_centerpoints.append((zcp[:-1] + zcp[1:]) / 2)
# The sectional lift coefficient cl is normalised by arc-length area (chord × DS),
# where DS = sqrt(DY² + DZ²) is the true spanwise arc-length of each strip. The
# Kutta–Joukowski theorem requires the actual bound circulation Γ = 0.5·V·c·cl, so
# when cl was formed using DS rather than the projected width DY the circulation must
# be scaled back by DS/|DY| to recover the physical value. For a flat, unswept wing
# DS = |DY| everywhere and the factor is exactly 1, so this correction is invisible
# in that limit. For dihedral or winglet panels the correction can be significant.
dy_w = np.diff(ycp)
dz_w = np.diff(zcp)
# Guard against division by zero for panels whose projected span is zero (e.g. winglets
# oriented entirely in the XZ plane). A floor of 1e-6 m is negligibly small relative
# to any physical panel width and does not meaningfully affect the circulation value.
dy_w_g = np.where(dy_w == 0, 1e-6, dy_w)
DS_w = np.sqrt(dy_w_g**2 + dz_w**2)
circ_w = 0.5 * ch_w * v_inf * cl_w * (DS_w / np.abs(dy_w_g))
sign_w = np.sign(ycp[1] - ycp[0])
sv_w = np.zeros(len(ycp))
sv_w[1:-1] = sign_w * np.diff(circ_w)
sv_w[-1] = -sign_w * circ_w[-1]
if not symmetric_wing_flags[w]:
sv_w[0] = -sign_w * circ_w[0]
shed_vortices.append(sv_w)
# --- Induced velocity (vectorised over all wings simultaneously) ---
yc = np.concatenate(y_centerpoints)
zc = np.concatenate(z_centerpoints)
yp = np.concatenate(y_control_points)
zp = np.concatenate(z_control_points)
gamma = np.concatenate(shed_vortices)
sdy = np.concatenate([np.diff(y_control_points[w]) for w in range(n_wings)])
sdz = np.concatenate([np.diff(z_control_points[w]) for w in range(n_wings)])
sds = np.sqrt(sdy**2 + sdz**2)
dy_mat = yc[:, None] - yp[None, :]
dz_mat = zc[:, None] - zp[None, :]
dist2 = dy_mat**2 + dz_mat**2
dist2[dist2 == 0] = np.inf
numerator = sdz[:, None] * dz_mat + sdy[:, None] * dy_mat
induced_velocity_flat = (
(gamma[None, :] * np.sign(sdy[:, None]) * numerator
/ (4.0 * np.pi * sds[:, None] * dist2))
.sum(axis=1)
)
alpha_induced_flat = np.arctan(induced_velocity_flat / v_inf)
cl_flat = np.concatenate(cl_split)
cd_induced_flat = alpha_induced_flat * cl_flat
# --- Per-wing CDi integration ---
CDi = 0
offset = 0
for w in range(n_wings):
n_sw_w = len(y_centerpoints[w])
ycp_w = y_control_points[w]
zcp_w = z_control_points[w]
cd_w = cd_induced_flat[offset:offset + n_sw_w]
ch_w = chord_split[w]
ld_w = np.cumsum(np.sqrt(np.diff(ycp_w)**2 + np.diff(zcp_w)**2))
CDi_w = trapezoid(cd_w * ch_w / wing_areas[0][w], ld_w)
CDi_wing[k][w] = CDi_w
CDi += CDi_w
offset += n_sw_w
CDi_total[k] = CDi
Cd_i_distribution[k] = cd_induced_flat
alpha_i[k] = alpha_induced_flat
# Package results
results = Data()
results.CDrag_induced = CDi_total[:,np.newaxis]
results.sectional_CDrag_induced = Cd_i_distribution
results.CDrag_induced_wing = CDi_wing
results.alpha_induced = alpha_i
return results
