Commit e9e06cf9 authored by Frederik Zahle's avatar Frederik Zahle
Browse files

Merge branch 'tait_bryan_revisited' into 'master'

Tait bryan revisited

See merge request frza/PGL!12
parents 14383d35 85d71b55
......@@ -2,10 +2,10 @@
import numpy as np
from scipy.interpolate import pchip
import time
from PGL.main.domain import Domain, Block
from PGL.main.curve import Curve
from PGL.main.geom_tools import RotMat, transformation_matrix, dotX
from PGL.main.geom_tools import RotMat, transformation_matrix, dotX,\
calculate_length
from PGL.components.airfoil import AirfoilShape, BlendAirfoilShapes
......@@ -98,9 +98,8 @@ class LoftedBladeSurface(object):
self.dist_LE = np.array([])
self.gf_heights = np.array([])
self.shear_sweep = False
self.rotorder_sweep = False
self.grad_axis = np.array([])
self.theta = np.array([])
self.rotorder_sweep = True
self.analytic_grad = True
self.X_bar_unit = np.array([])
self.Y_bar_unit = np.array([])
self.Z_bar_unit = np.array([])
......@@ -256,40 +255,48 @@ class LoftedBladeSurface(object):
# final X
X = np.zeros((self.ni_chord, self.ni_span, 3), dtype = float)
# create the curve and get the main axis points in body CS
main_axis = Curve(points=np.array([self.x, self.y, self.z]).T)
# store the gradients of the axis
self.grad_axis = main_axis.dp
axis = Curve(points=np.array([self.x, self.y, self.z]).T,
spline='pchip')
# get gradient
if self.analytic_grad:
self.grad = self._analytic_gradient(axis)
else:
self.grad = axis.dp
# create the rotation matrix based on order
# rot_order is in the order of actual transformation
# ex: Prebend-Sweep-Twist will have the rotation order:[2, 1, 0]
rot_order = self.rot_order
# specify rotation order with either prebend first or sweep first
if self.rotorder_sweep and not self.shear_sweep:
if self.rotorder_sweep:
# get V_theta_direction
Vtheta_unit = self._rotation_direction()
# make rot_order=[2, 0, 1]: twist, prebend, sweep
rot_order[2] = 1
rot_order[1] = 0
# obtain angles:
# prebend
rot_x = np.arcsin(-self.grad_axis[:, 1])
rot_x = np.arcsin(-self.grad[:, 1])
# sweep
rot_y = np.arctan2(self.grad_axis[:, 0],
self.grad_axis[:, 2] + 1.e-20)
rot_y = np.arctan2(self.grad[:, 0],
self.grad[:, 2] + 1.e-20)
#shear-sweep
if self.shear_sweep:
rot_y[:] = 0
else:
else:
# prebend
rot_x = np.arctan2(-self.grad[:, 1],
self.grad[:, 2] + 1.e-20)
# sweep
rot_y = np.arcsin(self.grad_axis[:, 0])
rot_y = np.arcsin(self.grad[:, 0])
# shear-sweep ie local x is parallel to global X before twisting
if self.shear_sweep:
rot_y[:] = 0
# prebend
rot_x = np.arctan2(-self.grad_axis[:, 1],
self.grad_axis[:, 2] + 1.e-20)
# twist
rot_z = self.rot_z*np.pi/180.
# axis for rotation
......@@ -309,6 +316,8 @@ class LoftedBladeSurface(object):
p = np.zeros(3*self.ni_chord, dtype = float)
# construct the profile array
ind_p = np.arange(0, 3*self.ni_chord, step = 3, dtype = int)
#
twist = np.zeros((self.ni_span), dtype = float)
for i in range(self.ni_span):
# Rx: prebend
......@@ -316,15 +325,31 @@ class LoftedBladeSurface(object):
# Ry: sweep
Rot[1, :, :] = RotMat(rot_axis[1, :], rot_y[i])
# Rz: twist
Rot[2, :, :] = RotMat(rot_axis[2, :], rot_z[i])
# Rotation matrix
R = np.dot(Rot[rot_order[2], :, :],
np.dot(Rot[rot_order[1], :, :], Rot[rot_order[0], :, :]))
R = np.dot(Rot[rot_order[2], :, :], Rot[rot_order[1], :, :])
#
if self.rotorder_sweep:
twist_corr = self._twist_correction(R, Vtheta_unit[i])
# Rotation matrix after applying twist correction
Rot[2, :, :] = RotMat(rot_axis[2, :], twist_corr)
# corect the twist
twist[i] = rot_z[i] + twist_corr
else:
twist[i] = rot_z[i]
# Rz:
Rot[2, :, :] = RotMat(rot_axis[2, :], twist[i])
#Rotation matrix final
R = np.dot(R, Rot[rot_order[0], :, :])
# store the local X, Y and Z axes
X_bar_unit[i, :] = R[:, 0]
Y_bar_unit[i, :] = R[:, 1]
Z_bar_unit[i, :] = R[:, 2]
#
# construct the transformation matrix for the cross-section
T = transformation_matrix(X_bar_unit[i, :], Y_bar_unit[i, :],
Z_bar_unit[i, :], self.ni_chord)
......@@ -346,13 +371,84 @@ class LoftedBladeSurface(object):
X[:, i, 1] = Pb[:, 1] + self.y[i]
X[:, i, 2] = Pb[:, 2] + self.z[i]
# save the twist_correction
self.twist_corrected = twist
# save the local cross-sectional axes
self.X_bar_unit = X_bar_unit
self.Y_bar_unit = Y_bar_unit
self.Z_bar_unit = Z_bar_unit
return X
return X
def _twist_correction(self, R, Vtheta_unit):
"""
Returns the angles by which the twist has to be corrected.
"""
# V_theta in local co-ordinates = inv(R)*V_theta_unit
Vtheta_local = np.linalg.solve(R, Vtheta_unit)
# twist correction angle
twist_corr = np.arctan2(Vtheta_local[1], Vtheta_local[0])
return twist_corr
def _rotation_direction(self):
"""
Returns the direction of the local velocity vector tangential to the
rotor undergoing pure rotation without external inflow.
"""
# position vectors
P = np.array([self.x, self.y, self.z]).T
P[0,:] += 1.e-20
# unit vectors
P_unit = np.zeros((P.shape), dtype = float)
P_norm = np.linalg.norm(P, axis = 1)
P_unit[:,0] = np.divide(P[:,0], P_norm) #x
P_unit[:,1] = np.divide(P[:,1], P_norm) #y
P_unit[:,2] = np.divide(P[:,2], P_norm) #z
# Rotation vector
R_unit = np.zeros((P.shape), dtype = float)
R_unit[:, 1] = 1.
# direction of V_theta
Vtheta = np.cross(R_unit, P_unit, axisa = 1, axisb = 1)
# Vtheta_unit vector
Vtheta_unit = np.zeros((Vtheta.shape), dtype = float)
Vtheta_norm = np.linalg.norm(Vtheta, axis = 1)
Vtheta_unit[:,0] = np.divide(Vtheta[:,0], Vtheta_norm)#x
Vtheta_unit[:,1] = np.divide(Vtheta[:,1], Vtheta_norm)#y
Vtheta_unit[:,2] = np.divide(Vtheta[:,2], Vtheta_norm)#z
return Vtheta_unit
def _analytic_gradient(self, axis):
"""
Peforms a pchip interpolation and rturns the uit direction vector
"""
s = axis.s
points = axis.points
splx = axis._splines[0] # store pchip spline function for x = f(s)
sply = axis._splines[1] # store pchip spline function for y = f(s)
splz = axis._splines[2] # store pchip spline function for z = f(s)
#
grad = np.zeros((points.shape), dtype= float)
grad_unit = np.zeros((points.shape), dtype= float)
#
grad[:,0] = splx(s, nu=1)
grad[:,1] = sply(s, nu=1)
grad[:,2] = splz(s, nu=1)
# calculate norm for each span section
grad_norm = np.linalg.norm(grad, axis=1)
# obtain the unit direction vector
grad_unit[:, 0] = np.divide(grad[:,0], grad_norm)
grad_unit[:, 1] = np.divide(grad[:,1], grad_norm)
grad_unit[:, 2] = np.divide(grad[:,2], grad_norm)
return grad_unit
def _redistribute(self, points, pos_z, i):
if self.redistribute_flag == False:
......
......@@ -28,10 +28,10 @@ class SlicedLoftedSurface(object):
surface: ndarray
The lofted blade surface as input having shape (chord_sections,
span_sections, 3)
Ns : int
ni_span : int
The required spanwise sections which are less than or equal to the
spanwise sections of the input surface
n_points: int
ni_slice: int
The number of points on the required slice. Or can be also regarded as
the number of chordwise sections, which are less than or equal to that
of the input surface.
......@@ -81,7 +81,6 @@ class SlicedLoftedSurface(object):
self.tol = 1.e-5 # tolerance for residual in [m]
self.include_tip = True # flag for including tip
self.blade_length = 1.0 # blade length
self.close_surface = False # Optional: flag to close the surface
self.alpha_tol = 1.e-6 # Optional : tolerance for min. alpha
self.span_low = 0 # Optional: starting span section for slicing
self.span_high = 0 # optional: end span section for slicing
......@@ -95,6 +94,7 @@ class SlicedLoftedSurface(object):
self.sliced_surface = np.array([])
self.flag_open = False
self.perturb = 1.e-5 # open the trailing edge
self.close_surface = False # Optional: flag to close the surface
def update(self):
"""
......@@ -109,15 +109,15 @@ class SlicedLoftedSurface(object):
t1 = time.time()
t_elapsed = t1 - t0
zc_vec = self.zc_vec
span_low = self.span_low
span_high = self.span_high
#delspan_ind = self.delspan_ind
print('\n Geometry built from Z = %3.2f m to Z = %3.2f m \n'\
%(zc_vec[span_low], zc_vec[span_high-1]))
print('\nGeometry built from Z = %3.4f m to Z = %3.4f m \n'\
%(self.zc_vec[self.span_low], self.zc_vec[self.span_high-1]))
print('Time taken = %3.4f s \n'%t_elapsed)
if self.delspan_ind:
print('Interpolated spanwise sections:\n')
for i,e in enumerate(self.delspan_ind):
print('%i) Section %i at Z = %3.4f m\n'%(i+1, e,
self.zc_vec[e]))
def main(self):
......@@ -139,8 +139,6 @@ class SlicedLoftedSurface(object):
n_points = int(self.ni_slice)
#blade length
blade_length = self.blade_length
# flag for endpoint
#flag_endpoint = self.include_tip
#-------------------------------------------------------------------
#--------------------------------------------------------------------
# pre-process
......@@ -276,16 +274,17 @@ class SlicedLoftedSurface(object):
R, Nc_orig,
n_points)
# print the health of the iteration
print('--------------------------------------------------------------')
print('\n Span location: S = %i, radius = %3.2f \n'%(i, zc))
print('\n Iteration = %i, dc = %3.7f, Main jac cond = %e'%(count, dc,
jac_main_cond))
print('\n Residual : R_max = %3.7f, R_norm = %3.7f,' \
'R_new/R_prev = %3.5f \n'%(R_max, R_norm, R_norm/R_norm_prev))
print('------------------------------------------------------')
print('\n Span location: S = %i, radius = %3.4f m\n'%(i, zc))
print('\n Iteration = %i, dc = %3.6f m, Main jac cond = %e'
%(count, dc, jac_main_cond))
print('\n Residual : R_max = %3.7f\n'\
' R_norm = %3.7f\n'\
' R_new/R_prev = %3.5f \n'
%(R_max, R_norm, R_norm/R_norm_prev))
print('\n Relaxation factor : alpha = %3.7f \n'%(alpha))
print('\n Delta vector : delta_norm= %3.7f \n'%(delta_norm))
print('--------------------------------------------------------------')
time.sleep(0.1)
print('------------------------------------------------------')
#increase count
count += 1
......@@ -299,8 +298,9 @@ class SlicedLoftedSurface(object):
# escape the section
delspan_ind.append(i)
# print
print('\n\n Skipping section %i, radius = %3.2f\n\n'%(i, zc))
time.sleep(0.1)
print('\n\n Skipping section %i, radius = %3.2f\n\n'
%(i, zc))
#time.sleep(0.1)
# set exit flag as False
exit_flag = 0
......@@ -602,8 +602,10 @@ class SlicedLoftedSurface(object):
param_map_in[i, j, 1] = t_ind[j]
return surface_in, param_map_in
def bilinear_surface(self, surface_orig, grid_s, grid_t, S, T):
"""
Performs bilinear interpolation on the input lofted surface in order
to obtain the slice in a Z = c plane.
......@@ -633,7 +635,6 @@ class SlicedLoftedSurface(object):
#---------------------- perform bilinear interpolation----------------------
#1) find the 4 known data points closest to the point to be interpolated
Ncs_desired= T.shape[0]
# Ns_desired= S.shape[0]
# stores the positions of the four neighbourhood points for each corresponding
# interpolant grid location
grid_map= np.empty((Ncs_desired), dtype= object)
......@@ -709,49 +710,49 @@ class SlicedLoftedSurface(object):
x2= grid_map[i][3][0]
y2= grid_map[i][3][1]
A= np.matrix([[1, x1, y1, x1*y1],
A= np.array([[1, x1, y1, x1*y1],
[1, x1, y2, x1*y2],
[1, x2, y1, x2*y1],
[1, x2, y2, x2*y2]])
#X- values
X= np.matrix([[val_map[i][0][0]],
X= np.array([[val_map[i][0][0]],
[val_map[i][1][0]],
[val_map[i][2][0]],
[val_map[i][3][0]]])
#Y-values
Y= np.matrix([[val_map[i][0][1]],
Y= np.array([[val_map[i][0][1]],
[val_map[i][1][1]],
[val_map[i][2][1]],
[val_map[i][3][1]]])
#Z-values
Z= np.matrix([[val_map[i][0][2]],
Z= np.array([[val_map[i][0][2]],
[val_map[i][1][2]],
[val_map[i][2][2]],
[val_map[i][3][2]]])
# Coefficient matrix for X-values
Bx= np.linalg.inv(A) * X
Bx= np.linalg.solve(A, X)
#Coefficient matrix for Y-values
By= np.linalg.inv(A) * Y
By= np.linalg.solve(A, Y)
#Coefficient matrix for Z-values
Bz= np.linalg.inv(A) * Z
Bz= np.linalg.solve(A, Z)
x_desired= S[i]
y_desired= T[i]
# X-value for the new cross-sectional surface
Q[i, 0]= (float(Bx[0]) + float(Bx[1])*x_desired +
float(Bx[2])*y_desired + float(Bx[3])*x_desired*y_desired)
Q[i, 0]= (Bx[0] + Bx[1]*x_desired +
Bx[2]*y_desired + Bx[3]*x_desired*y_desired)
# Y-value for the new cross-sectional surface
Q[i, 1]= (float(By[0]) + float(By[1])*x_desired +
float(By[2])*y_desired + float(By[3])*x_desired*y_desired)
Q[i, 1]= (By[0] + By[1]*x_desired +
By[2]*y_desired + By[3]*x_desired*y_desired)
# Z-value for the new cross-sectional surface
Q[i, 2]= (float(Bz[0]) + float(Bz[1])*x_desired +
float(Bz[2])*y_desired + float(Bz[3])*x_desired*y_desired)
Q[i, 2]= (Bz[0] + Bz[1]*x_desired +
Bz[2]*y_desired + Bz[3]*x_desired*y_desired)
return Q, grid_map, val_map
def _jacobian_Q(self, S, T, grid_map, val_map):
"""
Stores the jacobians of the obtained X,Y,Z points with respect to the
......@@ -830,7 +831,7 @@ class SlicedLoftedSurface(object):
dYdt= np.zeros(n_points, dtype= float)
dZds= np.zeros(n_points, dtype= float)
dZdt= np.zeros(n_points, dtype= float)
for i in range(n_points):
s= S[i]
t= T[i]
......@@ -925,11 +926,11 @@ class SlicedLoftedSurface(object):
Parameters
----------
Q: float
Numpy array of the x,y,z coordinate points for the sampled
surface in the s-t parametric space.
D: float
Numpy array of discrete lengths between consecutive points on the
slice boundary.
Floating point numpy array of the x,y,z coordinate points for the
sampled surface in the s-t parametric space.
D: ndarray
Floating point numpy array of discrete lengths between consecutive
points on the slice boundary.
n_D: int
Size of the discrete length array.
n_points: int
......@@ -937,10 +938,10 @@ class SlicedLoftedSurface(object):
Returns
-------
jac_d: float
A sparse array of size n_Dx3N where N is the number of points on
the slice, consisting of jacobians of the discrete lengths
of the slice boundary with respect to the cartesian space.
jac_d: ndarray
A floating point sparse array of size n_Dx3N where N is the number
of points on the slice, consisting of jacobians of the discrete
lengths of the slice boundary with respect to the cartesian space.
"""
# extract the x,y,z coordinates of the slice points
......@@ -1024,17 +1025,17 @@ class SlicedLoftedSurface(object):
Parameters
----------
dZds: float array
Gradient of cartesian co-ordinate Z wrt to the parametric space
variable S
dZds: ndarray
Floating point numpy array of gradient of cartesian co-ordinate
Z wrt to the parametric space variable S
dZdt: float array
Gradient of cartesian co-ordinate Z wrt to the parametric space
variable
dZdt: ndarray
Floating point numpy array of gradient of cartesian co-ordinate Z
wrt to the parametric space variable
jac_dp: sparse array
Jacobian matrix of the discrete lenghts between slice points wrt
the parametric space (S,T)
jac_dp: ndarray
Sparse array of Jacobian of the discrete lenghts between slice
points wrt the parametric space (S,T)
n_points: int
Number of points on the slice
......@@ -1048,9 +1049,9 @@ class SlicedLoftedSurface(object):
Returns
-------
jac_main: sparse array
Main jacobian matrix used in the Newton iteration of shape 2N X
(2N+1), where N is the number of points on the slice
jac_main: ndarray
Sparse array of main jacobian used in the Newton iteration of
shape 2N X (2N+1), where N is the number of points on the slice
"""
......@@ -1131,12 +1132,13 @@ class SlicedLoftedSurface(object):
Parameters
----------
D: float
A 1D numpy array of the distances in cartesian space.
Q: float
A 2D numpy array of the cartesian space.
T: float
A 1D numpy array of the t-coordinates in the parametric space.
D: ndarray
A floating point 1D numpy array of the distances in cartesian space.
Q: ndarray
A floating point 2D numpy array of the cartesian space.
T: ndarray
A floating point 1D numpy array of the t-coordinates in the
parametric space.
zc: float
Constant z plane where the cross-section is being constructed.
dc: float
......@@ -1149,8 +1151,8 @@ class SlicedLoftedSurface(object):
Returns
-------
R: float
Numpy vector of the residual terms.
R: ndarray
A floating point 1D numpy array of the residual terms.
"""
......@@ -1176,14 +1178,15 @@ class SlicedLoftedSurface(object):
Parameters
----------
Pk: float
A (2n+1) numpy array representing the state of the problem at the
kth stage.
jac_main: float
A (2n+1) X (2n+1) numpy array representing the jacobian of the
state Pk.
R: float
A (2n+1) numpy array representing the residual of the problem.
Pk: ndarray
A floating point (2n+1) numpy array representing the state of the
problem at the kth stage.
jac_main: ndarray
A floating point sparse array of size (2n+1) X (2n+1) representing
the jacobian of the state Pk.
R: ndarray
A floating point numpy array of size (2n+1) representing the
residual of the problem.
Nc_orig: int
The number of cross-sectional points in the input surface.
n_points: int
......@@ -1191,9 +1194,9 @@ class SlicedLoftedSurface(object):
Returns
-------
Pk1: float
A (2n+1) 1D numpy array representing the state at the (k+1)th
iteration.
Pk1: ndarray
A floating point numpy array of size (2n+1) representing the state
at the (k+1)th iteration.
R_norm : float
Second norm of the residual.
delta_norm : float
......@@ -1228,9 +1231,9 @@ class SlicedLoftedSurface(object):
Parameters
----------
Pk: float
Pk: ndarray
The state vector comprising of [Si,Ti, S(i+1), T(i+1)] values.
delta: float
delta: ndarray
The change in state vector for the newton step.
Returns
......@@ -1292,7 +1295,7 @@ class SlicedLoftedSurface(object):
Parameters
----------
surface_new: float
surface_new: ndarray
Cross-sectional surface definition of order NsXNcX3
span_ind: int
......@@ -1300,7 +1303,7 @@ class SlicedLoftedSurface(object):
Returns
-------
surf_concat: float
surf_concat: ndarray
Numpy array of interpolated cross-sectional co-ordinates at the
requested spanwise section
......@@ -1325,7 +1328,7 @@ class SlicedLoftedSurface(object):
indc_len = ind - ind_prev
indc_new = indc_prev + indc_len
# concatenate the surface
surf_concat[indc_prev : indc_new, :, :] = surface_new[ind_prev:ind, :, :]
surf_concat[indc_prev : indc_new, :, :] = surface_new[ind_prev:ind,:, :]
# update the previous indices
indc_prev = indc_new
ind_prev = ind + 1
......@@ -1342,7 +1345,7 @@ class SlicedLoftedSurface(object):
Parameters
----------
surface: float
surface: ndarray
Cross-sectional surface definition of order NsXNcX3
zc: float
The spanwise z-location where the x,y co-ordiantes are to be
......@@ -1354,7 +1357,7 @@ class SlicedLoftedSurface(object):
Returns
-------
surf: float
surf: ndarray
Interpolated cross-sectional co-ordinates at the requested section
"""
......
......@@ -59,7 +59,7 @@ m.blade_length = 1.0
# m.dist = dist
# Option 2: s (spanwise distribution of blade curve length)
# m.s = pf['s']
#m.s = pf['s']
#Option 3: Not setting a specific the above two, proceeds with linear distribution
#-----------------------------------------------------------------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment