# RCAIDE/Methods/Aerodynamics/Airfoil_Panel_Method/velocity_distribution.py
#
#
# Created: Dec 2023, M. Clarke
# ----------------------------------------------------------------------------------------------------------------------
# IMPORT
# ----------------------------------------------------------------------------------------------------------------------
# pacakge imports
import numpy as np
# ----------------------------------------------------------------------------------------------------------------------
# velocity_distribution
# ----------------------------------------------------------------------------------------------------------------------
[docs]
def velocity_distribution(qg,x,y,xbar,ybar,st,ct,alpha_2d,npanel,ncases,ncpts):
"""Compute the tangential velocity distribution at the
midpoint of each panel
Source:
None
Assumptions:
None
Inputs:
qg - Vector of source/sink and vortex strengths [unitless]
x - Vector of x coordinates of the surface nodes [unitless]
y - Vector of y coordinates of the surface nodes [unitless]
xbar - x-coordinate of the midpoint of each panel [unitless]
ybar - y-coordinate of the midpoint of each panel [unitless]
st - np.sin(theta) for each panel [radians]
ct - np.cos(theta) for each panel [radians]
al - Angle of attack in radians [radians]
npanel - Number of panels on the airfoil [unitless]
Outputs:
vt_2d - Vector of tangential velocities
Properties Used:
N/A
"""
# flow tangency boundary condition - source distribution
vt_2d = ct *np.cos(alpha_2d) + st*np.sin(alpha_2d)
gamma = np.repeat(qg[-1,:,:][np.newaxis,:,:],npanel, axis = 0)
# convert 1d matrices to 2d
qg_2d = np.repeat(qg[:-1,:,:][np.newaxis,:,:,:],npanel, axis = 0)
x_2d = np.repeat(np.swapaxes(np.swapaxes(x,0, 2),0,1)[:,:,np.newaxis,:],npanel, axis = 2)
y_2d = np.repeat(np.swapaxes(np.swapaxes(y,0, 2),0,1)[:,:,np.newaxis,:],npanel, axis = 2)
xbar_2d = np.repeat(np.swapaxes(np.swapaxes(xbar,0, 2),0,1)[:,:,:,np.newaxis],npanel, axis = 3)
ybar_2d = np.repeat(np.swapaxes(np.swapaxes(ybar,0, 2),0,1)[:,:,:,np.newaxis],npanel, axis = 3)
st_2d = np.repeat(np.swapaxes(np.swapaxes(st,0, 2),0,1)[:,:,:,np.newaxis],npanel, axis = 3)
ct_2d = np.repeat(np.swapaxes(np.swapaxes(ct,0, 2),0,1)[:,:,:,np.newaxis],npanel, axis = 3)
st_2d_T = np.swapaxes(st_2d,2,3)
ct_2d_T = np.swapaxes(ct_2d,2,3)
# Fill the elements of the matrix of aero influence coefficients
sti_minus_j = ct_2d_T*st_2d - st_2d_T*ct_2d
cti_minus_j = ct_2d_T*ct_2d + st_2d_T*st_2d
rij = np.sqrt((xbar_2d-x_2d[:,:,:,:-1])**2 + (ybar_2d-y_2d[:,:,:,:-1])**2)
rij_plus_1 = np.sqrt((xbar_2d-x_2d[:,:,:,1:])**2 + (ybar_2d-y_2d[:,:,:,1:])**2)
rij_dot_rij_plus_1 = (xbar_2d-x_2d[:,:,:,:-1])*(xbar_2d-x_2d[:,:,:,1:]) + (ybar_2d-y_2d[:,:,:,:-1])*(ybar_2d-y_2d[:,:,:,1:])
anglesign = np.sign((xbar_2d-x_2d[:,:,:,:-1])*(ybar_2d-y_2d[:,:,:,1:]) - (xbar_2d-x_2d[:,:,:,1:])*(ybar_2d-y_2d[:,:,:,:-1]))
r_ratio = rij_dot_rij_plus_1/rij/rij_plus_1
r_ratio[r_ratio>1.0] = 1.0 # numerical noise
betaij = np.real(anglesign*np.arccos(r_ratio))
# betaij[aoas,res,diag_indices,diag_indices] = np.pi
for i in range(ncases):
for j in range(ncpts):
np.fill_diagonal(betaij[i,j,:,:], np.pi)
# swap axes
sti_minus_j_2d = np.swapaxes(np.swapaxes(sti_minus_j,0,2),1,3)
betaij_2d = np.swapaxes(np.swapaxes(betaij,0,2),1,3)
cti_minus_j_2d = np.swapaxes(np.swapaxes(cti_minus_j,0,2),1,3)
rij_2d = np.swapaxes(np.swapaxes(rij,0,2),1,3)
rij_plus_1_2d = np.swapaxes(np.swapaxes(rij_plus_1,0,2),1,3)
vt_2d += np.sum(qg_2d/2/np.pi*(sti_minus_j_2d*betaij_2d - cti_minus_j_2d*np.log(rij_plus_1_2d/rij_2d)),1) + \
np.sum(gamma/2/np.pi*(sti_minus_j_2d*np.log(rij_plus_1_2d/rij_2d) + cti_minus_j_2d*betaij_2d),1)
return vt_2d
