From 025f62ca5b2fefb2eb2d2f0d5d43f99b6e8dbce6 Mon Sep 17 00:00:00 2001
From: Malo Rosemeier <malo.rosemeier@iwes.fraunhofer.de>
Date: Thu, 20 Jul 2017 17:22:24 +0200
Subject: [PATCH] wrapper works, verification of cross section props required.
---
finstrip_wrapper/core/blade.py | 330 ++++-
finstrip_wrapper/core/main.py | 24 +-
finstrip_wrapper/omdao/finstrip_blade.py | 1143 +++++++++--------
.../omdao/test/test_finstrip_blade.py | 130 +-
4 files changed, 1020 insertions(+), 607 deletions(-)
diff --git a/finstrip_wrapper/core/blade.py b/finstrip_wrapper/core/blade.py
index c2fcc38..6a15e97 100644
--- a/finstrip_wrapper/core/blade.py
+++ b/finstrip_wrapper/core/blade.py
@@ -1,5 +1,6 @@
import numpy as np
from finstrip_wrapper.core.main import FINSTRIPWrapper, status
+from PGL.components.airfoil import AirfoilShape
class FINSTRIPWrapperBlade(FINSTRIPWrapper):
@@ -12,18 +13,10 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
# data base output name
self.name = 'blade'
- # loads dummy
- self.loads = np.array([])
- #self.aligns = np.array([])
+ self.core_mats = [] # list of materials to be identified as core
+ self.spanpos = 0 # position along span of analysed cross-section
- # flags
- #self.post_rst = True
-
- # standard buckling analysis variables
- self.buckle_modes = 10 # number of bucling modes
- self.buckle_load_start = 0.0 # start boundary load
- self.buckle_load_end = 10.0 # end boundary load
- # self.buckle_method = 'LANB' # 'SUBSP'
+ self.use_becas_stiffness = False
for k, w in kwargs.iteritems():
try:
@@ -31,56 +24,248 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
except:
pass
+ self._init_half_wavelength_steps()
self._init_finstrip_input_files_names()
- def init_dicts(self, pf):
+ def init_dicts(self, cs2d):
# dictionaries written by FINSTRIPBladeStructure according to FUSED-Wind format
# NOTE: absolut dimensions are required!
- self.pf = pf
-
- def init_full_scale_load_case(self, loads, aligns):
- # TODO: here we need the conversion from load distro to point loads
- # and the calculation of the load introduction positions
- # By now it takes only external point loads
-
- # TODO: remove hardcoded loads
- # position of master nodes and force vector,and actuator positions
- # x,y,z,rotx,roty,rotz,Fx,Fy,Fz,Mx,My,Mz,ax,ay,az
- # actuator positions in the test_tig locinate system
- self.loads = loads
- self.nloads = self.loads.shape[0]
-
- # tilt, pitch, x_orig, y_orig, z_orig
- # tilt_angle rotation about root x-axis in deg
- # pitch_angle rotation about root z-axis in deg
- # test_rig_origin use distance in y-direction for vertical pull
- self.aligns = aligns
+ 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
+ '''
+
+ Mx_rcs = csloads_rcs[:, 3]
+ My_rcs = csloads_rcs[:, 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 = np.zeros_like(csloads_rcs)
+
+ self.csloads[:, 3] = M_x
+ self.csloads[:, 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 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):
- pass
+ 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
+ # deterimne material type
+ 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] + np.pi / 2
+ # TODO: check if sign is correct in I_xy!
+ # A = self.csprops_ref[15]
+
+ # Gere 12-30
+ theta = 0.5 * np.arctan(-2 * I_xy / (I_xx - I_yy))
+
+ # Gere 12-32
+ R = np.sqrt((0.5 * (I_xx - I_yy))**2 + I_xy**2)
+
+ # Gere 12-31a
+ theta1 = 0.5 * np.arccos((I_xx - I_yy) / (2 * R))
+
+ # Gere 12-31b
+ theta2 = 0.5 * np.arcsin(-I_xy / R)
+
+ # 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] * 1E-3 # to N/mm
+ else:
+ #hawc2_FPM = False
+ E = self.csprops[8]
+ An = self.csprops[15]
+ AE = E * An * 1E-3 # to N/mm
+
+ # 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')
+ for i, line in enumerate(out_str):
+ fid.write(line)
+ fid.close()
def _write_finstrip_analysis_files(self):
- self._write_file_loads()
+ self._write_file_lcd()
def _write_file_lcd(self):
- '''
- np.savetxt(self.file_loads, self.loads, fmt='%.10e', delimiter=',')
- np.savetxt(
- self.file_aligns, self.aligns, fmt='%.10e', delimiter=',')
- '''
+ 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'
+ # write header
+ out_str.append(
+ 'loads, %s, %i, %g, %g, %g, %g, %g' % (self.spanpos, fx, fy, fz, mx, my, mz) + n)
+
+ fid = open(self.file_load, 'w')
+ for i, line in enumerate(out_str):
+ fid.write(line)
+ fid.close()
def execute(self):
super(FINSTRIPWrapperBlade, self).execute()
@@ -91,16 +276,19 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
self._read_rst_mode_shape()
def _read_rst_eigenvalue(self):
- # proces_nr, section_nr, load_nr, wavelength, eigenvalue
- self.f_eig = np.loadtxt(self.file_rst_eigen, skiprows=1)
- # convert wavelength to m
- self.f_eig[:, -2] = self.f_eig[:, -2] * 1E-3
+ 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' +
- '%i' % (self.spanpos) +
+ '%05i' % (self.spanpos) +
'_LC' +
'%03d' % 0 +
'_mode'
@@ -108,41 +296,43 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
self.txt_sfx,
skiprows=2)[:, 0])
# get mode shapes
- self.mode_shape = np.zeros((self.nmode, self.nelm, 8))
+ 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):
- # transform to m
- self.mode_shape[mode, :, :] = 1E-3 * np.loadtxt('R' +
- '%i' % (self.spanpos) +
- '_LC' +
- '%03d' % (case) +
- '_mode'
- '%03d' % (mode) +
- self.txt_sfx,
- skiprows=2)
+ self.mode_shape[case, mode, :, :] = 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, elm, 2:4]
- elm2_node1 = self.mode_shape[0, elm + 1, :2]
+ 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]) - 1)
+ idxs.append(len(self.mode_shape[0, 0, :, 0]) - 1)
# obtain surface nodes
- self.nnode_loc_surf = idxs[0] + 1 + 1
- self.node_loc_surf = np.zeros((self.nnode_loc_surf, 2))
+ self.nnode_surf = idxs[0] + 1 + 1
+ self.node_loc_surf = np.zeros((self.nnode_surf, 2))
self.node_loc_surf[
- :-1, :] = self.mode_shape[0, :idxs[0] + 1, :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, idxs[0], 2:4]
+ :] = self.mode_shape[0, 0, idxs[0], 2:4]
# obtain web nodes
+ # 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 = []
@@ -154,29 +344,29 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
(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, idxs[w] + 1:idxs[w + 1] + 1, :2]
+ 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, idxs[w + 1], 2:4]
+ 0, 0, idxs[w + 1], 2:4]
# obtain surface displacements
self.node_u_surf = np.zeros(
- (self.nmode, self.nnode_loc_surf, 2))
+ (self.ncase, self.nmode, self.nnode_surf, 2))
self.node_u_surf[
- :, :-1, :] = self.mode_shape[:, :idxs[0] + 1, 4:6]
+ :, :, :-1, :] = self.mode_shape[:, :, :idxs[0] + 1, 4:6]
# add second node of last element
self.node_u_surf[
- :, -1, :] = self.mode_shape[:, idxs[0], 6:]
+ :, :, -1, :] = self.mode_shape[:, :, idxs[0], 6:]
# obtain web displacements
self.node_u_web = np.zeros(
- (self.nweb, self.nmode, self.nnode_web_max, 2))
+ (self.ncase, self.nweb, self.nmode, self.nnode_web_max, 2))
for w in range(self.nweb):
self.node_u_web[
- w, :, :self.nnode_web[w] - 1, :] = self.mode_shape[:, idxs[w] + 1:idxs[w + 1] + 1, 4:6]
+ :, w, :, :self.nnode_web[w] - 1, :] = self.mode_shape[:, :, idxs[w] + 1:idxs[w + 1] + 1, 4:6]
# add second node of last element
- self.node_u_web[w, :, self.nnode_web[w] - 1, :] = self.mode_shape[
- :, idxs[w + 1], 6:]
+ self.node_u_web[:, w, :, self.nnode_web[w] - 1, :] = self.mode_shape[
+ :, :, idxs[w + 1], 6:]
plot = False
if plot:
@@ -199,5 +389,3 @@ class FINSTRIPWrapperBlade(FINSTRIPWrapper):
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], '--')
-
- plt.show()
diff --git a/finstrip_wrapper/core/main.py b/finstrip_wrapper/core/main.py
index ee632ae..7275aec 100644
--- a/finstrip_wrapper/core/main.py
+++ b/finstrip_wrapper/core/main.py
@@ -34,23 +34,24 @@ class FINSTRIPWrapper(object):
def __init__(self, **kwargs):
super(FINSTRIPWrapper, self).__init__()
- # finstrip program parameters
- # maximum size of generated element in mm
- self.maxlength_element = 0.252E+03
- self.lambda_steps = 20 # number of half wavelength steps
- self.lambda_min = 0.01E+03 # minimum half wavelength in mm
- self.lambda_max = 15.126E+03 # maximum half wavelength in mm
+ # element length as fraction of arc length of airfoil surface
+ self.elsize_factor = 0.015
- self.spanpos = 0 # position along span of analysed cross-section
+ # finstrip program parameters
+ self.lambda_min = 0.01 # minimum half wavelength in m
+ self.lambda_max = 15.126 # maximum half wavelength in m
- self.nmode = self.lambda_steps + 1
- self.ncoord = 79
- self.ncase = 1
+ self.nmode = 21
+ #self.ncoord = 79
+ #self.ncase = 1
# flags
self.debug_mode = True
self.dry_run = False
+ def _init_half_wavelength_steps(self):
+ self.lambda_steps = self.nmode - 1 # number of half wavelength steps
+
def _init_finstrip_input_files_names(self):
''' Initializes standard input file names or extension
'''
@@ -59,10 +60,11 @@ class FINSTRIPWrapper(object):
# input file names
self.buc_sfx = '.buc'
- self.lcd_sfx = '.lcd'
+ self.lcd_sfx = '.ldc'
# self.fdb_prfx = 'db' # why no "." ??
self.txt_sfx = '.txt'
self.file_input = self.name + self.buc_sfx
+ self.file_load = self.name + self.lcd_sfx
#self.file_output = self.name + self.fcdb_prfx
#self.file_mat = 'MAT' + self.fd_sfx
#self.file_failmat = 'FAILMAT' + self.fd_sfx
diff --git a/finstrip_wrapper/omdao/finstrip_blade.py b/finstrip_wrapper/omdao/finstrip_blade.py
index 43005f4..ce2c380 100644
--- a/finstrip_wrapper/omdao/finstrip_blade.py
+++ b/finstrip_wrapper/omdao/finstrip_blade.py
@@ -1,549 +1,678 @@
import numpy as np
import os
-import shutil
+from openmdao.api import Component, Group, ParallelGroup, ExecComp
+from finstrip_wrapper.core.blade import FINSTRIPWrapperBlade
from PGL.components.airfoil import AirfoilShape
-
-from openmdao.api import Component, Group, ParallelGroup
from becas_wrapper.becas_bladestructure import set_cs_size
-from feproc_wrapper.omdao.blade import set_ncsnodes_nsets
-class FINSTRIPCrossSectionsCaseBase(Component):
+class FINSTRIPCSStructure(Component):
+ """
+ Component for computing beam structural properties
+ using the cross-sectional structure code FINSTRIP.
+
+ The code firstly calls CS2DtoFINSTRIP which is a wrapper around
+ shellexpander that comes with FINSTRIP, and second
+ calls FINSTRIP using a file interface.
+
+ parameters
+ ----------
+ config: dict
+ dictionary of model specific inputs
+ coords: array
+ cross-sectional shape. size ((ni_chord, 3))
+ matprops: array
+ material stiffness properties. Size ((10, nmat)).
+ DPs: array
+ vector of DPs. Size: (nDP)
+ r<xx><lname>T: float
+ layer thicknesses, e.g. r01triaxT.
+ r<xx><lname>A: float
+ layer angles, e.g. r01triaxA.
+ w<xx><lname>T: float
+ web thicknesses, e.g. r01triaxT.
+ w<xx><lname>A: float
+ web angles, e.g. r01triaxA.
+
+ returns
+ -------
+
+ """
- def __init__(self, name, config, caseid, sdim, ncase, nload):
- super(FINSTRIPCrossSectionsCaseBase, self).__init__()
+ def __init__(self, name, finstrip_hash, config, st3d, s, ni_chord, ncase):
+ """
+ parameters
+ ----------
+ config: dict
+ dictionary with inputs to CS2DtoFINSTRIP and FINSTRIPWrapper
+ st3d: dict
+ dictionary with blade structural definition
+ s: array
+ spanwise location of the cross-section
+ ni_chord: int
+ number of points defining the cross-section shape
+ """
+ super(FINSTRIPCSStructure, self).__init__()
+
+ self.basedir = os.getcwd()
+ self.finstrip_hash = finstrip_hash
+ self.nr = len(st3d['regions'])
+ self.ni_chord = ni_chord
+ # threshold for minimum layer thickness
+ try:
+ self.min_thickness = config['FINSTRIPWrapper']['min_thickness']
+ except:
+ self.min_thickness = 1.e-7
- self.name = name
- self.caseid = caseid
+ self.nmode = config['FINSTRIPWrapper']['nmode']
+ self.nnode_surf_max = nnode_surf_max = 150
+ self.nnode_web_max = nnode_web_max = 50
+ self.nweb = len(st3d['webs'])
self.ncase = ncase
- self.nsec = sdim[1]
- self.config = config
+ try:
+ self.use_becas_stiffness = config[
+ 'FINSTRIPWrapper']['use_becas_stiffness']
+ except:
+ self.use_becas_stiffness = False
+
+ # add materials properties array ((10, nmat))
+ self.add_param('matprops', st3d['matprops'])
+
+ # add materials strength properties array ((18, nmat))
+ # self.add_param('failmat', st3d['failmat'])
+
+ # add DPs array
+ self.add_param('stab:' + '%s:DPs' % name, np.zeros(self.nr + 1))
+ self.add_param('stab:' + '%s:rot_z_st' % name, 0.)
+
+ # add coords coords
+ self._varnames = []
+ self.add_param('stab:' + '%s:coords' % name, np.zeros((ni_chord, 3)))
+
+ self.cs2di = {}
+ self.cs2di['materials'] = st3d['materials']
+ self.cs2di['matprops'] = st3d['matprops']
+ #self.cs2di['failcrit'] = st3d['failcrit']
+ #self.cs2di['failmat'] = st3d['failmat']
+ self.cs2di['web_def'] = st3d['web_def']
+ #self.cs2di['dominant_regions'] = st3d['dominant_regions']
+ #self.cs2di['cap_DPs'] = st3d['cap_DPs']
+ #self.cs2di['le_DPs'] = st3d['le_DPs']
+ #self.cs2di['te_DPs'] = st3d['te_DPs']
+ self.cs2di['s'] = s
+ self.cs2di['DPs'] = np.zeros(self.nr + 1)
+ self.cs2di['regions'] = []
+ self.cs2di['webs'] = []
+ for ireg, reg in enumerate(st3d['regions']):
+ r = {}
+ r['layers'] = reg['layers']
+ nl = len(reg['layers'])
+ r['thicknesses'] = np.zeros(nl)
+ r['angles'] = np.zeros(nl)
+ self.cs2di['regions'].append(r)
+ for i, lname in enumerate(reg['layers']):
+ varname = '%s:r%02d%s' % (name, ireg, lname)
+ self._varnames.append(varname)
+ for ireg, reg in enumerate(st3d['webs']):
+ r = {}
+ r['layers'] = reg['layers']
+ nl = len(reg['layers'])
+ r['thicknesses'] = np.zeros(nl)
+ r['angles'] = np.zeros(nl)
+ self.cs2di['webs'].append(r)
+ for i, lname in enumerate(reg['layers']):
+ varname = '%s:w%02d%s' % (name, ireg, lname)
+ self._varnames.append(varname)
+ self.add_param(
+ 'stab:' + name + ':tvec', np.zeros(len(self._varnames) * 2))
+
+ if self.use_becas_stiffness:
+ cs_size, cs_size_ref = set_cs_size(config)
+ self.add_param('stab:' + '%s:csprops' % name, np.zeros((cs_size)))
+ self.add_param('stab:' + '%s:csprops_ref' %
+ name, np.zeros((cs_size_ref)))
+
+ self.add_param('load_cases_%s' % name, np.zeros((ncase, 6)))
+
+ # add outputs
+ self.add_output('%s:f_eig' %
+ name, np.zeros((self.ncase, self.nmode, 2)))
+ self.add_output('%s:nnode_surf' % name, 0)
+ self.add_output('%s:nnode_web' % name, np.zeros(self.nweb).astype(int))
+ self.add_output('%s:node_loc_surf' %
+ name, np.zeros((self.nnode_surf_max, 2)))
+ self.add_output('%s:node_loc_web' %
+ name, np.zeros((self.nweb, self.nnode_web_max, 2)))
+ self.add_output('%s:node_u_surf' %
+ name, np.zeros((self.ncase, self.nmode, self.nnode_surf_max, 2)))
+ self.add_output('%s:node_u_web' %
+ name, np.zeros((self.ncase, self.nweb, self.nmode, self.nnode_web_max, 2)))
+
+ #self.add_output('%s:cs_props' % name, np.zeros(cs_size))
+ #self.add_output('%s:csprops_ref' % name, np.zeros(cs_size_ref))
+ #self.add_output('%s:k_matrix' % name, shape=(6, 6))
+ #self.add_output('%s:m_matrix' % name, shape=(6, 6))
+ #self.cs_props_m1 = np.zeros(cs_size)
+ #self.csprops_ref_m1 = np.zeros(cs_size_ref)
+ #self.k_matrix_m1 = np.zeros((6, 6))
+ #self.m_matrix_m1 = np.zeros((6, 6))
+
+ #self.add_output('%s:DPcoords' % name, np.zeros((self.nr + 1, 3)))
+
+ self.workdir = 'finstrip_%s_%i' % (name, self.finstrip_hash)
+ # not so nice hack to ensure unique directory names when
+ # running parallel FD
+ self.add_output('stab:' + name + ':hash', float(self.finstrip_hash))
+
+ self.finstrip = FINSTRIPWrapperBlade(**config['FINSTRIPWrapper'])
+
+ def _params2dict(self, params):
+ """
+ convert the OpenMDAO params dictionary into
+ the dictionary format used in FINSTRIP.
+ """
+ tvec = params['stab:' + self.name + ':tvec']
+
+ self.cs2d = {}
+ # constants
+ self.cs2d['s'] = self.cs2di['s']
+ self.cs2d['web_def'] = self.cs2di['web_def']
+ #self.cs2d['failcrit'] = self.cs2di['failcrit']
+ self.cs2d['materials'] = self.cs2di['materials']
+ #self.cs2d['dominant_regions'] = self.cs2di['dominant_regions']
+ #self.cs2d['cap_DPs'] = self.cs2di['cap_DPs']
+ #self.cs2d['le_DPs'] = self.cs2di['le_DPs']
+ #self.cs2d['te_DPs'] = self.cs2di['te_DPs']
+
+ # params
+ self.cs2d['coords'] = params['stab:' + '%s:coords' % self.name][:, :2]
+ self.cs2d['matprops'] = params['matprops']
+ #self.cs2d['failmat'] = params['failmat']
+
+ self.cs2d['DPs'] = params['stab:' + '%s:DPs' % self.name]
+ self.cs2d['regions'] = []
+ self.cs2d['webs'] = []
+ counter = 0
+ nvar = len(self._varnames)
+ for ireg, reg in enumerate(self.cs2di['regions']):
+ self.cs2d['regions'].append({})
+ Ts = []
+ As = []
+ layers = []
+ for i, lname in enumerate(reg['layers']):
+ if tvec[counter] > self.min_thickness:
+ Ts.append(tvec[counter])
+ As.append(tvec[nvar + counter])
+ layers.append(lname)
+ counter += 1
+ self.cs2d['regions'][ireg]['thicknesses'] = np.asarray(Ts)
+ self.cs2d['regions'][ireg]['angles'] = np.asarray(As)
+ self.cs2d['regions'][ireg]['layers'] = layers
+ for ireg, reg in enumerate(self.cs2di['webs']):
+ self.cs2d['webs'].append({})
+ Ts = []
+ As = []
+ layers = []
+ for i, lname in enumerate(reg['layers']):
+ if tvec[counter] > self.min_thickness:
+ Ts.append(tvec[counter])
+ As.append(tvec[nvar + counter])
+ layers.append(lname)
+ counter += 1
+ self.cs2d['webs'][ireg]['thicknesses'] = np.asarray(Ts)
+ self.cs2d['webs'][ireg]['angles'] = np.asarray(As)
+ self.cs2d['webs'][ireg]['layers'] = layers
+
+ def solve_nonlinear(self, params, unknowns, resids):
+ """
+ calls FINSTRIP to compute the eigenvalues and mode shapes
+ """
+
+ try:
+ os.mkdir(self.workdir)
+ except:
+ pass
+ os.chdir(self.workdir)
+
+ self._params2dict(params)
- self.cs_size, self.cs_size_ref = set_cs_size(config)
- self.ncsnodes, _, _ = set_ncsnodes_nsets(config)
+ self.finstrip.init_dicts(self.cs2d)
- self.add_param('blade_length', 0., units='m', desc='Blade length')
- self.add_param('z_st', np.zeros(self.nsec), units='m',
- desc='non-dimensionalised y-coordinate of blade axis in structural grid')
- self.add_param('rot_z_st', np.zeros(self.nsec),
- desc='rotation about blade axis')
- self.add_param('chord_st', np.zeros(self.nsec), units='m',
- desc='blade chord distribution in structural grid')
- self.add_param('p_le_st', np.zeros(self.nsec), units='m',
- desc='blade pitch axis aft leading edge in structural grid')
- self.add_param('blade_beam_structure', np.zeros((self.nsec, self.cs_size)),
- desc='Beam properties of the blade')
- self.add_param('blade_beam_csprops_ref', np.zeros((self.nsec, self.cs_size_ref)),
- desc='Cross section properties of the blade')
+ twist = params['stab:' + '%s:rot_z_st' % self.name] # deg
+ cs_loads_rcs = params['load_cases_%s' % self.name]
+ self.finstrip.init_cs_load_cases(cs_loads_rcs, twist)
- #case_name = 'case_00'
- self.add_param('loc_cs_nodes', np.zeros((self.ncsnodes, 6, self.nsec)))
+ if self.use_becas_stiffness:
+ self.finstrip.init_becas_stiffness(params['stab:' + '%s:csprops' %
+ self.name],
+ params['stab:' + '%s:csprops_ref' %
+ self.name])
- self.add_param('n_cs_nodes', np.zeros((self.nsec)).astype(int))
+ self.finstrip.compute()
- self.add_output('%s:cs_strains' % name, np.zeros((self.ncsnodes, 4, self.nsec)),
- desc='Cross section strains of the blade')
- self.add_output('%s:loc_cs_nodes_local' % name, np.zeros((self.ncsnodes, 6, self.nsec)),
- desc='Local coordinates of cs nodes')
- self.add_output('%s:blade_beam_csprops_ref_global' % name, np.zeros((self.nsec, self.cs_size_ref)),
- desc='Global coordinates of cs nodes')
+ self.unknowns['%s:f_eig' % self.name] = self.finstrip.f_eig
+ nnode_surf = self.finstrip.nnode_surf
+ nnode_web_max = np.max(self.finstrip.nnode_web)
+ nweb_actual = self.finstrip.nweb
+ self.unknowns['%s:nnode_surf' % self.name] = nnode_surf
+ self.unknowns['%s:nnode_web' % self.name] = self.finstrip.nnode_web
+ self.unknowns['%s:node_loc_surf' % self.name][
+ :nnode_surf, :] = self.finstrip.node_loc_surf
+ self.unknowns['%s:node_loc_web' % self.name][
+ :nweb_actual, :nnode_web_max, :] = self.finstrip.node_loc_web
+ self.unknowns['%s:node_u_surf' % self.name][
+ :, :, :nnode_surf, :] = self.finstrip.node_u_surf
+ self.unknowns['%s:node_u_web' % self.name][
+ :, :nweb_actual, :, :nnode_web_max, :] = self.finstrip.node_u_web
+ os.chdir(self.basedir)
+
+
+class FINSTRIPCSStructureLoader(FINSTRIPCSStructure):
+
+ def __init__(self, name, finstrip_hash, config, st3d, s, ni_chord, cs_size, cs_size_ref, pickled_unkn,
+ grp_name=''):
+ super(FINSTRIPCSStructureLoader, self).__init__(
+ name, finstrip_hash, config, st3d, s, ni_chord, cs_size, cs_size_ref)
+
+ self.pickled_unkn = pickled_unkn
+ self.grp_name = grp_name
+ self.promoted = []
+
+ def solve_nonlinear(self, params, unknowns, resids):
+ # write data base pickled dict content into unknowns
+ for kp, vp in self.pickled_unkn.iteritems():
+ for ku in unknowns.iterkeys():
+ if kp.endswith(ku) and kp.startswith(self.grp_name):
+ unknowns[ku] = vp
+ elif kp.endswith(ku) and self.promoted:
+ unknowns[ku] = vp
-class FINSTRIPBladeCrossSectionsCase(FINSTRIPCrossSectionsCaseBase):
- """
- Component for analytical calculations of a blade cross section.
- Returns cross section angles and strains due to a full-scale blade load case.
+class Slice(Component):
"""
+ simple component for slicing arrays into vectors
+ for passing to sub-comps computing the eig values
+
+ parameters
+ ----------
+ DP<xx>: array
+ arrays of DPs along span. Size: (nsec)
+ blade_surface_norm_st: array
+ blade surface with structural discretization, no twist and prebend.
+ Size: ((ni_chord, nsec, 3))
+
+ returns
+ -------
+ sec<xxx>DPs: array
+ Vector of DPs along chord for each section. Size (nDP)
+ sec<xxx>coords: array
+ Array of cross section coords shapes. Size ((ni_chord, 3))
+ """
+
+ def __init__(self, st3d, sdim, config):
+ """
+ parameters
+ ----------
+ DPs: array
+ DPs array, size: ((nsec, nDP))
+ sdim: array
+ blade surface. Size: ((ni_chord, nsec, 3))
+ """
+ super(Slice, self).__init__()
- def __init__(self, name, config, caseid, sdim, ncase, nload):
- super(FINSTRIPBladeCrossSectionsCase, self).__init__(
- name, config, caseid, sdim, ncase, nload)
+ self.nsec = sdim[1]
+ DPs = st3d['DPs']
+ self.nDP = DPs.shape[1]
try:
- self.load_case_moment_distro = config[
- 'FINSTRIP']['load_case_moment_distro']
+ self.use_becas_stiffness = config[
+ 'FINSTRIPWrapper']['use_becas_stiffness']
except:
- self.load_case_moment_distro = False
+ self.use_becas_stiffness = False
+
+ for i in range(self.nDP):
+ self.add_param('DP%02d' % i, DPs[:, i])
+
+ self.add_param('blade_surface_norm_st', np.zeros(sdim))
+ #self.add_param('blade_surface_st', np.zeros(sdim))
+ self.add_param('blade_length', 0.)
+ self.add_param('rot_z_st', np.zeros(self.nsec))
+
+ if self.use_becas_stiffness:
+ cs_size, cs_size_ref = set_cs_size(config)
+ self.add_param('blade_beam_structure', np.zeros((self.nsec, cs_size)),
+ desc='Beam properties of the blade')
+ self.add_param('blade_beam_csprops_ref', np.zeros((self.nsec, cs_size_ref)),
+ desc='Cross section properties of the blade')
+
+ vsize = 0
+ self._varnames = []
+ for ireg, reg in enumerate(st3d['regions']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'r%02d%s' % (ireg, lname)
+ self.add_param(varname + 'T', np.zeros(self.nsec))
+ self.add_param(varname + 'A', np.zeros(self.nsec))
+ self._varnames.append(varname)
+ for ireg, reg in enumerate(st3d['webs']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'w%02d%s' % (ireg, lname)
+ self.add_param(varname + 'T', np.zeros(self.nsec))
+ self.add_param(varname + 'A', np.zeros(self.nsec))
+ self._varnames.append(varname)
+
+ for i in range(self.nsec):
+ self.add_output('stab:sec%03d:DPs' % i, DPs[i, :])
+ self.add_output('stab:sec%03d:coords' %
+ i, np.zeros((sdim[0], sdim[2])))
+ self.add_output('stab:sec%03d:tvec' %
+ i, np.zeros(len(self._varnames) * 2))
+ self.add_output('stab:sec%03d:rot_z_st' % i, 0.)
+
+ if self.use_becas_stiffness:
+ self.add_output('stab:sec%03d:csprops' %
+ i, np.zeros((cs_size)))
+ self.add_output('stab:sec%03d:csprops_ref' %
+ i, np.zeros((cs_size_ref)))
- self.add_param('load_cases', np.zeros((ncase, nload, 15)))
-
- self.add_output('%s:cs_angles' % name, np.zeros(self.nsec),
- desc='Cross section angles of the blade')
+ def solve_nonlinear(self, params, unknowns, resids):
- self.add_output('%s:reactions' % name, np.zeros((self.nsec, 6)),
- desc='Cross section reactions')
+ nvar = len(self._varnames)
+ for i in range(self.nsec):
+ #DPs = np.zeros(self.nDP)
+ # for j in range(self.nDP):
+ # DPs[j] = params['DP%02d' % j][i]
+
+ # convert DPs from s11 to s01 notation
+ af = AirfoilShape()
+ af.initialize(points=params['blade_surface_norm_st'][:, i, :] *
+ params['blade_length'])
+ DPs01 = np.zeros(self.nDP)
+ DPs11 = np.zeros(self.nDP)
+ for j in range(self.nDP):
+ DPs11[j] = params['DP%02d' % j][i]
+ if j == 0:
+ # fix to omit very small negative values (e.g.-1.0E-16)
+ DPs01[j] = 0.0
+ else:
+ DPs01[j] = af.s_to_01(DPs11[j])
+
+ unknowns['stab:sec%03d:DPs' % i] = DPs01
+ unknowns['stab:sec%03d:rot_z_st' % i] = params['rot_z_st'][i]
+ unknowns['stab:sec%03d:coords' % i] = params['blade_surface_norm_st'][:, i, :] * \
+ params['blade_length']
+ # unknowns['sec%03d:coords' % i] = params['blade_surface_st'][:, i, :] * \
+ # params['blade_length']
+
+ for ii, name in enumerate(self._varnames):
+ unknowns['stab:sec%03d:tvec' % i][ii] = params[name + 'T'][i]
+ unknowns['stab:sec%03d:tvec' % i][
+ nvar + ii] = params[name + 'A'][i]
+
+ if self.use_becas_stiffness:
+ unknowns['stab:sec%03d:csprops' % i] = params[
+ 'blade_beam_structure'][i, :]
+ unknowns['stab:sec%03d:csprops_ref' % i] = params[
+ 'blade_beam_csprops_ref'][i, :]
+
+
+class SliceLoader(Slice):
+
+ def __init__(self, st3d, sdim, pickled_unkn, grp_name=''):
+ super(SliceLoader, self).__init__(st3d, sdim)
+
+ self.pickled_unkn = pickled_unkn
+ self.grp_name = grp_name
+ self.promoted = []
def solve_nonlinear(self, params, unknowns, resids):
+ # write data base pickled dict content into unknowns
+ for kp, vp in self.pickled_unkn.iteritems():
+ for ku in unknowns.iterkeys():
+ if kp.endswith(ku) and kp.startswith(self.grp_name):
+ unknowns[ku] = vp
+ elif kp.endswith(ku) and self.promoted:
+ unknowns[ku] = vp
- nnodes = params['n_cs_nodes']
-
- # determine forces and moments
- # x,y,z, rotx,roty,rotz, Fx,Fy,Fz, Mx,My,Mz
- z_positions = params['load_cases'][self.caseid, :, 2]
- if self.load_case_moment_distro:
- moments = params['load_cases'][self.caseid, :, 9:12]
- else:
- forces = params['load_cases'][self.caseid, :, 6:9]
- #force_y = params['load_cases'][self.ncase, :, 7]
- #force_z = params['load_cases'][self.ncase, :, 8]
-
- for nsec in range(self.nsec):
- # center rotated along principal axis
- AE = params['blade_beam_structure'][nsec, 20]
- EI_x_prime = params['blade_beam_structure'][nsec, 24]
- EI_y_prime = params['blade_beam_structure'][nsec, 27]
- elastic_x = params['blade_beam_csprops_ref'][nsec, 2]
- elastic_y = params['blade_beam_csprops_ref'][nsec, 3]
- nu = - np.degrees(
- params['blade_beam_csprops_ref'][nsec, 17]) # deg (clock-wise positive)
- twist = params['rot_z_st'][nsec] # deg
- chord = params['chord_st'][nsec] * \
- params['blade_length'] # tip end
- le_x = params['p_le_st'][nsec] * chord
- te_x = -chord + le_x
- te_y = 0
-
- if self.load_case_moment_distro:
- Mx = moments[nsec, 0]
- My = moments[nsec, 1]
-
- else:
- # transform point loads to moment distro
- moment_position = params['z_st'][nsec] * params['blade_length']
- moments = convert_point_forces_to_moment_distro(
- moment_position, z_positions, forces)
-
- Mx = moments[0]
- My = moments[1]
-
- unknowns['%s:reactions' % self.name][nsec, 3] = Mx
- unknowns['%s:reactions' % self.name][nsec, 4] = My
-
- #My = -1.0E+06
- #Mx = 0.0
-
- locx = params['loc_cs_nodes'][:nnodes[nsec], 0, nsec]
- locy = params['loc_cs_nodes'][:nnodes[nsec], 1, nsec]
-
- # create component
- blade_an = CrossSection(AE,
- EI_x_prime,
- EI_y_prime,
- elastic_x,
- elastic_y,
- nu,
- twist,
- te_x,
- te_y)
- blade_an.apply_bending_moment(Mx, My)
-
- # neutral axis angle
- gamma = blade_an.gamma
- unknowns['%s:cs_angles' % self.name][nsec] = gamma
-
- # transform elastic center to rcs
- elastic_x_rcs_, elastic_y_rcs_ = blade_an.cs2root(
- [elastic_x], [elastic_y])
- unknowns['%s:blade_beam_csprops_ref_global' %
- self.name][nsec, 2] = elastic_x_rcs_
- unknowns['%s:blade_beam_csprops_ref_global' %
- self.name][nsec, 3] = elastic_y_rcs_
-
- # coordinates transformed to cross section cs
- loccsx, loccsy = blade_an.root2cs(locx, locy)
- unknowns['%s:loc_cs_nodes_local' %
- self.name][:nnodes[nsec], 0, nsec] = loccsx
- unknowns['%s:loc_cs_nodes_local' %
- self.name][:nnodes[nsec], 1, nsec] = loccsy
-
- # strains
- epto_blade_an, epto_blade_an_ax, epto_blade_an_rotx, epto_blade_an_roty = blade_an.recover_strain(
- locx, locy)
- unknowns['%s:cs_strains' % self.name][
- :nnodes[nsec], 0, nsec] = epto_blade_an
- unknowns['%s:cs_strains' % self.name][
- :nnodes[nsec], 1, nsec] = epto_blade_an_ax
- unknowns['%s:cs_strains' % self.name][
- :nnodes[nsec], 2, nsec] = epto_blade_an_rotx
- unknowns['%s:cs_strains' % self.name][
- :nnodes[nsec], 3, nsec] = epto_blade_an_roty
-
-
-class FINSTRIPComponentCrossSectionsCase(FINSTRIPCrossSectionsCaseBase):
- """
- Component for analytical calculations of a blade component cross section.
- Returns cross section angles and strains due to a full-scale blade load case.
+class PostprocessCS(Component):
"""
+ component for gathering cross section eigen values and wave lengths
+ into array as function of span
+
+ parameters
+ ----------
+ f_eig<xxx>: array
+ array of cross section eigen values. size (nmode,2)
+
+ returns
+ -------
+ blade_f_eig: array
+ array of buckling resistances. Size ((nsec,nmode,2)).
+ """
+
+ def __init__(self, nsec, nmode, ncase):
+ """
+ parameters
+ ----------
+ nsec: int
+ number of blade sections.
+ """
+ super(PostprocessCS, self).__init__()
+
+ self.nsec = nsec
- def __init__(self, name, config, caseid, sdim, ncase, nload, nsec_b):
+ for i in range(nsec):
+ self.add_param('f_eig%03d' % i, np.zeros((ncase, nmode, 2)))
- super(FINSTRIPComponentCrossSectionsCase, self).__init__(
- name, config, caseid, sdim, ncase, nload)
+ self.add_output('blade_f_eig', np.zeros((ncase, nsec, nmode, 2)))
+
+ def solve_nonlinear(self, params, unknowns, resids):
+ """
+ aggregate results
+ """
+ for i in range(self.nsec):
+ cname = 'f_eig%03d' % i
+ cs = params[cname]
+ unknowns['blade_f_eig'][:, i, :, :] = cs
- self.nsec_b = nsec_b
- # assume same initial number of nodes for blade and component
- ncsnodes_b = self.ncsnodes
- self.add_param('%s:fullblade_z_st' % name, np.zeros(self.nsec_b), units='m',
- desc='non-dimensionalised y-coordinate of blade axis in structural grid')
- self.add_param('%s:fullblade_cs_angles' % name, np.zeros(nsec_b),
- desc='Cross section angles of the blade')
- self.add_param('%s:fullblade_beam_csprops_ref' % name, np.zeros((nsec_b, self.cs_size_ref)),
- desc='Cross section properties of the blade')
- self.add_param('%s:fullblade_cs_strains' % name, np.zeros((ncsnodes_b, 4, self.nsec_b)),
- desc='Cross section strains of the blade')
- self.add_param('%s:fullblade_loc_cs_nodes' %
- name, np.zeros((ncsnodes_b, 6, self.nsec_b)))
- self.add_param('%s:fullblade_n_cs_nodes' %
- name, np.zeros((self.nsec_b)).astype(int))
+class PostprocessCSLoader(PostprocessCS):
- # x,y,z, rotx,roty,rotz, Fx,Fy,Fz, Mx,My,Mz
- self.add_output('%s:comp_load_case' % name, np.zeros((2, 12)),
- desc='Load case for component')
+ def __init__(self, nsec, cs_size, cs_size_ref, pickled_unkn, grp_name='', promoted=[]):
+ super(PostprocessCSLoader, self).__init__(nsec, cs_size, cs_size_ref)
+
+ self.pickled_unkn = pickled_unkn
+ self.grp_name = grp_name
+ self.promoted = promoted
def solve_nonlinear(self, params, unknowns, resids):
+ # write data base pickled dict content into unknowns
+ for kp, vp in self.pickled_unkn.iteritems():
+ for ku in unknowns.iterkeys():
+ if kp.endswith(ku) and kp.startswith(self.grp_name):
+ unknowns[ku] = vp
+ elif kp.endswith(ku) and self.promoted:
+ unknowns[ku] = vp
+
- nnodes_c = params['n_cs_nodes']
- nnodes_b = params['%s:fullblade_n_cs_nodes' % self.name]
-
- # component z coordinate
- comp_z = params['z_st'] * params['blade_length']
-
- # find nearest index in blade z
- blade_z = params['%s:fullblade_z_st' %
- self.name] * params['blade_length']
-
- # collect fullscale strains at TE
- node_idx = -1 # node idx on SS (TESS) at mid te?
- ref_blade_epto_an_te = []
- ref_te_x_rcs = []
- ref_te_y_rcs = []
- for nsec_bi in range(self.nsec_b):
- ref_blade_epto_an_ = params['%s:fullblade_cs_strains' % self.name][
- :nnodes_b[nsec_bi], 0, nsec_bi][node_idx]
- ref_te_x_ = params['%s:fullblade_loc_cs_nodes' % self.name][
- :nnodes_b[nsec_bi], 0, nsec_bi][node_idx]
- ref_te_y_ = params['%s:fullblade_loc_cs_nodes' % self.name][
- :nnodes_b[nsec_bi], 1, nsec_bi][node_idx]
-
- ref_blade_epto_an_te.append(ref_blade_epto_an_)
- ref_te_x_rcs.append(ref_te_x_)
- ref_te_y_rcs.append(ref_te_y_)
- ref_blade_epto_an_te = np.array(ref_blade_epto_an_te)
- ref_te_x_rcs = np.array(ref_te_x_rcs)
- ref_te_y_rcs = np.array(ref_te_y_rcs)
-
- blade_an_elastic_x = params[
- '%s:fullblade_beam_csprops_ref' % self.name][:, 2]
- blade_an_elastic_y = params[
- '%s:fullblade_beam_csprops_ref' % self.name][:, 3]
-
- blade_an_gamma = params['%s:fullblade_cs_angles' % self.name][:]
-
- # interpolate fullscale strains at TE
- ref_blade_epto_an_te_comp_z = np.interp(
- comp_z, blade_z, ref_blade_epto_an_te)
- ref_te_x_comp_z_rcs = np.interp(comp_z, blade_z, ref_te_x_rcs)
- ref_te_y_comp_z_rcs = np.interp(comp_z, blade_z, ref_te_y_rcs)
- blade_an_elastic_x_comp_z = np.interp(
- comp_z, blade_z, blade_an_elastic_x)
- blade_an_elastic_y_comp_z = np.interp(
- comp_z, blade_z, blade_an_elastic_y)
- blade_an_gamma_comp_z = np.interp(
- comp_z, blade_z, blade_an_gamma)
-
- # get component data
- AE_c = params['blade_beam_structure'][:, 20] # 9.9E+15
- EI_x_prime_c = params['blade_beam_structure'][:, 24]
- EI_y_prime_c = params['blade_beam_structure'][:, 27]
- elastic_x_c = params['blade_beam_csprops_ref'][:, 2]
- elastic_y_c = params['blade_beam_csprops_ref'][:, 3]
- nu_c = - np.degrees(
- params['blade_beam_csprops_ref'][:, 17]) # deg (clock-wise positive)
- twist_c = params['rot_z_st'][:] # deg
-
- # strain at the TE reference point at trailing edge
-
- beta = []
- l_fx = []
- l_fy = []
- Fz = []
- f_x = []
- f_y = []
- f_x_rcs = []
- f_y_rcs = []
- epto_an_te = []
- elastic_x_c_rcs = []
- elastic_y_c_rcs = []
- blade_an_elastic_x_comp_z_rcs = []
- blade_an_elastic_y_comp_z_rcs = []
- m_x = []
- m_y = []
-
- for nsec_c in range(self.nsec): # [nsec_cm]: #
-
- chord = params['chord_st'][nsec_c] * \
- params['blade_length'] # tip end
- le_x = params['p_le_st'][nsec_c] * chord
- te_x_c = -chord + le_x
- te_y_c = 0
-
- # TODO: remove dependency from trailing edge
- comp_an = CrossSection(AE_c[nsec_c],
- EI_x_prime_c[nsec_c],
- EI_y_prime_c[nsec_c],
- elastic_x_c[nsec_c],
- elastic_y_c[nsec_c],
- nu_c[nsec_c],
- twist_c[nsec_c],
- te_x_c, te_y_c)
-
- # convert elastic center from cs to root
- elastic_x_c_rcs_, elastic_y_c_rcs_ = comp_an.cs2root(
- [elastic_x_c[nsec_c]], [elastic_y_c[nsec_c]])
- elastic_x_c_rcs.append(elastic_x_c_rcs_[0])
- elastic_y_c_rcs.append(elastic_y_c_rcs_[0])
-
- unknowns['%s:blade_beam_csprops_ref_global' %
- self.name][nsec_c, 2] = elastic_x_c_rcs_
- unknowns['%s:blade_beam_csprops_ref_global' %
- self.name][nsec_c, 3] = elastic_y_c_rcs_
-
- # get coordinates from a "dummy" feproc mesh run
- # root coordinates
- locx_ = params['loc_cs_nodes'][:nnodes_c[nsec_c], 0, nsec_c]
- locy_ = params['loc_cs_nodes'][:nnodes_c[nsec_c], 1, nsec_c]
- # cs coordinates
- loccsx_, loccsy_ = comp_an.root2cs(locx_, locy_)
- #loccsx_, loccsy_ = comp_an.cs2root(locx_, locy_)
- unknowns['%s:loc_cs_nodes_local' %
- self.name][:nnodes_c[nsec_c], 0, nsec_c] = loccsx_
- unknowns['%s:loc_cs_nodes_local' %
- self.name][:nnodes_c[nsec_c], 1, nsec_c] = loccsy_
-
- # get blade gamma and elastic center
- beta_ = comp_an.required_load_angle(blade_an_gamma_comp_z[nsec_c])
- beta.append(beta_)
-
- # determine vector components of load intro point (from centroid in
- # cs cs)
- l_fx_, l_fy_ = comp_an.required_arm2(
- beta_, blade_an_elastic_x_comp_z[nsec_c], blade_an_elastic_y_comp_z[nsec_c])
- l_fx.append(l_fx_)
- l_fy.append(l_fy_)
-
- # transform blade elastic center from cs to root
- blade_an_elastic_x_comp_z_rcs_, blade_an_elastic_y_comp_z_rcs_ = comp_an.cs2root(
- [blade_an_elastic_x_comp_z[nsec_c]], [blade_an_elastic_y_comp_z[nsec_c]])
- blade_an_elastic_x_comp_z_rcs.append(
- blade_an_elastic_x_comp_z_rcs_[0])
- blade_an_elastic_y_comp_z_rcs.append(
- blade_an_elastic_y_comp_z_rcs_[0])
-
- # transform full cs coordinates from root to cs
- ref_te_x_comp_z_, ref_te_y_comp_z_ = comp_an.root2cs(
- [ref_te_x_comp_z_rcs[nsec_c]], [ref_te_y_comp_z_rcs[nsec_c]])
-
- # determine required load (if blade was prismatic)
- Fz_ = comp_an.required_load(
- l_fx_, beta_, ref_te_x_comp_z_[0], ref_te_y_comp_z_[0],
- ref_blade_epto_an_te_comp_z[nsec_c])
- Fz.append(Fz_)
-
- # determine required load introduction point in cs cs
- f_x_, f_y_ = comp_an.load_intro_intro_point2(beta_, l_fx_)
- f_x.append(f_x_)
- f_y.append(f_y_)
-
- # convert load intro point from cs to root
- f_x_rcs_, f_y_rcs_ = comp_an.cs2root(
- [f_x_], [f_y_])
- f_x_rcs.append(f_x_rcs_[0])
- f_y_rcs.append(f_y_rcs_[0])
-
- # cantilever arm for bending moment
- a_x = f_x_rcs_ - elastic_x_c_rcs_
- a_y = f_y_rcs_ - elastic_y_c_rcs_
- m_y_ = Fz_ * -a_x
- m_x_ = Fz_ * -a_y
- m_x.append(float(m_x_))
- m_y.append(float(m_y_))
-
- # apply required load
- comp_an.apply_axial_load2(Fz_, f_x_, f_y_)
-
- # recover strains of applied load
- _epto_an, _epto_an_ax, _epto_an_rotx, _epto_an_roty = comp_an.recover_strain(
- loccsx_, loccsy_, rcs=False)
-
- # get strains at node idx
- epto_an_te_ = _epto_an[node_idx]
- epto_an_te.append(epto_an_te_)
-
- # return strains of cross section
- unknowns['%s:cs_strains' % self.name][
- :nnodes_c[nsec_c], 0, nsec_c] = _epto_an
- unknowns['%s:cs_strains' % self.name][
- :nnodes_c[nsec_c], 1, nsec_c] = _epto_an_ax
- unknowns['%s:cs_strains' % self.name][
- :nnodes_c[nsec_c], 2, nsec_c] = _epto_an_rotx
- unknowns['%s:cs_strains' % self.name][
- :nnodes_c[nsec_c], 3, nsec_c] = _epto_an_roty
-
- # choose applied bending moment at mid section
- mean_z = comp_z[0] + 0.5 * (comp_z[-1] - comp_z[0])
- Fz_mean = np.interp(mean_z, comp_z, Fz)
- #Fz_applied = np.mean(Fz[4:6])
-
- m_x = np.array(m_x)
- m_y = np.array(m_y)
-
- # determine load intro positions with constant force
- a_y_applied = - m_x / Fz_mean
- a_x_applied = - m_y / Fz_mean
-
- # rcs = root coordinate system
- f_x_rcs_applied = a_x_applied + elastic_x_c_rcs
- f_y_rcs_applied = a_y_applied + elastic_y_c_rcs
-
- # fit load intro positions along specimen
- # cs_idx = range(self.nsec)
- # alternatively use two neighbouring plus mid cs to fit load intro
- # positions by linear extrapolation
- fit_mid = False
- if fit_mid:
- mid_idx = self.nsec / 2
- cs_idx = [mid_idx - 1, mid_idx, mid_idx + 1]
- else:
- cs_idx = range(len(comp_z))
-
- f_x_rcs_fit = np.polyfit(comp_z[cs_idx], f_x_rcs_applied[cs_idx], 1)
- f_x_rcs_fit_fn = np.poly1d(f_x_rcs_fit)
-
- f_y_rcs_fit = np.polyfit(comp_z[cs_idx], f_y_rcs_applied[cs_idx], 1)
- f_y_rcs_fit_fn = np.poly1d(f_y_rcs_fit)
-
- # TODO: remove zhinge hardcoding
- zhinge = 0.1
- f_z_tip = comp_z[-1] + zhinge
- f_z_root = comp_z[0] - zhinge
-
- f_x_tip = f_x_rcs_fit_fn(f_z_tip)
- f_y_tip = f_y_rcs_fit_fn(f_z_tip)
-
- # determine required load intro at root end
- f_x_root = f_x_rcs_fit_fn(f_z_root)
- f_y_root = f_y_rcs_fit_fn(f_z_root)
-
- # applied load (double the full-scale test load)
- # TODO: remove factor 2 hardcoding
- Fz_applied = 2 * Fz_mean
-
- # angles
- #unknowns['%s:cs_angles' % self.name][0, :] = beta
-
- import matplotlib.pylab as plt
- plt.figure()
- plt.plot(comp_z, beta)
- plt.plot(comp_z, nu_c)
- plt.plot(comp_z, twist_c)
- plt.plot(comp_z, beta + twist_c)
-
- plt.figure()
- plt.plot(comp_z, blade_an_elastic_x_comp_z_rcs)
- plt.plot(comp_z, elastic_x_c_rcs)
- plt.plot(comp_z, f_x_rcs)
- plt.plot(comp_z, f_x_rcs_fit_fn(comp_z), 'k-o')
- plt.plot(comp_z, f_x_rcs_applied)
-
- plt.figure()
- plt.plot(comp_z, blade_an_elastic_y_comp_z_rcs)
- plt.plot(comp_z, elastic_y_c_rcs)
- plt.plot(comp_z, f_y_rcs)
- plt.plot(comp_z, f_y_rcs_fit_fn(comp_z), 'k-o')
- plt.plot(comp_z, f_y_rcs_applied)
-
- plt.figure()
- plt.plot(blade_an_elastic_x_comp_z_rcs,
- blade_an_elastic_y_comp_z_rcs, 'bo')
- plt.plot(
- blade_an_elastic_x_comp_z_rcs[0], blade_an_elastic_y_comp_z_rcs[0], 'bs')
- plt.plot(elastic_x_c_rcs, elastic_y_c_rcs, 'go')
- plt.plot(elastic_x_c_rcs[0], elastic_y_c_rcs[0], 'gs')
- plt.plot(f_x_rcs, f_y_rcs, 'ro')
- plt.plot(f_x_rcs[0], f_y_rcs[0], 'rs')
- plt.plot(ref_te_x_comp_z_rcs, ref_te_y_comp_z_rcs, 'mo')
- plt.plot(ref_te_x_comp_z_rcs[0], ref_te_y_comp_z_rcs[0], 'ms')
- plt.plot(blade_an_elastic_x_comp_z[
- 0], blade_an_elastic_y_comp_z[0], 'b+')
- plt.plot(elastic_x_c, elastic_y_c, 'g+')
- plt.plot(elastic_x_c[0], elastic_y_c[0], 'gs')
- plt.plot(f_x, f_y, 'r+')
- plt.plot(f_x[0], f_y[0], 'rs')
- plt.plot(ref_te_x_comp_z_, ref_te_y_comp_z_, 'm+')
- plt.plot(ref_te_x_comp_z_[0], ref_te_y_comp_z_[0], 'ms')
- plt.axis('equal')
-
- plt.figure()
- plt.plot(comp_z, ref_blade_epto_an_te_comp_z)
- plt.plot(comp_z, epto_an_te)
-
- plt.figure()
- plt.plot(ref_te_x_comp_z_rcs[-1], ref_te_y_comp_z_rcs[-1], 'ro')
- plt.plot(ref_te_x_comp_z_, ref_te_y_comp_z_, 'bo')
- plt.plot(locx_, locy_, 'r')
- plt.plot(loccsx_, loccsy_, 'b')
-
- plt.figure()
- for nsec_c in range(self.nsec):
- plt.plot(params['loc_cs_nodes'][:nnodes_c[nsec_c], 0, nsec_c],
- params['loc_cs_nodes'][:nnodes_c[nsec_c], 1, nsec_c])
- plt.axis('equal')
-
- plt.figure()
- plt.plot(comp_z, m_x)
- plt.plot(comp_z, m_y)
- plt.plot(comp_z, Fz)
- # plt.show()
-
- unknowns['%s:comp_load_case' % self.name][0, 8] = Fz_applied
- unknowns['%s:comp_load_case' % self.name][0, 0] = f_x_tip
- unknowns['%s:comp_load_case' % self.name][0, 1] = f_y_tip
- unknowns['%s:comp_load_case' % self.name][0, 2] = f_z_tip
- unknowns['%s:comp_load_case' % self.name][1, 0] = f_x_root
- unknowns['%s:comp_load_case' % self.name][1, 1] = f_y_root
- unknowns['%s:comp_load_case' % self.name][1, 2] = f_z_root
-
-
-class FINSTRIPCrossSectionsBase(Group):
-
- def __init__(self):
- super(FINSTRIPCrossSectionsBase, self).__init__()
- #self.nsec = sdim[1]
- #self.ncase = ncase
-
- self.par = self.add('par', ParallelGroup(), promotes=['*'])
-
- def add_full_scale_analysis_cases(self, config, sdim, ncase, nload):
- for caseid in range(ncase):
- self.par.add('case_%02d' % caseid, FINSTRIPBladeCrossSectionsCase(
- 'case_%02d' % caseid, config, caseid, sdim, ncase, nload), promotes=['*'])
-
- def add_component_analysis_cases(self, config, sdim, ncase, nload, nsec_b):
- for caseid in range(ncase):
- self.par.add('case_%02d' % caseid, FINSTRIPComponentCrossSectionsCase(
- 'case_%02d' % caseid, config, caseid, sdim, ncase, nload, nsec_b), promotes=['*'])
-
-
-class FINSTRIPBladeCrossSections(FINSTRIPCrossSectionsBase):
+class FINSTRIPBeamStructure(Group):
+ """
+ Group for computing blade buckling resistance using FINSTRIP.
+
+ The geometric and structural inputs used are defined
+ in detail in FUSED-Wind.
+
+ parameters
+ ----------
+ blade_surface_norm_st: array
+ blade surface with structural discretization, no twist and prebend.
+ Size: ((ni_chord, nsec, 3))
+ matprops: array
+ material stiffness properties. Size (10, nmat).
+ sec<xx>DPs: array
+ 2D array of DPs. Size: ((nsec, nDP))
+ sec<xx>coords: array
+ blade surface. Size: ((ni_chord, nsec, 3))
+ sec<xx>r<yy><lname>T: array
+ region layer thicknesses, e.g. r01triaxT. Size (nsec)
+ sec<xx>r<yy><lname>A: array
+ region layer angles, e.g. r01triaxA. Size (nsec)
+ sec<xx>w<yy><lname>T: array
+ web layer thicknesses, e.g. r01triaxT. Size (nsec)
+ sec<xx>w<yy><lname>A: array
+ web layer angles, e.g. r01triaxA. Size (nsec)
+
+ returns
+ -------
+ blade_f_eig: array
+ array of buckling resistances. Size ((nsec,nmode,2)).
"""
- Group for analysing cross sections of different load cases in parallel.
+
+ def __init__(self, group, config, st3d, sdim, ncase):
+ """
+ initializes parameters and adds a csprops component
+ for each section
+
+ parameters
+ ----------
+ config: dict
+ dictionary of inputs for the cs_code class
+ st3d: dict
+ dictionary of blade structure properties
+ surface: array
+ blade surface with structural discretization.
+ Size: ((ni_chord, nsec, 3))
+ """
+ super(FINSTRIPBeamStructure, self).__init__()
+
+ # check that the config is ok
+ if not 'FINSTRIPWrapper' in config.keys():
+ raise RuntimeError('You need to supply a config dict',
+ 'for FINSTRIPWrapper')
+
+ self.st3d = st3d
+ nsec = st3d['s'].shape[0]
+ try:
+ nmode = config['FINSTRIPWrapper']['nmode']
+ except:
+ nmode = 10
+
+ # create a unique ID for this group so that FD's are not overwritten
+ self.add('hash_c', ExecComp('finstrip_hash=%f' %
+ float(self.__hash__())), promotes=['*'])
+
+ # add comp to slice the 2D arrays DPs and surface
+ promotions = []
+ self.add('slice', Slice(st3d, sdim, config), promotes=['*'])
+
+ self._varnames = []
+ for ireg, reg in enumerate(st3d['regions']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'r%02d%s' % (ireg, lname)
+ self._varnames.append(varname)
+ for ireg, reg in enumerate(st3d['webs']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'w%02d%s' % (ireg, lname)
+ self._varnames.append(varname)
+
+ # # now add a component for each section
+ par = self.add('par', ParallelGroup(), promotes=['*'])
+
+ for i in range(nsec):
+ secname = 'sec%03d' % i
+ par.add(secname, FINSTRIPCSStructure(secname, self.__hash__(), config, st3d,
+ st3d['s'][i], sdim[0], ncase), promotes=['*'])
+
+ promotions = ['blade_f_eig']
+ self.add('postpro', PostprocessCS(
+ nsec, nmode, ncase), promotes=promotions)
+ for i in range(nsec):
+ secname = 'sec%03d' % i
+ self.connect('%s:f_eig' % secname, 'postpro.f_eig%03d' % i)
+
+
+class FINSTRIPBeamStructureLoader(Group):
"""
+ Component for loading eigen values stored in a database by a
+ previous run.
+
+ This component can be replaced by FINSTRIPBeamStructure group in case the FINSTRIP
+ run shall be skipped, i.e. when the structure or load case does not change during
+ iterations.
- def __init__(self, config, sdim, ncase, nload):
- super(FINSTRIPBladeCrossSections, self).__init__()
+ Note: If parameters are promoted on the level higher than the group
+ FINSTRIPBeamStructure belongs to, they need to be listed in the promoted list.
+ """
- self.add_full_scale_analysis_cases(config, sdim, ncase, nload)
+ def __init__(self, group, config, st3d, sdim, pickled_unkn, grp_name=''):
+ """
+ initializes parameters
+
+ parameters
+ ----------
+ config: dict
+ dictionary of inputs for the cs_code class
+ st3d: dict
+ dictionary of blade structure properties
+ pickled_unkn: dict
+ unknowns dict read in via SqliteDict()
+ grp_name: str
+ name of the upper group
+ promoted: list
+ list of promoted variables in group above upper group
+
+ Note:
+ load pickled unkowns from data base
+ >>> import sqlitedict
+ >>> db = sqlitedict.SqliteDict('some.db', 'iterations')
+ >>> data = db['rank0:Driver|1']
+ >>> pickled_unkn = data['Unknowns']
+
+ """
+ super(FINSTRIPBeamStructureLoader, self).__init__()
+
+ # check that the config is ok
+ if not 'FINSTRIPWrapper' in config.keys():
+ raise RuntimeError('You need to supply a config dict',
+ 'for FINSTRIPWrapper')
+
+ self.st3d = st3d
+ nsec = st3d['s'].shape[0]
+ try:
+ nmode = config['FINSTRIPWrapper']['nmode']
+ except:
+ nmode = 10
+
+ # create a unique ID for this group so that FD's are not overwritten
+ self.add('hash_c', ExecComp('finstrip_hash=%f' %
+ float(self.__hash__())), promotes=['*'])
+
+ # add comp to slice the 2D arrays DPs and surface
+ self.add('slice', SliceLoader(
+ st3d, sdim, pickled_unkn, grp_name), promotes=['*'])
+
+ self._varnames = []
+ for ireg, reg in enumerate(st3d['regions']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'r%02d%s' % (ireg, lname)
+ self._varnames.append(varname)
+ for ireg, reg in enumerate(st3d['webs']):
+ for i, lname in enumerate(reg['layers']):
+ varname = 'w%02d%s' % (ireg, lname)
+ self._varnames.append(varname)
+
+ # # now add a component for each section
+ par = self.add('par', ParallelGroup(), promotes=['*'])
+
+ for i in range(nsec):
+ secname = 'sec%03d' % i
+ par.add(secname, FINSTRIPCSStructureLoader(secname, self.__hash__(), config, st3d,
+ st3d['s'][i], sdim[0], pickled_unkn, grp_name), promotes=['*'])
+
+ promotions = ['blade_f_eig']
+ self.add('postpro', PostprocessCS(nsec, nmode), promotes=promotions)
+ for i in range(nsec):
+ secname = 'sec%03d' % i
+ self.connect('%s:f_eig' % secname, 'postpro.f_eig%03d' % i)
diff --git a/finstrip_wrapper/omdao/test/test_finstrip_blade.py b/finstrip_wrapper/omdao/test/test_finstrip_blade.py
index 6cd0155..f01b045 100644
--- a/finstrip_wrapper/omdao/test/test_finstrip_blade.py
+++ b/finstrip_wrapper/omdao/test/test_finstrip_blade.py
@@ -23,6 +23,8 @@ from becas_wrapper.becas_bladestructure import BECASBeamStructureKM,\
from becas_wrapper.becas_stressrecovery import BECASStressRecovery
from distutils.spawn import find_executable
+from finstrip_wrapper.omdao.finstrip_blade import FINSTRIPBeamStructure
+from robe.lib.load import convert_point_forces_to_moment_distro
_matlab_installed = find_executable('matlab')
@@ -43,12 +45,12 @@ else:
DRYRUN = False
-def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=False, with_sr=False, param2=False):
+def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=False, with_sr=False, param2=False, with_becas=False):
p = Problem(impl=impl, root=Group())
p.root.add('blade_length_c', IndepVarComp(
- 'blade_length', 86.366), promotes=['*'])
+ 'blade_length', 86.366, units='m'), promotes=['*'])
pf = read_blade_planform(
os.path.join(path_data, 'DTU_10MW_RWT_blade_axis_prebend.dat'))
@@ -116,30 +118,79 @@ def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=
cfg['debug_mode'] = False
config['BECASWrapper'] = cfg
- p.root.add('stiffness', BECASBeamStructure(p.root, config, st3dn,
- (200, nsec_st, 3)), promotes=['*'])
- if with_sr:
- p.root.add('stress_recovery', BECASStressRecovery(
- config, s_st, 2), promotes=['*'])
+ if with_becas:
+ p.root.add('stiffness', BECASBeamStructure(p.root, config, st3dn,
+ (200, nsec_st, 3)), promotes=['*'])
+ # if with_sr:
+ # p.root.add('stress_recovery', BECASStressRecovery(
+ # config, s_st, 2), promotes=['*'])
+
+ cfg = {}
+ cfg['nmode'] = 20
+ cfg['core_mats'] = ['balsa']
+ cfg['elsize_factor'] = 0.015
+ cfg['lambda_min'] = 0.01 # minimum half wavelength in m
+ cfg['lambda_max'] = 15.126 # maximum half wavelength in m
+ cfg['use_becas_stiffness'] = with_becas
+ config['FINSTRIPWrapper'] = cfg
+
+ ncase = 2
+ p.root.add('stability', FINSTRIPBeamStructure(p.root, config, st3dn,
+ (200, nsec_st, 3), ncase), promotes=['*'])
p.setup()
- p['hub_radius'] = 2.8
for k, v in pf.iteritems():
if k in p.root.pf_splines.params.keys():
p.root.pf_splines.params[k] = v
- if with_sr:
- # set some arbitrary values in the load vectors used to compute strains
- for i, x in enumerate(s_st):
- try:
- p['load_cases_sec%03d' % i] = np.ones((2, 6))
- except:
- pass
-
+ # assuming we have an external shear force loading along the blade, which
+ # need to be transformed to internal section load case
+ nload = 9
+ cfg = {}
+ cfg['load_cases'] = np.zeros((ncase, nload, 15))
+ hr = 2.8
+ cfg['load_cases'][0, :, :] = np.array(
+ [[0.26658496, -0.064059925, 20.073 - hr, 0.0, 0.0, 0.0,
+ 0.0, 230000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.49336469, -0.04793707, 30.437 - hr, 0.0, 0.0,
+ 0.0, 0.0, 270000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.47505332, -0.053173575, 33.028 - hr, 0.0,
+ 0.0, 0.0, 290000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.37001989, -0.01194097, 47.71 - hr, 0.0, 0.0, 0.0,
+ 180000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.33257294, 0.00242491, 52.892 - hr, 0.0, 0.0, 0.0,
+ 0.0, 250000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.26347361, 0.01553866, 62.393 - hr, 0.0, 0.0, 0.0,
+ 130000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.23826286, 0.01686549, 65.847 - hr, 0.0, 0.0, 0.0,
+ 0.0, 220000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.16251535, 0.017204025, 76.211 - hr, 0.0, 0.0, 0.0,
+ 18000, 190000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ [0.09939912, 0.014232575, 84.848 - hr, 0.0, 0.0, 0.0,
+ 25000, 165000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
+ csmoments = np.zeros((len(s_st), 3))
+ moment_positions = s_st * p['blade_length']
+ for cs, _ in enumerate(csmoments):
+ csmoments[cs, :] = convert_point_forces_to_moment_distro(
+ moment_position=moment_positions[cs],
+ force_positions=cfg['load_cases'][0, :, 2],
+ forces=cfg['load_cases'][0, :, 6:9])
+ for i, _ in enumerate(s_st):
+ p['load_cases_sec%03d' % i] = np.zeros((ncase, 6))
+ # Fx,Fy,Fz,Mx,My,Mz
+ p['load_cases_sec%03d' % i][0, 3] = csmoments[i, 0]
+ p['load_cases_sec%03d' % i][0, 4] = csmoments[i, 1]
+
+ # make second dummy load case
+ p['load_cases_sec%03d' % i][1, 3] = -2 * csmoments[i, 0]
+ p['load_cases_sec%03d' % i][1, 4] = -2 * csmoments[i, 1]
+
+ if with_becas:
+ p['hub_radius'] = hr
return p
-class BECASBeamStructureTestCase(unittest.TestCase):
+class FINSTRIPBeamStructureTestCase(unittest.TestCase):
def setUp(self):
pass
@@ -147,6 +198,7 @@ class BECASBeamStructureTestCase(unittest.TestCase):
def tearDown(self):
pass
+ #@unittest.skip("skipping")
@unittest.skipIf(not _matlab_installed,
"Matlab not available on this system")
def test_standard_matlab(self):
@@ -155,8 +207,50 @@ class BECASBeamStructureTestCase(unittest.TestCase):
dry_run=DRYRUN,
FPM=False,
with_sr=True,
- param2=False)
+ param2=False,
+ with_becas=True)
+ p.run()
+
+ ncase = 2
+ nsec = len(p['s_st'])
+ f_eig_min = np.zeros((ncase, nsec))
+ for case in range(ncase):
+ for sec in range(nsec):
+ f_eig_min[case, sec] = np.min(
+ p['blade_f_eig'][case, sec, :, 1])
+
+ f_eig_min_exp = np.array([[0.727, 0.833, 1.385, 1.],
+ [4.734, 3.066, 1.694, 1.]])
+ for case in range(ncase):
+ for i, f in enumerate(f_eig_min[case, :]):
+ self.assertAlmostEqual(
+ f, f_eig_min_exp[case, i], places=3)
+
+ #@unittest.skip("skipping")
+ def test_finstrip(self):
+
+ p = configure_BECASBeamStructure(4, 'matlab', 'data',
+ dry_run=DRYRUN,
+ FPM=False,
+ with_sr=True,
+ param2=False,
+ with_becas=False)
+
p.run()
+ ncase = 2
+ nsec = len(p['s_st'])
+ f_eig_min = np.zeros((ncase, nsec))
+ for case in range(ncase):
+ for sec in range(nsec):
+ f_eig_min[case, sec] = np.min(
+ p['blade_f_eig'][case, sec, :, 1])
+
+ f_eig_min_exp = np.array([[0.996, 1.877, 2.09, 1.],
+ [6.484, 5.243, 2.011, 1.]])
+ for case in range(ncase):
+ for i, f in enumerate(f_eig_min[case, :]):
+ self.assertAlmostEqual(
+ f, f_eig_min_exp[case, i], places=3)
if __name__ == "__main__":
--
GitLab