Newer
Older
from finstrip_wrapper.core.main import FINSTRIPWrapper
from PGL.components.airfoil import AirfoilShape
'''
A basic wrapper for FINSTRIP that provides three primary functionalities:
file input generation, execution and file output reading.
This class extends FINSTRIPWrapper doing blade specific setup
The following keys should be passed as dictionary, otherwise the standard
values are used.
keys
----
name: str
job name
elsize_factor: float
element length as fraction of arc length of airfoil surface
core_mats: list
list of material name strings to be identified as core
spanpos: float
span wise position of cross section
use_becas_stiffness: bool
use BECAS cross section properties if True
'''
def __init__(self, **kwargs):
super(FINSTRIPWrapperBlade, self).__init__(**kwargs)
# data base output name
self.name = 'blade'
# element length as fraction of arc length of airfoil surface
self.elsize_factor = 0.015
# list of material name strings to be identified as core
self.core_mats = []
self.spanpos = 0 # position along span of analysed cross-section
self.use_becas_stiffness = False
for k, w in kwargs.iteritems():
try:
setattr(self, k, w)
except:
pass
self._init_half_wavelength_steps()
def init_dicts(self, cs2d):
# dictionaries written by FINSTRIPBladeStability according to FUSED-Wind format
self.cs2d = cs2d
def init_becas_stiffness(self, csprops, csprops_ref):
# stiffness properties as described in BECASWrapper
self.csprops = csprops
self.csprops_ref = csprops_ref
def init_cs_load_cases(self, csloads_rcs, twist):
''' csload is assumed to be provided in root cs
:param csloads: cross section loads (Fx,Fy,Fz,Mx,My,Mz) size ((ncase,6))
:param twist: geometric twist angle of cross section w.r.t root coordinate system size ((1))
note: so far only bending moments mx and my are used
'''
self.ncase = len(csloads_rcs[:, 0])
self.csloads = np.zeros_like(csloads_rcs)
for case in range(self.ncase):
Mx_rcs = csloads_rcs[case, 3]
My_rcs = csloads_rcs[case, 4]
# transformation matrix from root cs to cross section cs
sin_twist = np.sin(np.deg2rad(twist))
cos_twist = np.cos(np.deg2rad(twist))
self.T_twist = np.matrix(
[[cos_twist, sin_twist], [-sin_twist, cos_twist]])
# transform from root cs to cross section cs
v_Mrcs = np.array([[Mx_rcs], [My_rcs]])
v_Mcs = self.T_twist * v_Mrcs
M_x = np.squeeze(np.asarray(v_Mcs[0]))
M_y = np.squeeze(np.asarray(v_Mcs[1]))
self.csloads[case, 3] = M_x
self.csloads[case, 4] = M_y
self.csloads[:, 3:] = self.csloads[:, 3:] * 1E3 # convert to Nmm
self.ncase = len(self.csloads[:, 0])
def setup(self):
''' Sets up the FINSTRIP execution.
'''
super(FINSTRIPWrapperBlade, self).setup()
self._det_maxlength_element()
self._write_finstrip_input_files()
self._write_finstrip_analysis_files()
def _det_maxlength_element(self):
''' calculate element length based on fraction of arc length of airfoil
surface
'''
af = AirfoilShape()
af.initialize(points=self.cs2d['coords'])
# maximum size of generated element in m
self.maxlength_element = af.smax * self.elsize_factor
def _write_finstrip_input_files(self):
super(FINSTRIPWrapperBlade, self)._write_finstrip_input_files()
self._write_file_buc()
def _write_file_buc(self):
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
out_str = []
n = '\n'
# write header
out_str.append(
'// Finstrip input file, created by FINSTRIPWrapper' + n)
# write vars converted to mm
for var_name, var_value in zip(['maxlength_element',
'lambda_steps',
'lambda_min',
'lambda_max'],
[self.maxlength_element * 1E3,
self.lambda_steps,
self.lambda_min * 1E3,
self.lambda_max * 1E3,
]):
out_str.append('%s, %0.2f' % (var_name, var_value) + n)
# write materials
mat_names = self.cs2d['materials'].keys()
matprops = self.cs2d['matprops']
def _write_iso(out_str, n, label, E1, nu12, rho):
out_str.append('iso, %s, %g, %g, %g' %
(label, E1, nu12, rho) + n)
def _write_ortho(out_str, n, label, E1, E2, nu12, G12, rho):
out_str.append('ortho, %s, %g, %g, %g, %g, %g' % (
label, E1, E2, nu12, G12, rho) + n)
def _write_core(out_str, n, label, E1, E2, nu12, G12, G13, G23, rho):
out_str.append('core, %s, %g, %g, %g, %g, %g,%g, %g' % (
label, E1, E2, nu12, G12, G13, G23, rho) + n)
for i, label in enumerate(mat_names):
E1 = matprops[i, 0] * 1E-6 # to N/mm*2
E2 = matprops[i, 1] * 1E-6 # to N/mm*2
E3 = matprops[i, 2] * 1E-6 # to N/mm*2
nu12 = matprops[i, 3]
nu13 = matprops[i, 4]
nu23 = matprops[i, 5]
G12 = matprops[i, 6] * 1E-6 # to N/mm*2
G13 = matprops[i, 7] * 1E-6 # to N/mm*2
G23 = matprops[i, 8] * 1E-6 # to N/mm*2
rho = matprops[i, 9] * 1E-9 # to kg/mm*3
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
if label in self.core_mats:
# mat_type = 'core'
_write_core(
out_str, n, label, E1, E2, nu12, G12, G13, G23, rho)
elif E1 == E2 and nu12 == nu13:
# mat_type = 'iso'
_write_iso(out_str, n, label, E1, nu12, rho)
else:
# mat_type = 'ortho'
_write_ortho(out_str, n, label, E1, E2, nu12, G12, rho)
# write blade name and span pos
self.spanpos = self.cs2d['s'] * 1E3 # to mm
out_str.append('blade, %s, %g' % (self.name, self.spanpos) + n)
# write first (and only) section
out_str.append('Section, %g' % (self.spanpos) + n)
# write coords
coords = self.cs2d['coords'] * 1E3 # to mm
for i, c in enumerate(coords):
out_str.append('k, %g, %g' % (c[0], c[1]) + n)
# write region layup and dps
for ireg, _ in enumerate(self.cs2d['regions']):
layers = self.cs2d['regions'][ireg]['layers']
thicknesses = self.cs2d['regions'][
ireg]['thicknesses'] * 1E3 # to mm
angles = self.cs2d['regions'][ireg]['angles']
dp0 = self.cs2d['DPs'][ireg]
dp1 = self.cs2d['DPs'][ireg + 1]
for l, t, a in zip(layers, thicknesses, angles):
out_str.append('Layer, %s, %g, %g, %g, %g' %
(l[:-2], t, a, dp0, dp1) + n)
# write web layup and dps
for ireg, _ in enumerate(self.cs2d['webs']):
layers = self.cs2d['webs'][ireg]['layers']
thicknesses = self.cs2d['webs'][
ireg]['thicknesses'] * 1E3 # to mm
angles = self.cs2d['regions'][ireg]['angles']
# check if web exists at this cross section
if not np.sum(thicknesses) == 0.0:
dp0 = self.cs2d['DPs'][self.cs2d['web_def'][ireg][0]]
dp1 = self.cs2d['DPs'][self.cs2d['web_def'][ireg][1]]
out_str.append('Shearweb, %g, %g' % (dp0, dp1) + n)
wdp0 = 0.0 # web layup does not change along width
wdp1 = 1.0
for l, t, a in zip(layers, thicknesses, angles):
out_str.append('Layer, %s, %g, %g, %g, %g' %
(l[:-2], t, a, wdp0, wdp1) + n)
if self.use_becas_stiffness:
x_e = self.csprops_ref[2] * 1E3 # to mm
y_e = self.csprops_ref[3] * 1E3 # to mm
I_xx = self.csprops_ref[7]
I_yy = self.csprops_ref[8]
I_xy = self.csprops_ref[9]
nu = self.csprops_ref[17]
# TODO: check if sign is correctly taken for I_xy!
# A = self.csprops_ref[15]
# Gere 12-30
theta = 0.5 * np.arctan(-2 * I_xy / (I_xx - I_yy)) + np.pi / 2
# Gere 12-32
# Gere 12-31a
#theta1 = 0.5 * np.arccos((I_xx - I_yy) / (2 * R))
# Gere 12-31b
# if not np.sign(theta) == np.sign(nu):
# raise RuntimeError('BECAS input not consistent')
config = {}
config['BECASWrapper'] = {}
config['BECASWrapper']['hawc2_FPM'] = True
from becas_wrapper.becas_bladestructure import set_cs_size
cs_size, _ = set_cs_size(config)
if cs_size == len(self.csprops):
#hawc2_FPM = True
AE = self.csprops[20] # in N
else:
#hawc2_FPM = False
E = self.csprops[8]
An = self.csprops[15]
# AE = A * E # to N/mm
EI_xx = E * I_xx * 1E3 # to Nmm**2/mm
EI_yy = E * I_yy * 1E3 # to Nmm**2/mm
EI_xy = E * I_xy * 1E3 # to Nmm**2/mm
out_str.append('farob, %g, %g, %g, %g, %g, %g' %
(x_e, y_e, AE, EI_xx, EI_yy, EI_xy) + n)
fid = open(self.file_input, 'w')
fid.write(line)
fid.close()
self._write_file_lcd()
out_str = []
for case in range(self.ncase):
fx = self.csloads[case, 0]
fy = self.csloads[case, 1]
fz = self.csloads[case, 2]
mx = self.csloads[case, 3]
my = self.csloads[case, 4]
mz = self.csloads[case, 5]
n = '\n'
out_str.append(
'loads, %s, %g, %g, %g, %g, %g, %g' % (self.spanpos, fx, fy, fz, mx, my, mz) + n)
fid = open(self.file_load, 'w')
fid.write(line)
fid.close()
def execute(self):
super(FINSTRIPWrapperBlade, self).execute()
def post(self):
super(FINSTRIPWrapperBlade, self).post()
self._read_rst_eigenvalue()
self._read_rst_mode_shape()
def _read_rst_eigenvalue(self):
self.f_eig = np.zeros((self.ncase, self.nmode, 2))
for case in range(self.ncase):
# proces_nr, section_nr, load_nr, wavelength, eigenvalue
self.f_eig[case, :, :] = np.loadtxt(
self.file_rst_eigen, skiprows=1)[case * self.nmode:(case + 1) * self.nmode, -2:]
# convert wavelength to m
self.f_eig[:, :, 0] = self.f_eig[:, :, 0] * 1E-3
def _read_rst_mode_shape(self):
# determine number of elements
self.nelm = len(np.loadtxt('R' +
'%05i' % (self.spanpos) +
'_LC' +
'%03d' % 0 +
'_mode'
'%03d' % 0 +
self.txt_sfx,
skiprows=2)[:, 0])
# get mode shapes
self.mode_shape = np.zeros((self.ncase, self.nmode, self.nelm, 8))
# x-loc p1, y-loc p1, x-loc p2, y-loc p2,
# x-u p1, y-u p1, x-u p2, y-u p2
for case in range(self.ncase):
for mode in range(self.nmode):
# convert to m
self.mode_shape[case, mode, :, :] = 1E-3 * np.loadtxt('R' +
'%05i' % (self.spanpos) +
'_LC' +
'%03d' % (case) +
'_mode'
'%03d' % (mode) +
self.txt_sfx,
skiprows=2)
idxs = []
# identify surface and web nodes
for elm in range(self.nelm - 1):
elm1_node2 = self.mode_shape[0, 0, elm, 2:4]
elm2_node1 = self.mode_shape[0, 0, elm + 1, :2]
if elm1_node2[0] != elm2_node1[0] and elm1_node2[1] != elm2_node1[1]:
idxs.append(elm)
# append last elm index
idxs.append(len(self.mode_shape[0, 0, :, 0]) - 1)
self.nnode_surf = idxs[0] + 1 + 1
self.node_loc_surf = np.zeros((self.nnode_surf, 2))
:-1, :] = self.mode_shape[0, 0, :idxs[0] + 1, :2]
# add second node of last element
self.node_loc_surf[-1:,
:] = self.mode_shape[0, 0, idxs[0], 2:4]
# actual amount of webs in post processing files can differ from input
# NOTE: if the flatback web is oriented such that it closes the cross section,
# then it is identified as if the web was part of the surface
self.nweb = len(idxs[1:])
self.nnode_web = []
for w in range(self.nweb):
self.nnode_web.append(idxs[w + 1] - idxs[w] + 1)
self.nnode_web_max = np.max(np.array(self.nnode_web))
self.node_loc_web = np.zeros(
(self.nweb, self.nnode_web_max, 2))
for w in range(self.nweb):
self.node_loc_web[
w, :self.nnode_web[w] - 1, :] = self.mode_shape[0, 0, idxs[w] + 1:idxs[w + 1] + 1, :2]
# add second node of last element
self.node_loc_web[w, self.nnode_web[w] - 1, :] = self.mode_shape[
0, 0, idxs[w + 1], 2:4]
# obtain surface displacements
self.node_u_surf = np.zeros(
(self.ncase, self.nmode, self.nnode_surf, 2))
:, :, :-1, :] = self.mode_shape[:, :, :idxs[0] + 1, 4:6]
:, :, -1, :] = self.mode_shape[:, :, idxs[0], 6:]
# obtain web displacements
self.node_u_web = np.zeros(
(self.ncase, self.nweb, self.nmode, self.nnode_web_max, 2))
:, w, :, :self.nnode_web[w] - 1, :] = self.mode_shape[:, :, idxs[w] + 1:idxs[w + 1] + 1, 4:6]
self.node_u_web[:, w, :, self.nnode_web[w] - 1, :] = self.mode_shape[
:, :, idxs[w + 1], 6:]
plot = False
if plot:
import matplotlib.pylab as plt
plt.figure()
plt.plot(
self.node_loc_surf[:, 0], self.node_loc_surf[:, 1], '-o')
for w in range(self.nweb):
plt.plot(
self.node_loc_web[w, :self.nnode_web[w], 0],
self.node_loc_web[w, :self.nnode_web[w], 1], '-o')
mode = 15
sf = 1.0E+03
plt.plot(
self.node_loc_surf[:, 0] + sf *
self.node_u_surf[mode, :, 0],
self.node_loc_surf[:, 1] + sf * self.node_u_surf[mode, :, 1], '--')
for w in range(self.nweb):
plt.plot(
self.node_loc_web[
w, :self.nnode_web[w], 0] + sf * self.node_u_web[w, mode, :self.nnode_web[w], 0],
self.node_loc_web[w, :self.nnode_web[w], 1] + sf * self.node_u_web[w, mode, :self.nnode_web[w], 1], '--')