'''
Construct airfoil with CST
'''
from typing import Tuple
import numpy as np
from numpy.linalg import lstsq
from scipy.special import factorial
from .math import (rotate, interp_from_curve, find_circle_3p,
cst_foil, clustcos, dist_clustcos, cst_curve)
from .basic import BasicSection
#* ===========================================
#* CST sections
#* ===========================================
[docs]
class Section(BasicSection):
'''
Section 3D curve generated by CST foil (upper & lower surface)
Examples
--------
>>> sec = Section(thick=None, chord=1.0, twist=0.0, tail=0.0, lTwistAroundLE=True)
Attributes
------------
tail : float
tail thickness (m)
RLE : float
relative leading edge radius
te_angle : float
trailing edge angle (degree)
te_slope : float
slope of the mean camber line at trailing edge (dy/dx)
cst_u, cst_l : ndarray
cst coefficients of the upper and lower surfaces
refine_u, refine_l : {None, ndarray}
cst coefficients of the refinement curve on upper and lower surfaces
'''
def __init__(self, thick=None, chord=1.0, twist=0.0, tail=0.0, lTwistAroundLE=True):
super().__init__(thick=thick, chord=chord, twist=twist, lTwistAroundLE=lTwistAroundLE)
self.tail = tail
self.RLE = 0.0
self.te_angle = 0.0 # trailing edge angle (degree)
self.te_slope = 0.0 # slope of the mean camber line at trailing edge (dy/dx)
#* 2D unit airfoil
self.cst_u = np.zeros(1)
self.cst_l = np.zeros(1)
#* Refine airfoil
self.refine_u = None
self.refine_l = None
[docs]
def section(self, cst_u=None, cst_l=None, nn=1001, flip_x=False, projection=True) -> None:
'''
Generating the section (3D) by `cst_foil`.
First construct the 2D unit curve, then transform it to the 3D curve.
Parameters
------------
nn : int
number of points in `xx`, `yy`, `yu`, and `yl`.
cst_u, cst_l : ndarray
CST coefficients of upper and lower surfaces (optional)
flip_x : bool
whether flip `xx` in the reverse order, by default False.
projection : bool
whether keeps the projection length the same when rotating the section, by default True.
Examples
------------
>>> sec.section(cst_u=None, cst_l=None, nn=1001, flip_x=False, projection=True)
'''
#* Update CST parameters
if isinstance(cst_u, np.ndarray) and isinstance(cst_l, np.ndarray):
self.cst_u = cst_u.copy()
self.cst_l = cst_l.copy()
#* Construct airfoil with CST parameters
self.xx, self.yu, self.yl, self.thick, self.RLE = cst_foil(
nn, self.cst_u, self.cst_l, t=self.specified_thickness, tail=self.tail)
#* Trailing edge information
a1 = [self.xx[-1]-self.xx[-5], self.yu[-1]-self.yu[-5]]
a2 = [self.xx[-1]-self.xx[-5], self.yl[-1]-self.yl[-5]]
self.te_angle = np.arccos(np.dot(a1,a2)/np.linalg.norm(a1)/np.linalg.norm(a2))/np.pi*180.0
self.te_slope = 0.5*((self.yu[-1]+self.yl[-1])-(self.yu[-5]+self.yl[-5]))/(self.xx[-1]-self.xx[-5])
#* Refine the airfoil by incremental curves
yu_i = np.zeros(nn)
yl_i = np.zeros(nn)
if isinstance(self.refine_u, np.ndarray):
_, y_tmp = cst_curve(nn, self.refine_u, x=self.xx)
yu_i += y_tmp
if isinstance(self.refine_l, np.ndarray):
_, y_tmp = cst_curve(nn, self.refine_l, x=self.xx)
yl_i += y_tmp
self.yu, self.yl = foil_increment_curve(self.xx, self.yu, self.yl, yu_i=yu_i, yl_i=yl_i, t=self.specified_thickness)
#* Transform to 3D
super().section(flip_x=flip_x, projection=projection)
[docs]
class OpenSection(BasicSection):
'''
Section 3D curve generated by CST curve (open curve)
Examples
--------
>>> sec = OpenSection(thick=None, chord=1.0, twist=0.0, lTwistAroundLE=True)
Attributes
------------
cst : ndarray
cst coefficients of the curve
refine: {None, ndarray}
cst coefficients of the refinement curve
'''
def __init__(self, thick=None, chord=1.0, twist=0.0, lTwistAroundLE=True):
super().__init__(thick=thick, chord=chord, twist=twist, lTwistAroundLE=lTwistAroundLE)
#* 2D unit curve
self.cst = np.zeros(1)
#* Refine airfoil
self.refine = None
[docs]
def section(self, cst=None, nn=1001, flip_x=False, projection=True):
'''
Generating the section (3D) by `cst_curve`.
First construct the 2D unit curve, then transform it to the 3D curve.
Parameters
------------
nn : int
number of points in `xx`, `yy`, `yu`, and `yl`.
cst : ndarray
CST coefficients of upper and lower surfaces (optional)
flip_x : bool
whether flip `xx` in the reverse order, by default False.
projection : bool
whether keeps the projection length the same when rotating the section, by default True.
Examples
------------
>>> sec.section(cst=None, nn=1001, flip_x=False, projection=True)
'''
#* Update CST parameters
if isinstance(cst, np.ndarray):
self.cst = cst.copy()
#* Construct curve with CST parameters
self.xx, self.yy = cst_curve(nn, self.cst)
#* Refine the geometry with an incremental curve
if isinstance(self.refine, np.ndarray):
_, y_i = cst_curve(nn, self.refine, x=self.xx)
self.yy += y_i
#* Apply thickness
self.thick = np.max(self.yy, axis=0)
if isinstance(self.specified_thickness, float):
self.yy = self.yy/self.thick*self.specified_thickness
self.thick = self.specified_thickness
#* Transform to 3D
super().section(flip_x=flip_x, projection=projection)
[docs]
class RoundTipSection(BasicSection):
'''
Section 3D curve generated by CST foil and base shape (upper & lower surface)
Suitable for round trailing edge foils or blades, or plates.
Parameters
-----------
xLE, yLE, zLE: float,
coordinates of the leading edge
chord: float
chord length (m)
thick: float
maximum relative thickness
twist: float
twist angle (deg)
tail: float
actual thickness of TE (m)
cst_u, cst_l : ndarray
cst coefficients of the upper and lower surfaces
base_le_radius: float
relative radius of base shape function leading edge
base_te_radius: float
relative radius of base shape function trailing edge
base_abs_thick: float
actual thickness of the base shape
base_le_ratio: float
ratio of the leading edge region
base_te_ratio: float
ratio of the trailing edge region
aLE: float
angle (deg) of the slope at leading edge (a>0 => dy/dx>0)
aTE: float
angle (deg) of the slope at trailing edge (a<0 => dy/dx<0)
i_split: {None, int}
active when leading edge and trailing edge curves are intersected
nn : int
number of points in `xx`, `yy`, `yu`, and `yl`.
lTwistAroundLE : bool
whether the twist center is LE, otherwise TE, by default True.
'''
def __init__(self, xLE: float, yLE: float, zLE: float,
chord: float, thick: float, twist: float, tail: float,
cst_u: np.ndarray, cst_l: np.ndarray,
base_le_ratio: float, base_te_ratio: float, base_abs_thick: float,
base_le_radius: float, base_te_radius: float,
aLE=0.0, aTE=0.0, i_split=None, nn=501, lTwistAroundLE=False):
super().__init__(thick=thick, chord=chord, twist=twist, lTwistAroundLE=lTwistAroundLE)
self.xLE = xLE
self.yLE = yLE
self.zLE = zLE
self.thick = thick
self.tail = tail
#* --------------------------------------------------
#* Base shape (actual chord length)
x_ref = dist_clustcos(nn, a0=0.0079, a1=0.96, beta=2)
if base_abs_thick > 0.0:
x_, y_ = RoundTipSection.base_shape(x_ref, xLE, xLE+chord, base_le_ratio*chord, base_te_ratio*chord,
base_le_radius, base_te_radius, base_abs_thick/2.0, i_split=i_split)
dy_ = RoundTipSection.base_camber(x_, a_LE=aLE, a_TE=aTE)
# Scale to unit chord length
self.xx = (x_ - xLE)/chord
self.base_yu = (dy_ + y_)/chord
self.base_yl = (dy_ - y_)/chord
else:
dy_ = RoundTipSection.base_camber(xLE+x_ref*chord, a_LE=aLE, a_TE=aTE)
self.xx = x_ref
self.base_yu = dy_/chord
self.base_yl = dy_/chord
#* Base shape thickness check
base_tmax = np.max(self.base_yu-self.base_yl)
if base_tmax > thick:
print('Warning: base shape is thicker than tmax (%.3f > %.3f).'%(base_tmax*chord, thick*chord))
print(' The base shape is scaled to tmax.')
self.base_yu = dy_ + y_*thick/base_tmax
self.base_yl = dy_ - y_*thick/base_tmax
if base_tmax > tail/chord and tail>0.0:
print('Warning: base shape is thicker than the specified tail thickness (%.3f > %.3f).'%(base_tmax*chord, tail))
print(' The final tail thickness is about the base thickness.')
#* --------------------------------------------------
#* CST shape (actual chord length)
self.cst_yu = np.zeros_like(self.xx)
self.cst_yl = np.zeros_like(self.xx)
if np.max(np.abs(cst_u))>1E-6 or np.max(np.abs(cst_l))>1E-6 > 0:
self.cst_u = cst_u.copy()
self.cst_l = cst_l.copy()
cst_tail = max(0, tail/chord-base_tmax)
_, self.cst_yu, self.cst_yl, _, _ = cst_foil(self.xx.shape[0], cst_u, cst_l,
x=self.xx, t=thick-base_tmax, tail=cst_tail)
#* --------------------------------------------------
#* Unit chord shape
self.yu = self.base_yu + self.cst_yu
self.yl = self.base_yl + self.cst_yl
#* Calculate leading edge radius
x_RLE = 0.005
yu_RLE = interp_from_curve(x_RLE, self.xx, self.yu)
yl_RLE = interp_from_curve(x_RLE, self.xx, self.yl)
self.RLE, _ = find_circle_3p([0.0,0.0], [x_RLE,yu_RLE], [x_RLE,yl_RLE])
[docs]
@staticmethod
def base_shape(x_ref: np.ndarray, x_LE: float, x_TE: float, l_LE: float, l_TE: float,
r_LE: float, r_TE: float, h: float, i_split=None) -> Tuple[np.ndarray, np.ndarray]:
'''
Base shape function of wing sections.
Parameters
------------
x_ref: ndarray [nn]
reference point distribution in [0,1]
x_LE: float
leading edge location
x_TE: float
trailing edge location
l_LE: float
length of leading edge curve/ramp
l_TE: float
length of trailing edge curve/ramp
r_LE: float
relative radius of leading edge
r_TE: float
relative radius of trailing edge
h: float
height
i_split: {None, int}
active when leading edge and trailing edge curves are intersected
Notes
-------
Base shape: \n
`_______________________________________________` \n
`(______________________________________________)` \n
Examples
---------
>>> x, y = base_shape(x_ref, x_LE, x_TE, l_LE, l_TE, r_LE, r_TE, h, i_split=None)
'''
l0 = x_TE - x_LE
if abs(l0)<=1e-10:
return np.ones_like(x_ref)*x_LE, np.zeros_like(x_ref)
x = x_ref*l0+x_LE
y = np.ones_like(x)*h
ratio_LE = l_LE/l0
ratio_TE = l_TE/l0
if l_LE+l_TE<=l0:
i_LE = np.argmin(np.abs(x_ref-ratio_LE))+1
i_TE = np.argmin(np.abs(x_ref-1+ratio_TE))-1
y[:i_LE] = RoundTipSection.general_eqn(x_ref[:i_LE], ratio_LE, r_LE, h)
y[i_TE:] = RoundTipSection.general_eqn(1-x_ref[i_TE:], ratio_TE, r_TE, h)
return x, y
else:
print('Warning: the specified length of LE and TE region %.2f > chord length %2f.'%(l_LE+l_TE, l0))
print(' The LE and TE curves will intersect.')
y_le = RoundTipSection.general_eqn(x_ref, ratio_LE, r_LE, h)
y_te = RoundTipSection.general_eqn(1-x_ref, ratio_TE, r_TE, h)
i_IT = np.argmin(np.abs(y_le-y_te))
#* Locate intersection point
x_m = x_ref[i_IT]
x_l = x_ref[i_IT-1] # y_le(x_l) < y_te(x_l)
x_r = x_ref[i_IT+1] # y_le(x_r) > y_te(x_r)
for _ in range(10):
x_m = 0.5*(x_l+x_r)
d_m = RoundTipSection.general_eqn(np.array([x_m]), ratio_LE, r_LE, h) \
- RoundTipSection.general_eqn(np.array([1-x_m]), ratio_TE, r_TE, h)
if d_m < -1e-10:
x_l = x_m
elif d_m > 1e-10:
x_r = x_m
else:
break
y_m = RoundTipSection.general_eqn(np.array([x_m]), ratio_LE, r_LE, h)[0]
if i_split == None:
y = np.concatenate((y_le[:i_IT], y_te[i_IT:]))
y[i_IT] = y_m
return x, y
else:
print(' The intersect point is specified by i_split = %d'%(i_split))
nn = x_ref.shape[0]
x_le = dist_clustcos(i_split, a0=0.01, a1=0.96, beta=2)*x_m
x_te = dist_clustcos(nn-i_split+1, a0=0.05, a1=0.96)*(1.0-x_m)+x_m
y_le = RoundTipSection.general_eqn(x_le, ratio_LE, r_LE, h)
y_te = RoundTipSection.general_eqn(1-x_te, ratio_TE, r_TE, h)
xx = np.concatenate((x_le[:, np.newaxis], x_te[1:, np.newaxis]), axis=0)
yy = np.concatenate((y_le[:, np.newaxis], y_te[1:, np.newaxis]), axis=0)
xx = xx[:,0]*l0+x_LE
yy = yy[:,0]
xx[i_split-1] = x_m*l0+x_LE
yy[i_split-1] = y_m
return xx, yy
[docs]
@staticmethod
def general_eqn(x: np.ndarray, l: float, rr: float, h: float) -> np.ndarray:
'''
General equations to define the leading edge semi-thickness,
the flat plate semi-thickness, the trailing edge closure semi-thickness,
and the transverse radius of the sting fairing.
Experimental Surface Pressure Data Obtained on 65° Delta Wing Across Reynolds Number
and Mach Number Ranges (Volume 2—Small-Radius Leading Edges)
https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19960025648.pdf
Parameters
------------
x: ndarray
current location
l: float
range of x
rr: float
relative radius
h: float
height
Examples
-----------
>>> phi = general_eqn(x, l, rr, t) # phi >= 0
'''
a = np.sqrt(2*rr*h)
b = -15/8.*a + 3*h
c = 5/4.*a - 3*h
d = -3/8.*a + h
xi = x/l
for i in range(xi.shape[0]):
xi[i] = max(0.0, min(1.0, xi[i]))
phi = a*np.sqrt(xi)+b*xi+c*np.power(xi,2)+d*np.power(xi,3)
return phi
[docs]
@staticmethod
def base_camber(x: np.ndarray, a_LE=0.0, a_TE=0.0) -> np.ndarray:
'''
Camber curve of the base shape function (in 3rd order spline)
Parameters
-----------
x: ndarray
actual x distribution
a_LE: float
angle (deg) of the slope at leading edge (a>0 => dy/dx>0)
a_TE: float
angle (deg) of the slope at trailing edge (a<0 => dy/dx<0)
Examples
----------
>>> dy = base_camber(x, a_LE, a_TE)
'''
chord = x[-1]-x[0]
x_ref = (x-x[0])/chord # reference point distribution in [0,1]
dy_LE = np.tan(a_LE/180.0*np.pi)*chord
dy_TE = np.tan(a_TE/180.0*np.pi)*chord
a = dy_LE + dy_TE
b = -2*dy_LE-dy_TE
c = dy_LE
dy = a*x_ref**3+b*x_ref**2+c*x_ref
return dy
#* ===========================================
#* Foil functions
#* ===========================================
[docs]
def normalize_foil(xu: np.ndarray, yu: np.ndarray, xl: np.ndarray, yl: np.ndarray):
'''
Transform the airfoil to a unit airfoil.
Parameters
----------
xu, yu, xl, yl: ndarray
coordinates of the airfoil
Returns
--------
xu_, yu_, xl_, yl_: ndarray
coordinates of the airfoil
twist: float
twist angle (degree)
chord: float
chord length
tail: float
tail height relative to chord length
Examples
---------
>>> xu_, yu_, xl_, yl_, twist, chord, tail = normalize_foil(xu, yu, xl, yl)
'''
if abs(xu[0]-xl[0])>1e-6 or abs(yu[0]-yl[0])>1e-6:
raise Exception('Two curves do not have the same leading edge')
#* Transform
xu_ = xu - xu[0]
xl_ = xl - xl[0]
yu_ = yu - yu[0]
yl_ = yl - yl[0]
#* Twist
xTE = 0.5*(xu_[-1]+xl_[-1])
yTE = 0.5*(yu_[-1]+yl_[-1])
twist = np.arctan(yTE/xTE)*180/np.pi
chord = np.sqrt(xTE**2+yTE**2)
xu_, yu_, _ = rotate(xu_, yu_, np.zeros_like(xu_), angle=-twist, axis='Z')
xl_, yl_, _ = rotate(xl_, yl_, np.zeros_like(xu_), angle=-twist, axis='Z')
#* Scale
yu_ = yu_ / xu_[-1]
yl_ = yl_ / xl_[-1]
xu_ = xu_ / xu_[-1]
xl_ = xl_ / xl_[-1]
#* Removing tail
tail = abs(yu_[-1]) + abs(yl_[-1])
for ip in range(yu_.shape[0]):
yu_[ip] -= xu_[ip]*yu_[-1]
for ip in range(yl_.shape[0]):
yl_[ip] -= xl_[ip]*yl_[-1]
return xu_, yu_, xl_, yl_, twist, chord, tail
[docs]
def scale_cst(x: np.ndarray, yu: np.ndarray, yl: np.ndarray, cst_u, cst_l, t: float, tail=0.0):
'''
Scale CST coefficients, so that the airfoil has the maximum thickness of t.
Parameters
-----------
x, yu, yl: ndarray [nn]
baseline airfoil. `x`, `yu`, `yl` must be directly generated by CST, without scaling.
cst_u, cst_l: list or ndarray
CST coefficients of the upper and lower surfaces
t: float
specified relative maximum thickness
tail: float
relative tail thickness (optional)
Examples
---------
>>> cst_u_new, cst_l_new = scale_cst(yu, yl, cst_u, cst_l, t)
'''
thick = yu - yl
it = np.argmax(thick)
t0 = thick[it]
r = (t-tail*x[it])/t0
cst_u_new = np.array(cst_u) * r
cst_l_new = np.array(cst_l) * r
return cst_u_new, cst_l_new
#* ===========================================
#* Fitting a curve/foil with CST
#* ===========================================
[docs]
def cst_foil_fit(xu: np.ndarray, yu: np.ndarray, xl: np.ndarray, yl: np.ndarray,
n_cst=7, xn1=0.5, xn2=1.0) -> Tuple[np.ndarray, np.ndarray]:
'''
Using CST method to fit an airfoil
Parameters
----------
xu, yu, xl, yl: ndarray
coordinates
n_cst: int
number of CST coefficients
xn1, xn2: float
CST parameters
Returns
-------
cst_u, cst_l: ndarray
CST coefficients
Examples
---------
>>> cst_u, cst_l = cst_foil_fit(xu, yu, xl, yl, n_cst=7, xn1=0.5, xn2=1.0)
Notes
-----
This function allows the airfoil has non-zero tail thickness.
Also allows the airfoil chord length not equals to one.
But yu[0] yl[0] should be 0.
'''
cst_u = fit_curve(xu, yu, n_cst=n_cst, xn1=xn1, xn2=xn2)
cst_l = fit_curve(xl, yl, n_cst=n_cst, xn1=xn1, xn2=xn2)
return cst_u, cst_l
[docs]
def fit_curve(x: np.ndarray, y: np.ndarray, n_cst=7, xn1=0.5, xn2=1.0):
'''
Using least square method to fit a CST curve.
Parameters
----------
x, y: ndarray
coordinates
n_cst: int
number of CST coefficients
xn1, xn2: float
CST parameters
Returns
-------
coef: ndarray [n_cst]
CST coefficients
Notes
-----
y[0] should be 0.
Examples
---------
>>> coef = fit_curve(x, y, n_cst=7, xn1=0.5, xn2=1.0)
'''
# Array A: A[nn, n_cst], nn=len(x).
# Array b: b[nn].
nn = x.shape[0]
L = x[-1] - x[0] # type: float
x_ = (x-x[0])/L # scaling x to 0~1
y_ = (y-y[0])/L # scaling according to L #! This means y[0] should be 0
b = y_.copy()
for ip in range(nn):
b[ip] -= x_[ip]*y_[-1] # removing tail
A = np.zeros((nn, n_cst))
for ip in range(nn):
C_n1n2 = np.power(x_[ip],xn1) * np.power(1-x_[ip],xn2)
for i in range(n_cst):
xk_i_n = factorial(n_cst-1)/factorial(i)/factorial(n_cst-1-i)
A[ip][i] = xk_i_n * np.power(x_[ip],i) * np.power(1-x_[ip],n_cst-1-i) * C_n1n2
solution = lstsq(A, b, rcond=None)
return solution[0]
[docs]
def fit_curve_with_twist(x: np.ndarray, y: np.ndarray, n_cst=7,
xn1=0.5, xn2=1.0) -> Tuple[np.ndarray, float, float, float]:
'''
Using least square method to fit a CST curve
Parameters
----------
x, y: ndarray
coordinates
n_cst: int
number of CST coefficients
xn1, xn2: float
CST parameters
Returns
-------
coef: ndarray [n_cst]
CST coefficients
chord: float
distance between two ends of the curve
twist: float
degree, +z axis
thick: float
maximum relative thickness
Examples
---------
>>> coef, chord, twist, thick = fit_curve_with_twist(x, y, n_cst, xn1, xn2)
'''
chord = np.sqrt((x[0]-x[-1])**2+(y[0]-y[-1])**2)
twist = np.arctan((y[-1]-y[0])/(x[-1]-x[0]))*180/np.pi
x_ = (x - x[0])/chord
y_ = (y - y[0])/chord
x_, y_, _ = rotate(x_, y_, np.zeros_like(x_), angle=-twist, axis='Z')
thick = np.max(y_, axis=0)
coef = fit_curve(x_, y_, n_cst=n_cst, xn1=xn1, xn2=xn2)
return coef, chord, twist, thick
[docs]
def fit_curve_partial(x: np.ndarray, y: np.ndarray, n_cst=7, ip0=0, ip1=0,
ic0=0, ic1=0, xn1=0.5, xn2=1.0):
'''
Using least square method to fit a part of a unit curve
Parameters
----------
x, y: ndarray
coordinates
n_cst: int
number of CST coefficients
ip0, ip1: int
index of the partial curve x[ip0:ip1]
ic0, ic1: int
index of the CST parameters cst[ic0:ic1] that are not 0
xn1, xn2: float
CST parameters
Returns
-------
coef: ndarray [n_cst]
CST coefficients
Examples
---------
>>> coef = fit_curve_partial(x: np.array, y: np.array, ip0=0, ip1=0, n_cst=7, xn1=0.5, xn2=1.0)
'''
# Array A: A[nn, n_cst], nn=len(x).
# Array b: b[nn].
ip0 = max(0, ip0)
if ip1 <= ip0:
ip1 = x.shape[0]
ic0 = max(0, ic0)
if ic1 <= ic0:
ic1 = n_cst
#* Fit the partial curve
A = np.zeros((ip1-ip0, ic1-ic0))
for ip in range(ip0, ip1):
C_n1n2 = np.power(x[ip],xn1) * np.power(1-x[ip],xn2)
for i in range(ic0,ic1):
xk_i_n = factorial(n_cst-1)/factorial(i)/factorial(n_cst-1-i)
A[ip-ip0][i-ic0] = xk_i_n * np.power(x[ip],i) * np.power(1-x[ip],n_cst-1-i) * C_n1n2
solution = lstsq(A, y[ip0:ip1], rcond=None)
coef = np.zeros(n_cst)
for i in range(ic0, ic1):
coef[i] = solution[0][i-ic0]
return coef
#* ===========================================
#* Modification of a curve/foil
#* ===========================================
[docs]
def foil_bump_modify(x: np.ndarray, yu: np.ndarray, yl: np.ndarray,
xc: float, h: float, s: float, side=1, n_cst=0,
return_cst=False, keep_tmax=True):
'''
Add bumps on the airfoil
Parameters
----------
x, yu, yl: ndarray
coordinates of the airfoil
xc: float
x of the bump center
h: float
relative height of the bump (to maximum thickness)
s: float
span of the bump
side: int
+1/-1 upper/lower side of the airfoil
n_cst: int
if specified (>0), then use CST to fit the new foil
return_cst: bool
if True, also return cst_u, cst_l when n_cst > 0
keep_tmax: bool
if True, keep the maximum thickness unchanged scale the opposite side of 'side' to keep thickness
Returns
-------
yu_new, yl_new: ndarray
coordinates
cst_u, cst_l: ndarray
CST coefficients
Examples
---------
>>> yu_new, yl_new (, cst_u, cst_l) = foil_bump_modify(
>>> x: np.array, yu: np.array, yl: np.array,
>>> xc: float, h: float, s: float, side=1,
>>> n_cst=0, return_cst=False, keep_tmax=True)
'''
yu_new = yu.copy()
yl_new = yl.copy()
t0 = np.max(yu_new-yl_new)
if xc<0.1 or xc>0.9:
kind = 'H'
else:
kind = 'G'
if side > 0:
yu_new = yu_new + bump_function(x, xc, h*t0, s, kind=kind)
else:
yl_new = yl_new + bump_function(x, xc, h*t0, s, kind=kind)
if keep_tmax:
it = np.argmax(yu_new-yl_new)
tu = np.abs(yu_new[it])
tl = np.abs(yl_new[it])
#* Scale the opposite side
if side > 0:
rl = (t0-tu)/tl
yl_new = rl * np.array(yl_new)
else:
ru = (t0-tl)/tu
yu_new = ru * np.array(yu_new)
t0 = None
if n_cst > 0:
# CST reverse
tail = yu[-1] - yl[-1]
cst_u, cst_l = cst_foil_fit(x, yu_new, x, yl_new, n_cst=n_cst)
_, yu_new, yl_new, _, _ = cst_foil(x.shape[0], cst_u, cst_l, x=x, t=t0, tail=tail)
else:
cst_u = None
cst_l = None
if return_cst:
return yu_new, yl_new, cst_u, cst_l
else:
return yu_new, yl_new
[docs]
def foil_increment(x, yu, yl, cst_u=None, cst_l=None, t=None) -> Tuple[np.ndarray, np.ndarray]:
'''
Add cst curve by incremental curves
Parameters
----------
x, yu, yl : ndarray
coordinates of the baseline airfoil.
cst_u, cst_l : ndarray or None
if not None, CST parameters of the incremental curves.
The number of CST parameters is arbitrary.
t : float or None
relative maximum thickness.
Returns
-------
yu_new, yl_new : ndarray
coordinates of the new airfoil.
Examples
---------
>>> yu_new, yl_new = foil_increment(x, yu, yl, cst_u, cst_l, t=None)
'''
nn = len(x)
if cst_u is not None:
_, yu_i = cst_curve(nn, cst_u, x=x)
else:
yu_i = None
if cst_u is not None:
_, yl_i = cst_curve(nn, cst_l, x=x)
else:
yl_i = None
yu_new, yl_new = foil_increment_curve(x, yu, yl, yu_i=yu_i, yl_i=yl_i, t=t)
return yu_new, yl_new
[docs]
def foil_increment_curve(x: np.ndarray, yu: np.ndarray, yl: np.ndarray,
yu_i=None, yl_i=None, t=None) -> Tuple[np.ndarray, np.ndarray]:
'''
Add cst curve by incremental curves
Parameters
-----------
x, yu, yl : ndarray
coordinates of the baseline airfoil.
yu_i, yl_i : {None, ndarray}
if not None, coordinates of the incremental curves.
t : {None, float}
relative maximum thickness.
Returns
---------
yu_, yl_ : ndarray
coordinates of the new airfoil.
Examples
---------
>>> yu_, yl_ = foil_increment_curve(x, yu, yl, yu_i, yl_i, t=None)
'''
nn = len(x)
if not isinstance(yu_i, np.ndarray):
yu_i = np.zeros(nn)
if not isinstance(yl_i, np.ndarray):
yl_i = np.zeros(nn)
x_ = x.copy()
yu_ = yu.copy()
yl_ = yl.copy()
# Remove tail
tail = yu_[-1] - yl_[-1]
if tail > 0.0:
yu_ = yu_ - 0.5*tail*x_
yl_ = yl_ + 0.5*tail*x_
# Add incremental curves
yu_ = yu_ + yu_i
yl_ = yl_ + yl_i
thick = yu_-yl_
it = np.argmax(thick)
t0 = thick[it]
# Apply thickness constraint
if t is not None:
r = (t-tail*x[it])/t0
yu_ = yu_ * r
yl_ = yl_ * r
# Add tail
if tail > 0.0:
yu_ = yu_ + 0.5*tail*x_
yl_ = yl_ - 0.5*tail*x_
return yu_, yl_
[docs]
def bump_function(x: np.ndarray, xc: float, h: float, s: float, kind='G') -> np.ndarray:
'''
A bump distribution [x, y].
Parameters
-----------
x, y: ndarray
current curve, x in [0,1]
xc: float
x of the bump center
h: float
height of the bump
s: float
span of the bump
kind: str
bump function type.
'G': Gaussian, less cpu cost
'H': Hicks-Henne, better when near leading edge
Returns
---------
y_bump : ndarray
coordinates
Examples
---------
>>> y_bump = bump_function(x, xc, h, s, kind)
'''
y_bump = np.zeros_like(x)
if xc<=0 or xc>=1:
print('Bump location not valid (0,1): xc = %.3f'%(xc))
return y_bump
if 'G' in kind:
sigma = np.ones_like(x) * (s / 6.)
if xc - s < 0.:
sigma = np.where(x < xc, xc / 3.5, sigma)
if xc + s > 1.:
sigma = np.where(x > xc, (1.0 - xc) / 3.5, sigma)
aa = -np.power(x - xc, 2) / 2.0 / sigma**2
y_bump += h * np.exp(aa)
else:
s0 = np.log(0.5)/np.log(xc)
Pow = 1
span = 1.0
hm = np.abs(h)
while Pow < 100 and span > s:
xx = np.linspace(0, 1, 201)
rr = np.pi * np.power(xx, s0)
yy = hm * np.power(np.sin(rr), Pow)
bump_range = xx[yy > 0.01 * hm]
span = bump_range[-1] - bump_range[0]
Pow = Pow + 1
rr = np.pi * np.power(x, s0)
y_bump += h * np.power(np.sin(rr), Pow)
return y_bump