diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 7140104a..5d62337c 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -18,7 +18,7 @@ jobs: fail-fast: False matrix: os: [ubuntu-latest, macos-latest] - python-version: ["3.9", "3.10", "3.11"] + python-version: ["3.13", "3.14"] steps: - name: checkout repository @@ -36,10 +36,10 @@ jobs: activate-environment: test auto-activate-base: false - - name: Download and install anba + - name: Download and install b3_secfem run: | - git clone https://github.com/ANBA4/anba4.git - cd anba4 + git clone git@github.com:gbarter/b3_secfem.git + cd b3_secfem pip install -e . - name: Install SONATA diff --git a/README.md b/README.md index 8ab92f1d..e7d40782 100644 --- a/README.md +++ b/README.md @@ -36,11 +36,11 @@ conda activate sonata-env cd .. ``` -Next, in the same conda environment compile ANBA4 (open-source) +Next, in the same conda environment compile b3_secfem4 (open-source) ``` -git clone git@github.com:ANBA4/anba4.git # (or git clone https://github.com/ANBA4/anba4.git) -cd anba4 +git clone git@github.com:gbarter/b3_secfem4.git # (or git clone https://github.com/gbarter/b3_secfem4.git) +cd b3_secfem4 pip install -e . cd .. ``` @@ -100,4 +100,4 @@ Detailed instructions can be found [here](examples/2_IEA22MW/README.md). Further **Pflumm, T., Rex, W., Hajek, M.:** Propagation of Material and Manufacturing Uncertainties in Composite Helicopter Rotor Blades, 45th European Rotorcraft Forum, Warsaw, Poland, 2019 [[more…\]](https://mediatum.ub.tum.de/1520025) [BibTeX\]](https://mediatum.ub.tum.de/export/1520025/bibtex) -**Porter, J. H., Mace, T., Bortolotti, P., et al.:** Prediction of structural damping in a composite structure from coupon tests, 2025, Preprint. [[link]](https://doi.org/10.2139/ssrn.5408061) +**Porter, J. H., Mace, T., Bortolotti, P., et al.:** Novel application of the modal strain energy technique for state-of-the-art damping predictions, 2026, Composite Structures. [[link]](https://doi.org/10.1016/j.compstruct.2026.120394) diff --git a/SONATA/anbax/anbax_utl.py b/SONATA/b3_secfem/b3_secfem_utl.py similarity index 86% rename from SONATA/anbax/anbax_utl.py rename to SONATA/b3_secfem/b3_secfem_utl.py index 543f6eea..863d1b41 100644 --- a/SONATA/anbax/anbax_utl.py +++ b/SONATA/b3_secfem/b3_secfem_utl.py @@ -10,7 +10,7 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): """function to generate the dolfin.Mesh from a SONATA-CBM definition to run - with anbax + with b3_secfem Parameters ---------- @@ -24,7 +24,7 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): Returns ---------- mesh : dolfin.Mesh - matLibrary : vector of anbax materials + matLibrary : vector of b3_secfem materials materials : dolfin.MeshFunction definign cell materials plane_orientations : dolfin.MeshFunction defining cell plane orientations fiber_orientations : dolfin.MeshFunction defining cell material fiber orientation @@ -33,7 +33,7 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): Notes ---------- the cells of cbm_mesh already contain the nodes. So the information is - currently passed twice. But consistent with the export_cells_for_vabs. + currently passed twice. """ # Would like to avoid writing mesh to a file in the future, but for now just give dummy name @@ -47,7 +47,7 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): coords = np.zeros((len(cbm_nodes), 2)) for n in cbm_nodes: coords[n.id-1,:] = n.coordinates - + n_cells = len(cbm_mesh) tris = np.zeros((n_cells, 3), dtype=np.int32) fiber_orientations = np.zeros(n_cells) @@ -69,7 +69,7 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): materials_vec = [materials_vec[m] for m in oci] fiber_orientations = fiber_orientations[oci] plane_orientations = plane_orientations[oci] - + # This meshtag doesn't do anything if not using the "region_mat" approach, but code complains if no tags are given indices = np.arange(n_cells, dtype=np.int32) values = np.array([c.id-1 for c in cbm_mesh], dtype=np.int32) @@ -94,32 +94,32 @@ def dolfin_solve(cbm_mesh, cbm_nodes, cbm_materials): # Keep oci tracker consistent result.oci = oci[result.oci] - + return result -def anbax_unit_recovery(anba, T=None): +def b3_secfem_unit_recovery(b3_secfem_in, T=None): """ Function to recover unit stresses and strains from an applied loading Results are generated in global coordinates with dimensions: ([3 unit F + 3 unit M, num elements, 6 voigt]) INPUTS: - anba - dolfin construct from anbax - T - Transformation matrix to convert results from ANBA to SONATA/VABS coordinates + b3_secfem - dolfin construct from b3_secfem + T - Transformation matrix to convert results from b3_secfem to SONATA coordinates OUTPUTS: - elem_stress_tran - global stress field in SONATA/VABS coordinates - elem_strain_tran - global strain field in SONATA/VABS coordinates + elem_stress_tran - global stress field in SONATA coordinates + elem_strain_tran - global strain field in SONATA coordinates """ - fields = b3_secfem.recover_unit_load_strains(anba) #unit_load_ + fields = b3_secfem.recover_unit_load_strains(b3_secfem_in) #unit_load_ stress = fields.sigma.copy() # dim: ([3 unit F + 3 unitM, num elements, 6 voigt]) stressM = fields.sigma_mat.copy() # dim: ([3 unit F + 3 unitM, num elements, 6 voigt]) strain = fields.epsilon.copy() # dim: ([3 unit F + 3 unitM, num elements, 6 voigt]) if T is None: T = np.eye(3) - + # Initialize the outputs: elem_stress_tran = np.zeros(stress.shape) elem_stressM_tran = np.zeros(stress.shape) @@ -149,7 +149,7 @@ def anbax_unit_recovery(anba, T=None): [strain[ii,k,4], strain[ii,k,3], strain[ii,k,2]], ]) - # Rotate the matrix to SONATA/VABS coordinates + # Rotate the matrix to SONATA coordinates istress = T.T @ elem_stress_mat[ ii,k,:,:] @ T istressM = T.T @ elem_stressM_mat[ii,k,:,:] @ T istrain = T.T @ elem_strain_mat[ ii,k,:,:] @ T @@ -158,32 +158,32 @@ def anbax_unit_recovery(anba, T=None): elem_stress_tran[ ii,k,:] = np.r_[np.diag(istress ), istress[ 1,2], istress[ 0,2], istress[ 0,1]] elem_stressM_tran[ii,k,:] = np.r_[np.diag(istressM), istressM[1,2], istressM[0,2], istressM[0,1]] elem_strain_tran[ ii,k,:] = np.r_[np.diag(istrain ), istrain[ 1,2], istrain[ 0,2], istrain[ 0,1]] - - oci2orig = np.argsort(anba.oci) + + oci2orig = np.argsort(b3_secfem_in.oci) elem_stress_tran = elem_stress_tran[ :,oci2orig,:] elem_stressM_tran = elem_stressM_tran[:,oci2orig,:] elem_strain_tran = elem_strain_tran[ :,oci2orig,:] - + return elem_stress_tran, elem_stressM_tran, elem_strain_tran -def anbax_recovery(anba, force, moment, T=None): +def b3_secfem_recovery(b3_secfem, force, moment, T=None): """ Function to recover total stresses and strains from an applied loading Results are generated in global coordinates with dimensions: ([num elements, 6 voigt]) INPUTS: - anba - dolfin construct from anbax - force - Forces in anbax coordinates, [F1, F2, F3], e.g. force = [2.2, 3.4, 1.1] - moment - Moments in anbax coordinates, [M1, M2, M3], e.g. moment = [4.2, 5.7, 6.2] - T - Transformation matrix to convert results from ANBA to SONATA/VABS coordinates + b3_secfem - dolfin construct from b3_secfem + force - Forces in b3_secfem coordinates, [F1, F2, F3], e.g. force = [2.2, 3.4, 1.1] + moment - Moments in b3_secfem coordinates, [M1, M2, M3], e.g. moment = [4.2, 5.7, 6.2] + T - Transformation matrix to convert results from b3_secfem to SONATA coordinates OUTPUTS: - stress_sum - global total stress field in SONATA/VABS coordinates - strain_sum - global total strain field in SONATA/VABS coordinates + stress_sum - global total stress field in SONATA coordinates + strain_sum - global total strain field in SONATA coordinates """ - - stress, stressM, strain = anbax_unit_recovery(anba, T=T) + + stress, stressM, strain = b3_secfem_unit_recovery(b3_secfem, T=T) n_el = stress.shape[1] # Rotate forces and moments @@ -192,7 +192,7 @@ def anbax_recovery(anba, force, moment, T=None): else: myforce = T.T @ force @ T mymoment = T.T @ moment @ T - + # Add up total stress and strain through super-position linear combo of unit forces and moments stress_sum = np.zeros((n_el,6)) stressM_sum = np.zeros((n_el,6)) @@ -200,10 +200,10 @@ def anbax_recovery(anba, force, moment, T=None): for k in range(3): stress_sum += myforce[ k] * stress[ k,:,:] stress_sum += mymoment[k] * stress[3+k,:,:] - + stressM_sum += myforce[ k] * stressM[ k,:,:] stressM_sum += mymoment[k] * stressM[3+k,:,:] - + strain_sum += myforce[ k] * strain[ k,:,:] strain_sum += mymoment[k] * strain[3+k,:,:] diff --git a/SONATA/anbax/classANBAXConfig.py b/SONATA/b3_secfem/classb3_secfemConfig.py similarity index 82% rename from SONATA/anbax/classANBAXConfig.py rename to SONATA/b3_secfem/classb3_secfemConfig.py index 6bf235b4..2c5f6690 100644 --- a/SONATA/anbax/classANBAXConfig.py +++ b/SONATA/b3_secfem/classb3_secfemConfig.py @@ -7,10 +7,10 @@ import numpy as np -class ANBAXConfig(object): +class b3_secfemConfig(object): """ - this class contains the Configuration for a anbax + this class contains the Configuration for a b3_secfem Attributes: ---------- @@ -29,7 +29,7 @@ class ANBAXConfig(object): def __init__(self, **kw): self.recover_flag = 0 self.ref_sys = "global" - self.voigt_convention = "anba" + self.voigt_convention = "b3_secfem" if "recover_flag" in kw: self.recover_flag = kw["recover_flag"] @@ -40,4 +40,4 @@ def __init__(self, **kw): if __name__ == "__main__": - test = ANBAXConfig(recover_flag=1) + test = b3_secfemConfig(recover_flag=1) diff --git a/SONATA/cbm/classBeamSectionalProps.py b/SONATA/cbm/classBeamSectionalProps.py index 64fba765..bad6abe7 100644 --- a/SONATA/cbm/classBeamSectionalProps.py +++ b/SONATA/cbm/classBeamSectionalProps.py @@ -20,14 +20,7 @@ class BeamSectionalProps(object): """ - this class stores the beam cross-sectional data and has methods to - read the result files from vabs - - The coordinate system follows the definition of VABS and DYMORE. - For more information see the the vabs and dymore user manual. - the x/x1 axis is in the direction of the beam. the y/x2 axis points towards - the leading edge, and z/x3 points accordingly upwards. - + this class stores the beam cross-sectional data The 6x6 sectional stiffness matrix, TS (Timoshenko Stiffness Matrix) (1-extension; 2,3-shear, 4-twist; 5,6-bending) relates the sectional axial strain, epsilon1, transverse shearing strains, epsilon2 and epsilon3, @@ -153,9 +146,6 @@ def __init__(self, fname=None): self.ELE = None self.U = None - if fname: - self.read_vabs_K(fname) - @property def m00(self): return self.MM[0, 0] diff --git a/SONATA/cbm/classCBM.py b/SONATA/cbm/classCBM.py index 96182eed..027e55d6 100644 --- a/SONATA/cbm/classCBM.py +++ b/SONATA/cbm/classCBM.py @@ -22,7 +22,6 @@ # First party modules from SONATA.cbm.cbm_utl import trsf_sixbysix from SONATA.cbm.classBeamSectionalProps import BeamSectionalProps -from SONATA.cbm.classCBMConfig import CBMConfig from SONATA.cbm.display.display_mesh import plot_cells from SONATA.cbm.mesh.cell import Cell @@ -39,8 +38,8 @@ from SONATA.cbm.topo.utils import getID from SONATA.cbm.topo.web import Web from SONATA.cbm.topo.weight import Weight -from SONATA.vabs.classStrain import Strain -from SONATA.vabs.classStress import Stress +from SONATA.cbm.classStrain import Strain +from SONATA.cbm.classStress import Stress from SONATA.cbm.topo.projection import ( chop_interval_from_layup, sort_layup_projection,) # from SONATA.vabs.classVABSConfig import VABSConfig @@ -50,11 +49,11 @@ from OCC.Core.gp import gp_Pnt2d try: - from SONATA.anbax.anbax_utl import dolfin_solve, anbax_unit_recovery, anbax_recovery + from SONATA.b3_secfem.b3_secfem_utl import dolfin_solve, b3_secfem_unit_recovery, b3_secfem_recovery except: - print("dolfin and anbax could not be imported!") + print("dolfin and b3_secfem could not be imported!") raise @@ -99,11 +98,8 @@ class CBM(object): cbm_review_mesh() prints a summary of the mesh properties to the screen - cbm_run_vabs() - runs the solver VABS (Variational Asymptotic Beam Sectional Analysis) - - cbm_run_anbax() - runs the solver anbax from macro morandini + cbm_run_b3_secfem() + runs the solver b3_secfem cbm_post_2dmesh(attribute='MatID', title='NOTITLE', **kw) displays the mesh with specified attributes with matplotlib @@ -138,7 +134,6 @@ class CBM(object): >>> job.cbm_gen_topo() >>> job.cbm_gen_mesh() >>> job.cbm_review_mesh() - >>> job.cbm_run_vabs() >>> job.cbm_post_2dmesh(title='Hello World!') """ @@ -254,11 +249,11 @@ def check_for_bspline_intersections(self, segment): [intersected, intersection_pnt] = self.check_bspline_intersections(Boundary_BSplineLst) if intersected: print("WARNING: There is an intersection in the structure.") - plt.figure() - self.display_bsplinelst(self.SegmentLst[0].BSplineLst, 'black') - self.display_bsplinelst(Boundary_BSplineLst, 'blue') - for points in intersection_pnt: - plt.plot(points.X(), points.Y(), 'x', color = 'red', linewidth = 4, markersize = 10) + # plt.figure() + # self.display_bsplinelst(self.SegmentLst[0].BSplineLst, 'black') + # self.display_bsplinelst(Boundary_BSplineLst, 'blue') + # for points in intersection_pnt: + # plt.plot(points.X(), points.Y(), 'x', color = 'red', linewidth = 4, markersize = 10) def cbm_gen_topo(self, **kwargs): """ @@ -405,7 +400,7 @@ def cbm_gen_mesh(self, **kwargs): self.mesh.extend(bw_cells) - # invert nodes list of all cell to make sure they are counterclockwise for vabs in the right coordinate system! + # invert nodes list of all cell to make sure they are counterclockwise for c in self.mesh: if not c.orientation: c.invert_nodes() @@ -509,7 +504,7 @@ def cbm_custom_mesh(self, nodes, cells, materials, split_quads=True, reading the custom mesh. theta_3 : float, optional Value for fiber orientation angle to be passed down into SONATA - and ANBA. If None, then zero is passed down. Units are degrees. + and b3_secfem. If None, then zero is passed down. Units are degrees. The default value is None. Returns @@ -564,8 +559,8 @@ def cbm_custom_mesh(self, nodes, cells, materials, split_quads=True, return None - def cbm_run_anbax(self): - """interface method to run the solver anbax from marco.morandini + def cbm_run_b3_secfem(self): + """interface method to run the solver b3_secfem Notes ---------- @@ -589,21 +584,21 @@ def cbm_run_anbax(self): try: - anba = dolfin_solve(self.mesh, nodes, self.materials) + b3_secfem = dolfin_solve(self.mesh, nodes, self.materials) except: print('\n') print('==========================================\n\n') - print('Error, Anba4 wrapper called, likely ') - print('Anba4 _or_ Dolfin are not installed\n\n') + print('Error, b3_secfem4 wrapper called, likely ') + print('b3_secfem4 _or_ Dolfin are not installed\n\n') print('==========================================\n\n') raise - #TBD: pass it to anbax and run it! - tmp_TS = anba.K # get stiffness matrix - tmp_MM = anba.M # get mass matrix + #TBD: pass it to b3_secfem and run it! + tmp_TS = b3_secfem.K # get stiffness matrix + tmp_MM = b3_secfem.M # get mass matrix - # Define transformation T (from ANBA to SONATA/VABS coordinates) + # Define transformation T (from b3_secfem to SONATA coordinates) B = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) T = np.dot(np.identity(3), np.linalg.inv(B)) @@ -611,15 +606,15 @@ def cbm_run_anbax(self): self.BeamProperties.TS = trsf_sixbysix(tmp_TS, T) self.BeamProperties.MM = trsf_sixbysix(tmp_MM, T) - # self.BeamProperties.Xm = anba.mass_center # mass center - is already allocated from mass matrix - self.BeamProperties.Xt = anba.tension_center - self.BeamProperties.Xs = anba.shear_center + # self.BeamProperties.Xm = b3_secfem.mass_center # mass center - is already allocated from mass matrix + self.BeamProperties.Xt = b3_secfem.tension_center + self.BeamProperties.Xs = b3_secfem.shear_center # --- Stress & Strain recovery --- # - if self.config.anbax_cfg.recover_flag: - print("STATUS:\t Running ANBAX Stress & Strain Recovery:") - elem_stress, elem_stressM, elem_strain = anbax_recovery(anba, self.config.anbax_cfg.F.tolist(), self.config.anbax_cfg.M.tolist(), T) + if self.config.b3_secfem_cfg.recover_flag: + print("STATUS:\t Running b3_secfem Stress & Strain Recovery:") + elem_stress, elem_stressM, elem_strain = b3_secfem_recovery(b3_secfem, self.config.b3_secfem_cfg.F.tolist(), self.config.b3_secfem_cfg.M.tolist(), T) # ASSIGN stresses and strains to mesh elements in np.triu_indices_from order for i,c in enumerate(self.mesh): @@ -645,21 +640,21 @@ def cbm_run_viscoelastic(self, test_elastic=True, test_tau0=True): self.mesh, nodes = sort_and_reassignID(self.mesh) try: - anba = dolfin_solve(self.mesh, nodes, self.materials) + b3_secfem = dolfin_solve(self.mesh, nodes, self.materials) except: print('\n') print('==========================================\n\n') - print('Error, Anba4 wrapper called, likely ') - print('Anba4 _or_ Dolfin are not installed\n\n') + print('Error, b3_secfem4 wrapper called, likely ') + print('b3_secfem4 _or_ Dolfin are not installed\n\n') print('==========================================\n\n') raise - # Call ANBAX with baseline properties - tmp_TS = anba.K # get stiffness matrix - tmp_MM = anba.M # get mass matrix + # Call b3_secfem with baseline properties + tmp_TS = b3_secfem.K # get stiffness matrix + tmp_MM = b3_secfem.M # get mass matrix - # Define transformation T (from ANBA to SONATA/VABS coordinates) + # Define transformation T (from b3_secfem to SONATA coordinates) B = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) T = np.dot(np.identity(3), np.linalg.inv(B)) @@ -667,20 +662,20 @@ def cbm_run_viscoelastic(self, test_elastic=True, test_tau0=True): self.BeamProperties.TS = trsf_sixbysix(tmp_TS, T) self.BeamProperties.MM = trsf_sixbysix(tmp_MM, T) - # self.BeamProperties.Xm = anba.mass_center # mass center - is already allocated from mass matrix - self.BeamProperties.Xt = anba.tension_center - self.BeamProperties.Xs = anba.shear_center + # self.BeamProperties.Xm = b3_secfem.mass_center # mass center - is already allocated from mass matrix + self.BeamProperties.Xt = b3_secfem.tension_center + self.BeamProperties.Xs = b3_secfem.shear_center #print(self.BeamProperties.TS) #print(self.BeamProperties.MM) #print(self.BeamProperties.Xt) #print(self.BeamProperties.Xs) - + # Recover the mapping from sectional Forces and Moments to Strain - # in a non-invasive way from ANBA - elem_stress, _, elem_strain = anbax_unit_recovery(anba, T) + # in a non-invasive way from b3_secfem + elem_stress, _, elem_strain = b3_secfem_unit_recovery(b3_secfem, T) # Reordering is necessary from the conventional stress/strain order - # to match the constitutive tensor built from ANBAX. + # to match the constitutive tensor built from b3_secfem. # Both the strains and stresses need reordering so setting it here. # All other orders here were tested and verified to give wrong results. reorder_stress_strain = np.array([1, 2, 0, 4, 5, 3]) @@ -898,16 +893,16 @@ def cbm_exp_stress_strain_map(self, station_ind, station_pos, twist, self.mesh, nodes = sort_and_reassignID(self.mesh) try: - anba = dolfin_solve(self.mesh, nodes, self.materials) + b3_secfem = dolfin_solve(self.mesh, nodes, self.materials) except: print('\n') print('==========================================\n\n') - print('Error, Anba4 wrapper called, likely ') - print('Anba4 _or_ Dolfin are not installed\n\n') + print('Error, b3_secfem4 wrapper called, likely ') + print('b3_secfem4 _or_ Dolfin are not installed\n\n') print('==========================================\n\n') raise - # Define transformation T (from ANBA to SONATA/VABS coordinates) + # Define transformation T (from b3_secfem to SONATA coordinates) B = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) T = np.dot(np.identity(3), np.linalg.inv(B)) @@ -925,7 +920,7 @@ def cbm_exp_stress_strain_map(self, station_ind, station_pos, twist, elem_materials[i] = c.MatID elem_cxy[i, :] = c.center - # 2. Call ANBA looping over unit forces/moments + # 2. Call b3_secfem looping over unit forces/moments # Mapping should match `beam_struct_eval` external_to_internal_ind = [2, 0, 1] @@ -953,7 +948,7 @@ def cbm_exp_stress_strain_map(self, station_ind, station_pos, twist, # This ends up being potentially excessively slow since # it does a new calculation for each stress and strain field (4), # but only need the material coordinate strain field. - elem_stress, elem_stressM, elem_strain = anbax_recovery(anba, F, M, T) + elem_stress, elem_stressM, elem_strain = b3_secfem_recovery(b3_secfem, F, M, T) # 3. Store Strain results in each case / element fc_to_strain_m[:, i, :] = elem_strain.T @@ -1015,7 +1010,7 @@ def cbm_exp_stress_strain_map(self, station_ind, station_pos, twist, return - def cbm_exp_BeamDyn_beamprops(self, Theta=0, solver="vabs"): + def cbm_exp_BeamDyn_beamprops(self, Theta=0): """ Converts the Beam Properties of CBM to the correct coordinate System of BeamDyn and returns the 6x6 Stiffness matrix, the 6x6 MassMatrix. @@ -1031,7 +1026,6 @@ def cbm_exp_BeamDyn_beamprops(self, Theta=0, solver="vabs"): ---------- Theta: float, optional is the angle of rotation of the coordinate system in "radians" - solver: str, optional Returns ---------- @@ -1051,14 +1045,10 @@ def cbm_exp_BeamDyn_beamprops(self, Theta=0, solver="vabs"): """ - if solver == "vabs" or solver == "anbax": - if Theta != 0: - tmp_bp = self.BeamProperties.rotate(Theta) - else: - tmp_bp = self.BeamProperties - + if Theta != 0: + tmp_bp = self.BeamProperties.rotate(Theta) else: - print("Check solver for BeamDyn Beam Property input.") + tmp_bp = self.BeamProperties tmp_bp = copy.deepcopy(tmp_bp) @@ -1070,53 +1060,6 @@ def cbm_exp_BeamDyn_beamprops(self, Theta=0, solver="vabs"): tmp_MM = trsf_sixbysix(tmp_bp.MM, T) return (tmp_TS, tmp_MM) - - - def cbm_exp_dymore_beamprops(self, eta, Theta=0, solver="vabs", units={"mass": "kg", "length": "m", "force": "N"}): - """ - Converts the Units of CBM to DYMORE/PYMORE/MARC units and returns the - array of the beamproperties with Massterms(6), Stiffness(21), - damping(1) and curvilinear coordinate(1) - - Parameters - ---------- - - eta : float, - is the beam curvilinear coordinate of the beam from 0 to 1. - - Theta: float - is the angle of rotation of the coordinate system in "radians" - - Returns - ---------- - arr : ndarray - [Massterms(6) (m00, mEta2, mEta3, m33, m23, m22) - Stiffness(21) (k11, k12, k22, k13, k23, k33,... k16, k26, ...k66) - Viscous Damping(1) mu, Curvilinear coordinate(1) eta] - - - Notes - ---------- - - Unit Convertion takes sooo much time. Commented out for now! - - """ - if solver == "vabs" or solver == "anbax": - if Theta != 0: - tmp_bp = self.BeamProperties.rotate(Theta) - else: - tmp_bp = self.BeamProperties - - else: - print("Check solver for Dymore Beam Property input.") - - - MM = tmp_bp.MM - MASS = np.array([MM[0, 0], MM[2, 3], MM[0, 4], MM[5, 5], MM[4, 5], MM[4, 4]]) - STIFF = tmp_bp.TS[np.tril_indices(6)[1], np.tril_indices(6)[0]] - mu = 0.0 - return np.hstack((MASS, STIFF, mu, eta)) - - def cbm_post_2dmesh(self, attribute="MatID", title="NOTITLE", **kw): """ CBM Postprocessing method that displays the mesh with matplotlib. @@ -1127,13 +1070,12 @@ def cbm_post_2dmesh(self, attribute="MatID", title="NOTITLE", **kw): Uses the string to look for the cell attributes. The default attribute is MatID. Possible other attributes can be fiber orientation (theta_3) or strains and stresses. - If BeamProperties are already calculated by VABS or something - similar, elastic-axis, center-of-gravity... are displayed. + If BeamProperties are calculated, elastic-axis, center-of-gravity... are displayed. title : string, optional Title to be placed over the plot. **kw : keyword arguments, optional are passed to the lower "plot_cells" function. Such options are: - VABSProperties=None, title='None', plotTheta11=False, + BeamProperties=None, title='None', plotTheta11=False, plotDisplacement=False, savepath Returns @@ -1149,28 +1091,3 @@ def cbm_post_2dmesh(self, attribute="MatID", title="NOTITLE", **kw): mesh, nodes = sort_and_reassignID(self.mesh) fig, ax = plot_cells(self.mesh, nodes, attribute, self.materials, self.BeamProperties, title, **kw) return fig, ax - -#%%############################################################################ -# M A I N # -############################################################################### -if __name__ == "__main__": - plt.close("all") - fname = "jobs/debug/issue20/sec_config.yml" - fname = "jobs/VariSpeed/uh60a_cbm_advanced/sec_config_R2000.yml" - # fname = 'jobs/AREA/R250/sec_config.yml' - # fname = 'jobs/PBortolotti/sec_config.yml' - config = CBMConfig(fname) - - job = CBM(config) - - job.cbm_gen_topo() - job.cbm_gen_mesh(split_quads=True) - - job.cbm_review_mesh() - job.cbm_run_vabs(rm_vabfiles=False) - # AnbaBeamProperties = job.cbm_run_anbax() - - # job.cbm_post_2dmesh(title='Hello World!') - job.cbm_post_3dtopo() -# job.config.vabs_cfg.recover_flag = 1 -# job.config.vabs_cfg.M = [0,2000e4,0] diff --git a/SONATA/cbm/classCBMConfig.py b/SONATA/cbm/classCBMConfig.py index 1fd19565..963f9882 100644 --- a/SONATA/cbm/classCBMConfig.py +++ b/SONATA/cbm/classCBMConfig.py @@ -13,7 +13,7 @@ # First party modules from SONATA.classMaterial import find_material -from SONATA.anbax.classANBAXConfig import ANBAXConfig +from SONATA.b3_secfem.classb3_secfemConfig import b3_secfemConfig if __name__ == "__main__": os.chdir("/media/gu32kij/work/TPflumm/SONATA") @@ -44,13 +44,12 @@ class CBMConfig(object): flags: dict - vabs_cfg: VABSConfig - anbax_cfg: ANBAXConfig + b3_secfem_cfg: b3_secfemConfig """ - __slots__ = ("filename", "setup", "webs", "segments", "bw", "flags", "vabs_cfg", "anbax_cfg") + __slots__ = ("filename", "setup", "webs", "segments", "bw", "flags", "b3_secfem_cfg") def __init__(self, inputdata=None, materials=None): self.setup, self.webs, self.segments, self.bw = {}, {}, {}, {} @@ -62,7 +61,7 @@ def __init__(self, inputdata=None, materials=None): else: print("Input data is not a dictionary. Check yaml input file.") - self.anbax_cfg = ANBAXConfig() + self.b3_secfem_cfg = b3_secfemConfig() self.flags = {"mesh_core": True} def read_yaml_cbm(self, yml, materials): diff --git a/SONATA/vabs/classStrain.py b/SONATA/cbm/classStrain.py similarity index 96% rename from SONATA/vabs/classStrain.py rename to SONATA/cbm/classStrain.py index a7742c15..239b88d5 100644 --- a/SONATA/vabs/classStrain.py +++ b/SONATA/cbm/classStrain.py @@ -11,7 +11,7 @@ import numpy as np # First party modules -from SONATA.vabs.classTensorException import TensorException +from SONATA.cbm.classTensorException import TensorException if __name__ == "__main__": os.chdir("../..") diff --git a/SONATA/vabs/classStress.py b/SONATA/cbm/classStress.py similarity index 96% rename from SONATA/vabs/classStress.py rename to SONATA/cbm/classStress.py index d86d6746..47d48a0a 100644 --- a/SONATA/vabs/classStress.py +++ b/SONATA/cbm/classStress.py @@ -12,7 +12,7 @@ import numpy as np # First party modules -from SONATA.vabs.classTensorException import TensorException +from SONATA.cbm.classTensorException import TensorException if __name__ == "__main__": os.chdir("../..") diff --git a/SONATA/vabs/classTensorException.py b/SONATA/cbm/classTensorException.py similarity index 100% rename from SONATA/vabs/classTensorException.py rename to SONATA/cbm/classTensorException.py diff --git a/SONATA/cbm/display/display_mesh.py b/SONATA/cbm/display/display_mesh.py index 40e33061..719382ab 100755 --- a/SONATA/cbm/display/display_mesh.py +++ b/SONATA/cbm/display/display_mesh.py @@ -31,7 +31,7 @@ def centroid(points): return centroid def plot_mesh(nodes, elements, theta_11, data, data_name, materials, - title=None, VABSProperties=None, + title=None, Properties=None, show_element_number=False, show_node_number=False, invert_xaxis = True, lfactor=0.5e-2, **kw): @@ -142,27 +142,27 @@ def plot_mesh(nodes, elements, theta_11, data, data_name, materials, for i, item in enumerate(nodes): ax.annotate(i + 1, (item[0], item[1]), color="red") - if (VABSProperties is not None) and show_beam_props: + if (Properties is not None) and show_beam_props: pass - (CG,) = plt.plot(VABSProperties.Xm[0], VABSProperties.Xm[1], "ro", label="CG: Mass Center") - # ax.annotate('CG', (VABSProperties.Xm2,VABSProperties.Xm3),fontsize=20) - (NA,) = plt.plot(VABSProperties.Xt[0], VABSProperties.Xt[1], "gs", label="NA: Neutral Axes") - # ax.annotate('NA', (VABSProperties.Xt2,VABSProperties.Xt3),fontsize=20) + (CG,) = plt.plot(Properties.Xm[0], Properties.Xm[1], "ro", label="CG: Mass Center") + # ax.annotate('CG', (Properties.Xm2,Properties.Xm3),fontsize=20) + (NA,) = plt.plot(Properties.Xt[0], Properties.Xt[1], "gs", label="NA: Neutral Axes") + # ax.annotate('NA', (Properties.Xt2,Properties.Xt3),fontsize=20) GC_defined = False - if VABSProperties.Xg is not None: - (GC,) = plt.plot(VABSProperties.Xg[0], VABSProperties.Xg[1], "b^", + if Properties.Xg is not None: + (GC,) = plt.plot(Properties.Xg[0], Properties.Xg[1], "b^", label="GC: Geometric Center") - # ax.annotate('GC', (VABSProperties.Xg2,VABSProperties.Xg3),fontsize=20) + # ax.annotate('GC', (Properties.Xg2,Properties.Xg3),fontsize=20) plt.legend(handles=[CG, GC, NA]) GC_defined = True else: plt.legend(handles=[CG, NA]) - if isinstance(VABSProperties.Xs, np.ndarray): - (SC,) = plt.plot(VABSProperties.Xs[0], VABSProperties.Xs[1], "kD", + if isinstance(Properties.Xs, np.ndarray): + (SC,) = plt.plot(Properties.Xs[0], Properties.Xs[1], "kD", label="SC: Generalized Shear Center") - # ax.annotate('SC', (VABSProperties.Xs2,VABSProperties.Xs3),fontsize=20) + # ax.annotate('SC', (Properties.Xs2,Properties.Xs3),fontsize=20) if GC_defined: plt.legend(handles=[CG, GC, NA, SC]) @@ -178,7 +178,7 @@ def plot_mesh(nodes, elements, theta_11, data, data_name, materials, return (fig,ax) -def plot_cells(cells,nodes, attr1, materials, VABSProperties=None, +def plot_cells(cells,nodes, attr1, materials, Properties=None, title='None', plotTheta11=False, plotDisplacement=False, **kw): """ @@ -193,7 +193,7 @@ def plot_cells(cells,nodes, attr1, materials, VABSProperties=None, DESCRIPTION. materials : TYPE DESCRIPTION. - VABSProperties : TYPE, optional + Properties : TYPE, optional DESCRIPTION. The default is None. title : TYPE, optional DESCRIPTION. The default is 'None'. @@ -250,7 +250,7 @@ def plot_cells(cells,nodes, attr1, materials, VABSProperties=None, fig,ax = plot_mesh(nodes_array, element_array, theta_11, data, data_name, - materials, title, VABSProperties, **kw) + materials, title, Properties, **kw) if 'savepath' in kw: diff --git a/SONATA/cbm/mesh/cell.py b/SONATA/cbm/mesh/cell.py index 7b9469cf..c1b4266d 100755 --- a/SONATA/cbm/mesh/cell.py +++ b/SONATA/cbm/mesh/cell.py @@ -21,7 +21,7 @@ class Cell(object): __slots__ = ("id", "nodes", "theta_1", "theta_3", "MatID", "structured", "interior_nodes", "strain", "strainM", "stress", "stressM", "sf", "failure_mode","fm_to_strain") - # fm_to_strain is the the mapping calculated with ANBA from the + # fm_to_strain is the the mapping calculated with b3_secfem from the # Force/Moment at section to the strain (using engineering shear strain) # within the element. class_counter = 1 @@ -103,7 +103,6 @@ def calc_theta_1(self): (theta_11). Afterwards the ply coordinate system us ritated avizt y3 in the right-hand sens by the amount of Theta_3 to form the material system. - For a detailed description see docs/man/VABS-Manual.pdf Figure 4. theta_11 is calculated as the angle between the x-axis and the Vector from Node 1 to Node 2. diff --git a/SONATA/classBlade.py b/SONATA/classBlade.py index f9847693..c144f512 100644 --- a/SONATA/classBlade.py +++ b/SONATA/classBlade.py @@ -24,8 +24,8 @@ from SONATA.classComponent import Component from SONATA.classMaterial import read_materials -# SONATA anbax module -from SONATA.anbax.classANBAXConfig import ANBAXConfig +# SONATA b3_secfem module +from SONATA.b3_secfem.classb3_secfemConfig import b3_secfemConfig # SONATA cbm module from SONATA.cbm.classCBM import CBM @@ -127,7 +127,7 @@ class Blade(Component): nparray([[grid, CBM instance]],dtype = object) beam_properties : ndarray - array of grid location and VABSSectionalProp instance + array of grid location and SectionalProp instance nparray([[grid, beam_properties]],dtype = object) @@ -165,7 +165,7 @@ class Blade(Component): >>> job.read_yaml(yml.get('components').get('blade'), airfoils, materials) >>> job.blade_gen_section() - >>> job.blade_run_anbax() + >>> job.blade_run_b3_secfem() >>> job.blade_plot_sections() >>> job.blade_post_3dtopo(flag_lft = True, flag_topo = True) @@ -190,7 +190,7 @@ class Blade(Component): "f_beam_ref_axis", "f_soy", "f_curvature_k1", - "anba_beam_properties", + "b3_secfem_beam_properties", "wopwop_bsplinelst", "wopwop_pnts", "wopwop_vecs", @@ -606,7 +606,7 @@ def blade_custom_mesh(self, nodes, cells, materials, split_quads=True, and order for element. theta_3 : float, optional Value for fiber orientation angle to be passed down into SONATA - and ANBA. If None, then zero is passed down. + and b3_secfem. If None, then zero is passed down. Units are degrees. The default value is None. @@ -630,9 +630,9 @@ def blade_custom_mesh(self, nodes, cells, materials, split_quads=True, return None - def blade_run_anbax(self, loads=None, **kwargs): + def blade_run_b3_secfem(self, loads=None, **kwargs): """ - runs anbax for every section + runs b3_secfem for every section Parameters ---------- @@ -643,7 +643,7 @@ def blade_run_anbax(self, loads=None, **kwargs): """ - ac = ANBAXConfig() + ac = b3_secfemConfig() lst = [] for (x, cs) in self.sections: if loads: @@ -652,34 +652,34 @@ def blade_run_anbax(self, loads=None, **kwargs): for k,v in load.items(): setattr(ac,k,v) - cs.config.anbax_cfg = ac - print("STATUS:\t Running ANBAX at grid location %s" % (x)) - cs.cbm_run_anbax(**kwargs) + cs.config.b3_secfem_cfg = ac + print("STATUS:\t Running b3_secfem at grid location %s" % (x)) + cs.cbm_run_b3_secfem(**kwargs) lst.append([x, cs.BeamProperties]) - # self.anba_beam_properties = np.asarray(lst) + # self.b3_secfem_beam_properties = np.asarray(lst) self.beam_properties = np.asarray(lst) return None def blade_run_viscoelastic(self, **kwargs): """ - Runs anbax for every section to evaluate viscoelastic 6x6 matrices. + Runs b3_secfem for every section to evaluate viscoelastic 6x6 matrices. """ - print('Running viscoelastic analysis. This requires calling ANBAX' + print('Running viscoelastic analysis. This requires calling b3_secfem' + ' multiple times per section.') - ac = ANBAXConfig() + ac = b3_secfemConfig() lst = [] for (x, cs) in self.sections: - cs.config.anbax_cfg = ac + cs.config.b3_secfem_cfg = ac print("STATUS:\t Running Viscoelastic Analysis at grid location %s" % (x)) cs.cbm_run_viscoelastic(**kwargs) lst.append([x, cs.BeamProperties]) - # self.anba_beam_properties = np.asarray(lst) + # self.b3_secfem_beam_properties = np.asarray(lst) self.beam_properties = np.asarray(lst) return None @@ -724,10 +724,10 @@ def blade_exp_stress_strain_map(self, flag_output_zero_twist=False, """ - ac = ANBAXConfig() + ac = b3_secfemConfig() for ind,(x, cs) in enumerate(self.sections): - cs.config.anbax_cfg = ac + cs.config.b3_secfem_cfg = ac print("STATUS:\t Running Stress and Strain Maps at %s" % (x)) @@ -947,7 +947,7 @@ def blade_post_3dtopo(self, flag_wf=True, flag_lft=False, flag_topo=False, self.start_display() - def blade_exp_beam_props(self, cosy='local', style='DYMORE', eta_offset=0, solver='anbax', filename = None): + def blade_exp_beam_props(self, cosy='local', style='DYMORE', eta_offset=0): """ Exports the beam_properties in the @@ -971,10 +971,6 @@ def blade_exp_beam_props(self, cosy='local', style='DYMORE', eta_offset=0, solve coincide with the global coorinate system of the blade. The unit is in nondimensional r coordinates (x/Radius) - solver : str, optional - solver : if multiple or other solvers than vabs were applied, use - this option - filename : str, optional if the user wants to write the output to a file. @@ -991,10 +987,10 @@ def blade_exp_beam_props(self, cosy='local', style='DYMORE', eta_offset=0, solve # eta = -eta_offset/(1-eta_offset) + (1/(1-eta_offset))*cs[0] eta = (cs[0] * R) - (eta_offset * R) if style == "DYMORE": - lst.append(cs[1].cbm_exp_dymore_beamprops(eta=eta, solver=solver)) + lst.append(cs[1].cbm_exp_dymore_beamprops(eta=eta)) elif style == "BeamDyn": - lst.append(cs[1].cbm_exp_BeamDyn_beamprops(eta=eta, solver=solver)) + lst.append(cs[1].cbm_exp_BeamDyn_beamprops(eta=eta)) elif style == "CAMRADII": pass diff --git a/SONATA/classMaterial.py b/SONATA/classMaterial.py index 6c280120..3a081987 100644 --- a/SONATA/classMaterial.py +++ b/SONATA/classMaterial.py @@ -36,8 +36,7 @@ class Material: orth : int orth is the flag to indicate whether the material is isotropic (0), - orthotropic (1) or general anisotropic (2) in consitency with VABS - Manual for Users (2011) + orthotropic (1) or general anisotropic (2) rho : float density in kg/m**3 @@ -126,7 +125,7 @@ def rotated_constitutive_tensor(self, plane_orientation, fiber_orientation): (e.g., gamma_ij = 2*eps_ij = eps_ij + eps_ji) The implementation is heavily copied from: - anba4/anba4/material/material.cpp/TransformationMatrix + b3_secfem4/b3_secfem4/material/material.cpp/TransformationMatrix This is added to the python implementation for easier access when doing extra calculations for viscoelastic materials. @@ -233,7 +232,7 @@ def __init__(self, **kw): self.alpha = None self.YS = None self.UTS = None - + if kw.get("E") is not None: self.E = float(kw.get("E")) @@ -278,7 +277,7 @@ def constitutive_tensor(self): (e.g., gamma_ij = 2*eps_ij = eps_ij + eps_ji) The implementation is heavily copied from: - anba4/anba4/material/material.cpp/IsotropicMaterial + b3_secfem4/b3_secfem4/material/material.cpp/IsotropicMaterial This is added to the python implementation for easier access when doing extra calculations for viscoelastic materials. @@ -397,7 +396,7 @@ def __init__(self, flag_mat, **kw): if all(k in kw for k in ('alpha_11', 'alpha_22', 'alpha_33')): self.alpha = np.array([kw.get('alpha_11'), kw.get('alpha_22'), kw.get('alpha_33')]).astype(float) - if flag_mat: # wisdem includes vectors for the following material properties that are to be converted in order to comply with SONATA and VABS/anbax + if flag_mat: # wisdem includes vectors for the following material properties that are to be converted in order to comply with SONATA and b3_secfem if kw.get('Xt') is not None: self.Xt = float(kw.get('Xt')[0]) # retrieve axial tensile strength in [MPa] from provided 3D vector @@ -432,25 +431,26 @@ def __init__(self, flag_mat, **kw): # self.S23 = float(kw.get('S23')) self.viscoelastic = {} - - # Default b3_secfem coordinate system with fiber axis on E1 - E1, E2, E3 = self.E #[0], m.E[1], m.E[2] - G12, G13, G23 = self.G #[0], m.G[1], m.G[2] - nu12, nu13, nu23 = self.nu #[0], m.nu[1], m.nu[2] - - # If trying to trick b3_secfem to give results like Anba that has fiber axis on E3: - #E1, E2, E3 = self.E[1], self.E[2], self.E[0] - #G12, G13, G23 = self.G[2], self.G[0], self.G[1] - #nu12, nu31, nu32 = self.nu[2], self.nu[0], self.nu[1] - #nu13 = nu31*E1/E3 - #nu23 = nu32*E2/E3 - - self.b3mat = b3_secfem.OrthotropicMaterial(E1=E1, # Young's modulus, fibre [Pa] (Ezz along beam axis) - E2=E2, # Young's modulus, transverse-2 [Pa] - E3=E3, # Young's modulus, transverse-3 [Pa] - G12=G12, G13=G13, G23=G23, - nu12=nu12, nu13=nu13, nu23=nu23, - rho=self.rho, name=self.name) + + # Map SONATA/b3_secfem4 material axes to b3_secfem's (x, y, z) principal + # frame the same way classMaterial.constitutive_tensor() does for + # b3_secfem4: at zero rotation, SONATA's fiber direction (E[0]) must sit + # on the beam axis (b3 axis 3 / +z), not on b3 axis 1 (+x), which is + # what a plain "forward the card unchanged" mapping would give. + # This mirrors constitutive_tensor(): e_xx=E[1], e_yy=E[2], e_zz=E[0]. + E1_b3, E2_b3, E3_b3 = self.E[1], self.E[2], self.E[0] + G12_b3, G13_b3, G23_b3 = self.G[2], self.G[0], self.G[1] + nu12_b3, nu31_b3, nu32_b3 = self.nu[2], self.nu[0], self.nu[1] + # Convert the minor Poisson ratios (loading along z=fiber) to the + # major ratios b3_secfem expects, same relation as constitutive_tensor(). + nu13_b3 = nu31_b3 * E1_b3 / E3_b3 + nu23_b3 = nu32_b3 * E2_b3 / E3_b3 + + self.b3mat = b3_secfem.OrthotropicMaterial( + E1=E1_b3, E2=E2_b3, E3=E3_b3, + G12=G12_b3, G13=G13_b3, G23=G23_b3, + nu12=nu12_b3, nu13=nu13_b3, nu23=nu23_b3, + rho=self.rho, name=self.name) ''' # Using the abstract version that calls the b3_secfem materials class instead @@ -478,7 +478,7 @@ def constitutive_tensor(self): (e.g., gamma_ij = 2*eps_ij = eps_ij + eps_ji) The implementation is heavily copied from: - anba4/anba4/material/material.cpp/OrthotropicMaterial + b3_secfem4/b3_secfem4/material/material.cpp/OrthotropicMaterial This is added to the python implementation for easier access when doing extra calculations for viscoelastic materials. @@ -488,8 +488,8 @@ def constitutive_tensor(self): """ # using tensor direction indices z=1, x=2, y=3 - # (anba uses this direction ordering to return stress/strain) - # Implementation should match anbax_util.py > build_mat_library. + # (b3_secfem uses this direction ordering to return stress/strain) + # Implementation should match b3_secfem_util.py > build_mat_library. # # This means elastic modulus order should be: # [Along Beam, Along Perimeter, Through Thickness] diff --git a/SONATA/utl/beam_struct_eval.py b/SONATA/utl/beam_struct_eval.py index eea61f66..5409d456 100644 --- a/SONATA/utl/beam_struct_eval.py +++ b/SONATA/utl/beam_struct_eval.py @@ -24,15 +24,15 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, job_str, mu): """ - Analyse, transform, evaluate and plot structural results from VABS and/or ANBAX + Analyse, transform, evaluate and plot structural results from b3_secfem Functions: beam_struct_eval - parent function (retrieve & transform data, structure of data evaluation) plot_beam_props_6by6 - function to plot the 6x6 stiffness and mass matrices plot_beam_mass_distribution - function to plot the beam mass per unit length distribution - plot_vabs_anbax - function to plot the 6x6 stiffness and mass matrices from both VABS and ANBAX for code-to-code verification - vabs_export_beam_struct_properties - csv export of structural beam properties - anbax_export_beam_struct_properties - csv export of structural beam properties + plot_b3_secfem - function to plot the 6x6 stiffness and mass matrices from both b3_secfem for code-to-code verification + export_beam_struct_properties - csv export of structural beam properties + b3_secfem_export_beam_struct_properties - csv export of structural beam properties Inputs: flags_dict - dictionary containing relevant flags @@ -71,14 +71,14 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, flags_dict[key] = False - # --- ANBAX --- # + # --- b3_secfem --- # # --------------------------------------- # - # clear var before initializing anbax + # clear var before initializing b3_secfem job.beam_properties = None # --------------------------------------- # - # --- ANBAX --- # + # --- b3_secfem --- # if flags_dict['flag_recovery'] and not flags_dict['viscoelastic']: if np.asarray(loads_dict['Forces']).shape == (3,): @@ -99,7 +99,7 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, float(loads_dict["Moments"][2]), float(loads_dict["Moments"][0])] - loads = {# sonata coord system input converted to anbax coordinates + loads = {# sonata coord system input converted to b3_secfem coordinates "F": np.array([[0] + Forces, [1] + Forces]), "M": np.array([[0] + Moments, @@ -114,9 +114,9 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, "M" : loads_dict["Moments"][:, [0, 2, 3, 1]] } - job.blade_run_anbax(loads) # run anbax + job.blade_run_b3_secfem(loads) # run b3_secfem elif not flags_dict['viscoelastic']: - job.blade_run_anbax() # run anbax + job.blade_run_b3_secfem() # run b3_secfem else: # flags_dict['viscoelastic'] == True @@ -127,30 +127,30 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, job.blade_run_viscoelastic() # init used matrices and arrays - anbax_beam_stiff_init = np.zeros([len(cs_pos), 6, 6]) - anbax_beam_inertia_init = np.zeros([len(cs_pos), 6, 6]) - anbax_beam_stiff = np.zeros([len(cs_pos), 6, 6]) - anbax_beam_inertia = np.zeros([len(cs_pos), 6, 6]) - anbax_beam_section_mass = np.zeros([len(cs_pos), 1]) + b3_secfem_beam_stiff_init = np.zeros([len(cs_pos), 6, 6]) + b3_secfem_beam_inertia_init = np.zeros([len(cs_pos), 6, 6]) + b3_secfem_beam_stiff = np.zeros([len(cs_pos), 6, 6]) + b3_secfem_beam_inertia = np.zeros([len(cs_pos), 6, 6]) + b3_secfem_beam_section_mass = np.zeros([len(cs_pos), 1]) if flags_dict['viscoelastic']: - anbax_beam_viscoelastic \ + b3_secfem_beam_viscoelastic \ = np.zeros((len(cs_pos), len(job.beam_properties[0][1].tau), 6, 6)) # --------------------------------------- # - # retrieve & allocate ANBAX results + # retrieve & allocate b3_secfem results for i in range(len(job.beam_properties)): - anbax_beam_section_mass[i] = job.beam_properties[i, 1].m00 # receive mass per unit span + b3_secfem_beam_section_mass[i] = job.beam_properties[i, 1].m00 # receive mass per unit span for j in range(6): - anbax_beam_stiff_init[i, j, :] = np.array(job.beam_properties[i, 1].TS[j, :]) # receive 6x6 timoshenko stiffness matrix - anbax_beam_inertia_init[i, j, :] = np.array(job.beam_properties[i, 1].MM[j, :]) # receive 6x6 mass matrix + b3_secfem_beam_stiff_init[i, j, :] = np.array(job.beam_properties[i, 1].TS[j, :]) # receive 6x6 timoshenko stiffness matrix + b3_secfem_beam_inertia_init[i, j, :] = np.array(job.beam_properties[i, 1].MM[j, :]) # receive 6x6 mass matrix if flags_dict['viscoelastic']: for k in range(len(job.beam_properties[0][1].tau)): - anbax_beam_viscoelastic[i, k, :, :] = job.beam_properties[i, 1].TSv[k] + b3_secfem_beam_viscoelastic[i, k, :, :] = job.beam_properties[i, 1].TSv[k] # --------------------------------------- # - # rotate anbax results from SONATA/VABS def to BeamDyn def coordinate + # rotate b3_secfem results from SONATA/b3_secfem def to BeamDyn def coordinate # system (for flag_DeamDyn_def_transform = True) # OpenTurbine transform is from BeamDyn, so rotate in that case as well. if flags_dict['flag_DeamDyn_def_transform'] \ @@ -161,49 +161,40 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, B = np.array([[0, 0, 1], [0, -1, 0], [1, 0, 0]]) # transformation matrix T = np.dot(np.identity(3), np.linalg.inv(B)) for n_sec in range(len(cs_pos)): - anbax_beam_stiff[n_sec, :, :] = trsf_sixbysix(anbax_beam_stiff_init[n_sec, :, :], T) - anbax_beam_inertia[n_sec, :, :] = trsf_sixbysix(anbax_beam_inertia_init[n_sec, :, :], T) + b3_secfem_beam_stiff[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_stiff_init[n_sec, :, :], T) + b3_secfem_beam_inertia[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_inertia_init[n_sec, :, :], T) if flags_dict['viscoelastic']: for k in range(len(job.beam_properties[0][1].tau)): - anbax_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( - anbax_beam_viscoelastic[n_sec, k, :, :], T) + b3_secfem_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( + b3_secfem_beam_viscoelastic[n_sec, k, :, :], T) # str_ext = '_BeamDyn_def' coordsys = 'BeamDyn' - print('STATUS:\t Structural characteristics of ANBAX converted from SONATA/VABS to BeamDyn coordinate system definition!') + print('STATUS:\t Structural characteristics of b3_secfem converted from SONATA to BeamDyn coordinate system definition!') else: - anbax_beam_stiff = anbax_beam_stiff_init - anbax_beam_inertia = anbax_beam_inertia_init + b3_secfem_beam_stiff = b3_secfem_beam_stiff_init + b3_secfem_beam_inertia = b3_secfem_beam_inertia_init # str_ext = '' - coordsys = 'VABS/SONATA' + coordsys = 'SONATA' # --------------------------------------- # # Export beam structural properties to csv file if flags_dict['flag_csv_export']: - print('STATUS:\t Export csv files with structural blade characeristics from ANBAX to: ' + folder_str + 'csv_export/') - anbax_export_beam_struct_properties(folder_str, job_str, cs_pos, coordsys=coordsys, solver='anbax', beam_stiff=anbax_beam_stiff, - beam_inertia=anbax_beam_inertia, beam_mass_per_length=anbax_beam_section_mass) + print('STATUS:\t Export csv files with structural blade characeristics from b3_secfem to: ' + folder_str + 'csv_export/') + b3_secfem_export_beam_struct_properties(folder_str, job_str, cs_pos, coordsys=coordsys, beam_stiff=b3_secfem_beam_stiff, + beam_inertia=b3_secfem_beam_inertia, beam_mass_per_length=b3_secfem_beam_section_mass) - - - # --------------------------------------- # - # Export beam structural properties to csv file - # if flags_dict['flag_csv_export']: - # print('STATUS:\t Export csv files with structural blade characeristics from ANBAX to: ' + folder_str + 'csv_export/') - # export_beam_struct_properties(folder_str, job_str, cs_pos, solver='anbax', beam_stiff=anbax_beam_stiff, beam_inertia=anbax_beam_inertia, beam_mass_per_length=anbax_beam_section_mass) - - # ToDo: also export BeamDyn files for results from anbax as soon as the verification is completed # --------------------------------------- # # write BeamDyn input files - np.savetxt('anbax_BAR00.txt', np.array([cs_pos, anbax_beam_stiff[:, 3, 3], - anbax_beam_stiff[:, 4, 4], - anbax_beam_stiff[:, 5, 5], - anbax_beam_stiff[:, 2, 2], - anbax_beam_inertia[:, 0, 0]]).T) + np.savetxt('b3_secfem_BAR00.txt', np.array([cs_pos, b3_secfem_beam_stiff[:, 3, 3], + b3_secfem_beam_stiff[:, 4, 4], + b3_secfem_beam_stiff[:, 5, 5], + b3_secfem_beam_stiff[:, 2, 2], + b3_secfem_beam_inertia[:, 0, 0]]).T) if flags_dict['flag_DeamDyn_def_transform'] \ and flags_dict['flag_output_zero_twist']: @@ -231,13 +222,13 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, [-np.sin(alpha), np.cos(alpha), 0], [0, 0, 1]]) - anbax_beam_stiff[n_sec, :, :] = trsf_sixbysix(anbax_beam_stiff[n_sec, :, :], rot_mat) - anbax_beam_inertia[n_sec, :, :] = trsf_sixbysix(anbax_beam_inertia[n_sec, :, :], rot_mat) + b3_secfem_beam_stiff[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_stiff[n_sec, :, :], rot_mat) + b3_secfem_beam_inertia[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_inertia[n_sec, :, :], rot_mat) if flags_dict['viscoelastic']: for k in range(len(job.beam_properties[0][1].tau)): - anbax_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( - anbax_beam_viscoelastic[n_sec, k, :, :], rot_mat) + b3_secfem_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( + b3_secfem_beam_viscoelastic[n_sec, k, :, :], rot_mat) job.true_twist = np.copy(job.twist[:, 1]) @@ -249,7 +240,7 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, print('STATUS:\t Write BeamDyn input files') # refine = int(30/len(cs_pos)) # initiate node refinement parameter write_beamdyn_axis(folder_str, flags_dict, job_name, job.blade_ref_axis, job.twist) - write_beamdyn_prop(folder_str, flags_dict, job_name, cs_pos, anbax_beam_stiff, anbax_beam_inertia, mu) + write_beamdyn_prop(folder_str, flags_dict, job_name, cs_pos, b3_secfem_beam_stiff, b3_secfem_beam_inertia, mu) if flags_dict['viscoelastic']: @@ -257,7 +248,7 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, write_beamdyn_viscoelastic(folder_str, flags_dict, job_name, cs_pos, job.beam_properties[0][1].tau, - anbax_beam_viscoelastic) + b3_secfem_beam_viscoelastic) if flags_dict['flag_OpenTurbine_transform']: @@ -266,13 +257,13 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, T = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) for n_sec in range(len(cs_pos)): - anbax_beam_stiff[n_sec, :, :] = trsf_sixbysix(anbax_beam_stiff[n_sec, :, :], T) - anbax_beam_inertia[n_sec, :, :] = trsf_sixbysix(anbax_beam_inertia[n_sec, :, :], T) + b3_secfem_beam_stiff[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_stiff[n_sec, :, :], T) + b3_secfem_beam_inertia[n_sec, :, :] = trsf_sixbysix(b3_secfem_beam_inertia[n_sec, :, :], T) if flags_dict['viscoelastic']: for k in range(len(job.beam_properties[0][1].tau)): - anbax_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( - anbax_beam_viscoelastic[n_sec, k, :, :], T) + b3_secfem_beam_viscoelastic[n_sec, k, :, :] = trsf_sixbysix( + b3_secfem_beam_viscoelastic[n_sec, k, :, :], T) # str_ext = '_OpenTurbine_def' @@ -286,7 +277,7 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, print('STATUS:\t Write OpenTurbine input files') write_beamdyn_prop(folder_str, flags_dict, job_name, - cs_pos, anbax_beam_stiff, anbax_beam_inertia, mu, + cs_pos, b3_secfem_beam_stiff, b3_secfem_beam_inertia, mu, format_name='OpenTurbine') if flags_dict['viscoelastic']: @@ -295,7 +286,7 @@ def beam_struct_eval(job_name, flags_dict, loads_dict, cs_pos, job, folder_str, write_beamdyn_viscoelastic(folder_str, flags_dict, job_name, cs_pos, job.beam_properties[0][1].tau, - anbax_beam_viscoelastic, + b3_secfem_beam_viscoelastic, format_name='OpenTurbine') # ============================================= # @@ -324,15 +315,15 @@ def plot_beam_props_6by6(cs_pos, data, fig_title, save_path): return None -def plot_beam_axes(cs_pos, vabs_beam_mass_center, vabs_beam_neutral_axes, vabs_beam_geometric_center, - vabs_beam_shear_center, save_path): +def plot_beam_axes(cs_pos, beam_mass_center, beam_neutral_axes, beam_geometric_center, + beam_shear_center, save_path): fig = plt.figure(tight_layout=True, figsize=(8, 5), dpi=80, facecolor='w', edgecolor='k') ax = plt.axes(projection='3d') - ax.plot3D(cs_pos, vabs_beam_mass_center[:,0], vabs_beam_mass_center[:,1], label='Mass center') - ax.plot3D(cs_pos, vabs_beam_neutral_axes[:,0], vabs_beam_neutral_axes[:,1], label='Neutral axes') - ax.plot3D(cs_pos, vabs_beam_geometric_center[:,0], vabs_beam_geometric_center[:,1], label='Geometric center') - ax.plot3D(cs_pos, vabs_beam_shear_center[:,0], vabs_beam_shear_center[:,1], label='Shear Center') + ax.plot3D(cs_pos, beam_mass_center[:,0], beam_mass_center[:,1], label='Mass center') + ax.plot3D(cs_pos, beam_neutral_axes[:,0], beam_neutral_axes[:,1], label='Neutral axes') + ax.plot3D(cs_pos, beam_geometric_center[:,0], beam_geometric_center[:,1], label='Geometric center') + ax.plot3D(cs_pos, beam_shear_center[:,0], beam_shear_center[:,1], label='Shear Center') ax.set_xlabel('r/R') ax.set_ylabel('chordwise location, m') ax.set_zlabel('thickness location, m') @@ -343,11 +334,11 @@ def plot_beam_axes(cs_pos, vabs_beam_mass_center, vabs_beam_neutral_axes, vabs_b return None -def plot_beam_mass_distribution(cs_pos, vabs_beam_mass_distribution, save_path): +def plot_beam_mass_distribution(cs_pos, beam_mass_distribution, save_path): fig = plt.figure(tight_layout=True, figsize=(8, 5), dpi=80, facecolor='w', edgecolor='k') - plt.plot(cs_pos, vabs_beam_mass_distribution[:,0]) + plt.plot(cs_pos, beam_mass_distribution[:,0]) plt.xlabel('r/R') plt.ylabel('Mass per unit length, N/m') # plt.ylim([0, 1000]) @@ -361,18 +352,18 @@ def plot_beam_mass_distribution(cs_pos, vabs_beam_mass_distribution, save_path): # ============================================= # -def plot_vabs_anbax(cs_pos, vabs_data, anbax_data, fig_title, save_path): +def plot_b3_secfem(cs_pos, data, b3_secfem_data, fig_title, save_path): # plots 6x6 matrix k = 1 fig = plt.figure(tight_layout=True, figsize=(14, 10), dpi=80, facecolor='w', edgecolor='k') # fig.suptitle(fig_title) - for i in range(len(vabs_data[0, :, 0])): - for j in range(len(vabs_data[0, 0, :])): + for i in range(len(data[0, :, 0])): + for j in range(len(data[0, 0, :])): if j >= i: - ax = fig.add_subplot(len(vabs_data[0, :, 0]), len(vabs_data[0, 0, :]), k) - ax.plot(cs_pos, vabs_data[:, i, j], '--k') - ax.plot(cs_pos, anbax_data[:, i, j], ':r') - plt.ylim(1.1 * min(-1, min(vabs_data[:, i, j]), min(anbax_data[:, i, j])), 1.1 * max(1, max(vabs_data[:, i, j]), max(anbax_data[:, i, j]))) + ax = fig.add_subplot(len(data[0, :, 0]), len(data[0, 0, :]), k) + ax.plot(cs_pos, data[:, i, j], '--k') + ax.plot(cs_pos, b3_secfem_data[:, i, j], ':r') + plt.ylim(1.1 * min(-1, min(data[:, i, j]), min(b3_secfem_data[:, i, j])), 1.1 * max(1, max(data[:, i, j]), max(b3_secfem_data[:, i, j]))) ax.set_xlabel('r/R') if fig_title == 'Mass matrix': ax.set_title('$m_{%i %i}$' % ((i + 1), (j + 1))) @@ -387,84 +378,11 @@ def plot_vabs_anbax(cs_pos, vabs_data, anbax_data, fig_title, save_path): -# ============================================= # -def vabs_export_beam_struct_properties(folder_str, job_str, radial_stations, coordsys, solver, beam_stiff, beam_inertia, beam_mass_per_length, - beam_mass_center, beam_neutral_axes, beam_geometric_center, beam_shear_center): - - if solver=='vabs': - export_name_general = 'vabs_beam_properties_general.csv' - export_name_stiff = 'vabs_beam_properties_stiff_matrices.csv' - export_name_mass = 'vabs_beam_properties_mass_matrices.csv' - else: - print('Define correct solver name (vabs or anbax) when calling export_beam_struct_properties') - - # -------------------------------------------------- # - # Export mass per unit length for the defined radial stations - if not os.path.isdir(folder_str + 'csv_export/'): - os.mkdir(folder_str + 'csv_export/') - with open(''.join([folder_str + 'csv_export/' + job_str[0:-5] + '_' + export_name_general]), mode='w') as csv_file: - beam_prop_writer = csv.writer(csv_file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) - if coordsys == 'BeamDyn': - beam_prop_writer.writerow(['Coordinate system:', 'Beamdyn coordinates']) - elif coordsys == 'VABS/SONATA': - beam_prop_writer.writerow(['Coordinate system:', 'VABS/SONATA coordinates']) - else: - beam_prop_writer.writerow(['Coordinate system:', 'to be verified']) - - beam_prop_writer.writerow(['section in r/R', 'Mass per unit length [kg/m]', - 'Mass center (chordwise), m', 'Mass center (thickness), m', - 'Neutral axes (chordwise), m', 'Neutral axes (thickness), m', - 'Geometric center (chordwise), m', 'Geometric center (thickness), m', - 'Shear center (chordwise), m', 'Shear center (thickness), m']) - for i in range(len(beam_mass_per_length)): # receive number of radial sections - beam_prop_writer.writerow([str(radial_stations[i]), str(beam_mass_per_length[i,0]), - str(beam_mass_center[i,0]), str(beam_mass_center[i,1]), - str(beam_neutral_axes[i,0]), str(beam_neutral_axes[i,1]), - str(beam_geometric_center[i,0]), str(beam_geometric_center[i,1]), - str(beam_shear_center[i,0]), str(beam_shear_center[i,1])]) - - csv_file.close() - # -------------------------------------------------- # - - # Export stiffness matrices for the defined radial stations - with open(''.join([folder_str + 'csv_export/' + job_str[0:-5] + '_' + export_name_stiff]), mode='w') as csv_file: - beam_prop_writer = csv.writer(csv_file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) - - for i in range(len(beam_stiff)): # receive number of radial sections - beam_prop_writer.writerow(' ') - beam_prop_writer.writerow(['section in r/R', str(radial_stations[i])]) - - for j in range(6): # number of rows for each matrix - beam_prop_writer.writerow(beam_stiff[i, j, :]) - # beam_prop_writer.writerow(job.beam_properties[i, 1].TS[j, :]) # can eventually be called as a standalone via the job.beam_properties object - csv_file.close() - - # -------------------------------------------------- # - # Export mass matrices for the defined radial stations - with open(''.join([folder_str + 'csv_export/' + job_str[0:-5] + '_' + export_name_mass]), mode='w') as csv_file: - beam_prop_writer = csv.writer(csv_file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) - - for i in range(len(beam_inertia)): # receive number of radial sections - beam_prop_writer.writerow(' ') - beam_prop_writer.writerow(['section in r/R', str(radial_stations[i])]) - - for j in range(6): # number of rows for each matrix - beam_prop_writer.writerow(beam_inertia[i, j, :]) - # beam_prop_writer.writerow(job.beam_properties[i, 1].TS[j, :]) # can eventually be called as a standalone via the job.beam_properties object - csv_file.close() - # -------------------------------------------------- # - return None - - - -def anbax_export_beam_struct_properties(folder_str, job_str, radial_stations, coordsys, solver, beam_stiff, beam_inertia, beam_mass_per_length): +def b3_secfem_export_beam_struct_properties(folder_str, job_str, radial_stations, coordsys, beam_stiff, beam_inertia, beam_mass_per_length): - if solver=='anbax': - export_name_general = 'anbax_beam_properties_general.csv' - export_name_stiff = 'anbax_beam_properties_stiff_matrices.csv' - export_name_mass = 'anbax_beam_properties_mass_matrices.csv' - else: - print('Define correct solver name (vabs or anbax) when calling export_beam_struct_properties') + export_name_general = 'b3_secfem_beam_properties_general.csv' + export_name_stiff = 'b3_secfem_beam_properties_stiff_matrices.csv' + export_name_mass = 'b3_secfem_beam_properties_mass_matrices.csv' # -------------------------------------------------- # # Export mass per unit length for the defined radial stations @@ -474,10 +392,10 @@ def anbax_export_beam_struct_properties(folder_str, job_str, radial_stations, co beam_prop_writer = csv.writer(csv_file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) if coordsys == 'BeamDyn': beam_prop_writer.writerow(['Coordinate system:', 'Beamdyn coordinates']) - elif coordsys == 'VABS/SONATA': - beam_prop_writer.writerow(['Coordinate system:', 'VABS/SONATA coordinates']) - elif coordsys == 'ANBAX': - beam_prop_writer.writerow(['Coordinate system:', 'VABS/SONATA coordinates']) + elif coordsys == 'SONATA': + beam_prop_writer.writerow(['Coordinate system:', 'SONATA coordinates']) + elif coordsys == 'b3_secfem': + beam_prop_writer.writerow(['Coordinate system:', 'SONATA coordinates']) else: beam_prop_writer.writerow(['Coordinate system:', 'to be verified']) diff --git a/SONATA/utl/blade_utl.py b/SONATA/utl/blade_utl.py index f0bfe9b7..bb954729 100644 --- a/SONATA/utl/blade_utl.py +++ b/SONATA/utl/blade_utl.py @@ -46,8 +46,6 @@ def interp_loads(loads, grid_loc): ---------- loads : dict dictionary of the following keys and values, (default=None) - for detailed information see the VABSConfig documentation or the - VABS user manual F : nparray([[grid, F1, F2, F3]]) M : nparray([[grid, M1, M2, M3]]) f : nparray([[grid, f1, f2, f2]]) @@ -63,8 +61,6 @@ def interp_loads(loads, grid_loc): ---------- sectional_load : dict dictionary of the following keys and values, (default=None) - for detailed information see the VABSConfig documentation or the - VABS user manual F : nparray([F1, F2, F3]) M : nparray([M1, M2, M3]) f : nparray([f1, f2, f2]) diff --git a/SONATA/utl_openfast/utl_sonata2beamdyn.py b/SONATA/utl_openfast/utl_sonata2beamdyn.py index e1099564..c675808f 100644 --- a/SONATA/utl_openfast/utl_sonata2beamdyn.py +++ b/SONATA/utl_openfast/utl_sonata2beamdyn.py @@ -26,14 +26,14 @@ if __name__ == '__main__': os.chdir('../..') -# --- Convert from SONATA/VABS coordinates to BeamDyn coordinates ---# +# --- Convert from SONATA coordinates to BeamDyn coordinates ---# def convert_structdef_SONATA_to_beamdyn(cs_pos, SONATA_beam_prop): """ Convert structural characteristics from SONATA definition to BeamDyn definition Inputs: cs_pos - array of radial stations along the span - SONATA_beam_prop - data struct containing the direct results from VABS in SONATA/VABS definition; equiv to job.beam_properties + SONATA_beam_prop - data struct containing the direct results from b3secfem in SONATA definition; equiv to job.beam_properties Outputs: BeamDyn_beam_prop - converted data struct in BeamDyn definition @@ -57,7 +57,7 @@ def convert_structdef_SONATA_to_beamdyn(cs_pos, SONATA_beam_prop): BeamDyn_beam_prop['beam_shear_center'] = np.zeros([len(cs_pos), 2]) # --------------------------------------- # - # retrieve & allocate VABS results + # retrieve & allocate results for i in range(len(SONATA_beam_prop)): if SONATA_beam_prop[i, 1] is not None: BeamDyn_beam_prop['beam_section_mass'][i] = SONATA_beam_prop[i, 1].m00 # mass per unit span (absolute - no transform needed) @@ -71,12 +71,12 @@ def convert_structdef_SONATA_to_beamdyn(cs_pos, SONATA_beam_prop): beam_stiff_init[i, j, :] = np.array(SONATA_beam_prop[i, 1].TS[j, :]) # receive 6x6 timoshenko stiffness matrix beam_inertia_init[i, j, :] = np.array(SONATA_beam_prop[i, 1].MM[j, :]) # receive 6x6 mass matrix # beam_inertia_init[i, j, :] = np.array(SONATA_beam_prop[i, 1].MMatMC[j, :]) # receive 6x6 mass matrix at mass center - else: # If no solution from VABS available - print('Radial station ' + str(cs_pos[i]) + ' did not run successfully in VABS. Check yaml input or change location!') + else: # If no solution from available + print('Radial station ' + str(cs_pos[i]) + ' did not run successfully in b3_secfem. Check yaml input or change location!') # ToDo: instead of break -> skip or interpolate in between working radial stations # --------------------------------------- # - # rotate VABS results from SONATA/VABS def to BeamDyn def coordinate system + # rotate results from SONATA def to BeamDyn def coordinate system # B = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]) # transformation matrix B = np.array([[0, 0, 1], [0, -1, 0], [1, 0, 0]]) # NEW transformation matrix @@ -90,7 +90,7 @@ def convert_structdef_SONATA_to_beamdyn(cs_pos, SONATA_beam_prop): BeamDyn_beam_prop['beam_geometric_center'][i, :] = [-beam_geometric_center_init[i, 0], beam_geometric_center_init[i, 1]] BeamDyn_beam_prop['beam_shear_center'][i, :] = [-beam_shear_center_init[i, 0], beam_shear_center_init[i, 1]] - print('STATUS:\t Structural characteristics of VABS converted from SONATA/VABS to BeamDyn coordinate system definition!') + print('STATUS:\t Structural characteristics converted from SONATA to BeamDyn coordinate system definition!') return BeamDyn_beam_prop @@ -223,7 +223,7 @@ def write_beamdyn_prop(folder, flags_dict, wt_name, radial_stations, beam_inertia[i, j, 0], beam_inertia[i, j, 1], beam_inertia[i, j, 2], beam_inertia[i, j, 3], beam_inertia[i, j, 4], beam_inertia[i, j, 5])) file.write('\n') - # ToDO: check correct translation of stiffness and mass matrices from VABS and anbax !!! + # ToDO: check correct translation of stiffness and mass matrices from SONATA and b3_secfem !!! file.close() print('Finished writing BeamDyn_Blade File') @@ -274,7 +274,7 @@ def write_beamdyn_viscoelastic(folder, flags_dict, wt_name, radial_stations, curr_stiff[j, 3], curr_stiff[j, 4], curr_stiff[j, 5])) file.write('\n') - # ToDO: check correct translation of stiffness and mass matrices from VABS and anbax !!! + # ToDO: check correct translation of stiffness and mass matrices from SONATA and b3_secfem !!! file.close() print('STATUS:\t Finished writing BeamDyn_Blade_Viscoelastic file.') diff --git a/SONATA/vabs/__init__.py b/SONATA/vabs/__init__.py deleted file mode 100755 index 01a6b2b3..00000000 --- a/SONATA/vabs/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Mar 24 10:51:07 2017 - -@author: TPflumm -""" diff --git a/docs/.keep b/docs/.keep deleted file mode 100644 index e69de29b..00000000 diff --git a/docs/Makefile b/docs/Makefile deleted file mode 100644 index ab77832e..00000000 --- a/docs/Makefile +++ /dev/null @@ -1,20 +0,0 @@ -# Minimal makefile for Sphinx documentation -# - -# You can set these variables from the command line. -SPHINXOPTS = -SPHINXBUILD = sphinx-build -SPHINXPROJ = SONATA -SOURCEDIR = source -BUILDDIR = build - -# Put it first so that "make" without argument is like "make help". -help: - @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) - -.PHONY: help Makefile - -# Catch-all target: route all unknown targets to Sphinx using the new -# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -%: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/Pflumm, T. - A Preprocessor for Parametric Composite Rotor Blade Cross-Sections (2018, ERF).pdf b/docs/Pflumm, T. - A Preprocessor for Parametric Composite Rotor Blade Cross-Sections (2018, ERF).pdf deleted file mode 100644 index 97b15445..00000000 Binary files a/docs/Pflumm, T. - A Preprocessor for Parametric Composite Rotor Blade Cross-Sections (2018, ERF).pdf and /dev/null differ diff --git a/docs/developer-guide.md b/docs/developer-guide.md deleted file mode 100644 index 648b73e3..00000000 --- a/docs/developer-guide.md +++ /dev/null @@ -1,80 +0,0 @@ -# Developer Guide - -This document provides information for developers working on SONATA. - - - -## Git Best Practices - -### Commit early and often - -Git only takes full responsibility for your data when you commit. If you fail to commit and then do something poorly thought out, you can run into trouble. Additionally, having periodic checkpoints means that you can understand how you broke something. - -### Do make useful commit messages - -Creating insightful and descriptive commit messages is one of the best things you can do for others who use the repository. It lets people quickly understand changes without having to read code. When doing history archeology to answer some question, good commit messages likewise become very important. - -The normal git rule of using the first line to provide a short (50-72 character) summary of the change is also very good. - -### Use issue tracking - -Issues are a great way to keep track of tasks, enhancements, and bugs for our project. They’re kind of like email—except they can be shared and discussed with the rest of your team. - -If you create an Issue make sure, one can reproduce your bug/feature/findings by adding at least your input files /configuration and a screenshot. - - - - - -## Code Documentation - -Our SONATA project uses the Syntax and best practices for docstrings to be used with the numpydoc extension for [Sphinx](http://sphinx-doc.org/). Allowing our tool to produce well-formatted reference guides. Our docstring standard uses [re-structured text (reST)](http://docutils.sourceforge.net/rst.html) syntax and is rendered using [Sphinx](http://sphinx.pocoo.org/) (a pre-processor that understands the particular documentation style we are using). Keep in mind that the length of docstring lines should be kept to 75 characters to facilitate reading the docstrings in text terminals. - -The complete Styleguide and can be found here: https://numpydoc.readthedocs.io/en/latest/format.html with an of snippets with [example.py](example.py) - -The following example of the CBM method cbm_run_vabs illustrates the idea: - -```python - def cbm_run_vabs(self, jobid=None, rm_vabfiles=True, ramdisk=False): - '''CBM method to run the solver VABS (Variational Asymptotic Beam - Sectional Analysis). Note that this method is designed to work if - VABSIII is set in the PATH variable. For Users at the TUM-HT please load - the vabs module beforehand. - - Parameters - ---------- - jobid : string, optional - assign a unique ID for the job. If no jobid is assigned the - isoformat of datetime with microseconds is used - rm_vabfiles : bool, optional - removes VABS files after the calculation is completed and - the results are stored. - ramdisk : bool, optional, - Instead of storing the writing and reading the vabs job directory, - the ramdisk "/tmpfs/username" is used. This options is currently - designed for linux users make sure to mount it beforehand with to - assign 200MB of Memory to the virtual drive. - >>> sudo mount -t tmpfs -o size=200M none /tmpfs/username - - Returns - ---------- - None : everything is stored within the CBM instance - - Examples - ---------- - >>> job.cbm_run_vabs(rm_vabfiles=True, ramdisk=True) - - ''' -``` - -## Useful References: -* [Numpydoc docstring guide](https://numpydoc.readthedocs.io/) - -* [PythonOCC](http://www.pythonocc.org/) -* [PythonOCC API Documentation](http://api.pythonocc.org/) -* [OpenCascadeTechnology Documentation](https://www.opencascade.com/doc/occt-6.9.1/refman/html/index.html) - -* [openMDAO](http://openmdao.org/) -* [OpenMDAO Documentation](http://openmdao.org/twodocs/versions/latest/) - -* [Dymoresolutions User's Manual](http://www.dymoresolutions.com/UsersManual/UsersManual.html) diff --git a/docs/getting_started.md b/docs/getting_started.md deleted file mode 100644 index bad55622..00000000 --- a/docs/getting_started.md +++ /dev/null @@ -1 +0,0 @@ -# Getting Started diff --git a/docs/img/2dmesh.png b/docs/img/2dmesh.png deleted file mode 100644 index 14a755d3..00000000 Binary files a/docs/img/2dmesh.png and /dev/null differ diff --git a/docs/img/2dtopo.png b/docs/img/2dtopo.png deleted file mode 100644 index f7c650ab..00000000 Binary files a/docs/img/2dtopo.png and /dev/null differ diff --git a/docs/img/3dtopo.png b/docs/img/3dtopo.png deleted file mode 100644 index 82ae7118..00000000 Binary files a/docs/img/3dtopo.png and /dev/null differ diff --git a/docs/img/CS0.png b/docs/img/CS0.png deleted file mode 100644 index b78f9e9e..00000000 Binary files a/docs/img/CS0.png and /dev/null differ diff --git a/docs/img/CS2.png b/docs/img/CS2.png deleted file mode 100644 index 593fedb3..00000000 Binary files a/docs/img/CS2.png and /dev/null differ diff --git a/docs/img/CS4.png b/docs/img/CS4.png deleted file mode 100644 index a217539b..00000000 Binary files a/docs/img/CS4.png and /dev/null differ diff --git a/docs/img/LE.png b/docs/img/LE.png deleted file mode 100644 index 87bdd6f3..00000000 Binary files a/docs/img/LE.png and /dev/null differ diff --git a/docs/img/acknowledgment.png b/docs/img/acknowledgment.png deleted file mode 100644 index d263167e..00000000 Binary files a/docs/img/acknowledgment.png and /dev/null differ diff --git a/docs/img/beam_definition.pdf b/docs/img/beam_definition.pdf deleted file mode 100644 index 5418fa6b..00000000 Binary files a/docs/img/beam_definition.pdf and /dev/null differ diff --git a/docs/img/environment.png b/docs/img/environment.png deleted file mode 100644 index b90f1d1b..00000000 Binary files a/docs/img/environment.png and /dev/null differ diff --git a/docs/img/layup.png b/docs/img/layup.png deleted file mode 100644 index da6065d4..00000000 Binary files a/docs/img/layup.png and /dev/null differ diff --git a/docs/img/logo.png b/docs/img/logo.png deleted file mode 100644 index fefacc22..00000000 Binary files a/docs/img/logo.png and /dev/null differ diff --git a/docs/img/logo_wframe.png b/docs/img/logo_wframe.png deleted file mode 100644 index c0f9af3e..00000000 Binary files a/docs/img/logo_wframe.png and /dev/null differ diff --git a/docs/img/mapping.png b/docs/img/mapping.png deleted file mode 100644 index 1059e7a9..00000000 Binary files a/docs/img/mapping.png and /dev/null differ diff --git a/docs/img/plot_sections.png b/docs/img/plot_sections.png deleted file mode 100644 index b9dfed3d..00000000 Binary files a/docs/img/plot_sections.png and /dev/null differ diff --git a/docs/img/post_3dtopo.png b/docs/img/post_3dtopo.png deleted file mode 100644 index 0540da9a..00000000 Binary files a/docs/img/post_3dtopo.png and /dev/null differ diff --git a/docs/intro.md b/docs/intro.md deleted file mode 100644 index dad0cd1b..00000000 --- a/docs/intro.md +++ /dev/null @@ -1,213 +0,0 @@ - -## Introduction: - -The large number of constraints and design drivers from various disciplines makes the helicopter rotor blade development process difficult, time consuming and costly. -The entire design process represents a classical aeroelastic problem, where the aerodynamic behavior, the structural elasticity and vibrational dynamics have to be studied simultaneously. The behavior can therefore not be examined with separate analysis of the different disciplines [1]. The integration of all the appropriate disciplines in the design process implies not only limitations on the design from various disciplines, but also defining and accounting for interactions so that the disciplines influence design decisions simultaneously rather than sequentially [2]. -Historically, the design and development of improved or entirely new rotor blades is conducted by departments in a company that maintain their separate simulation codes for performing their specific tasks. The aerodynamics department is responsible for performance calculations, aero-acoustics, rotor-wake interaction, unsteady airload prediction and computational fluid dynamics while the dynamics department focuses on rotor vibratory loads, stability and aeroelastic models [1]. The structural department determines the elastic properties as well as strength and fatigue characteristics. A Blade and Rotor Design Department often bundles the different aspects while considering materials, manufacturability, maintainability and safety requirements. [1] -This modular approach narrows the scope of solutions, because each department focuses on individual objectives satisfied by individual design parameters. Mutual interactions can only be covered by numerous iterations. -In contrast to that, a multidisciplinary approach offers a more systematic development process that is able to design a better helicopter rotor [2]. Because of the impact the rotor behavior has on the overall performance of the helicopter and on customer noticeable vibratory characteristics, rotor aeroelastic effects should be considered in the earliest stages of the design process [3]. -In the last 25 years, researchers have repeatedly stated the need for a design methodology and optimization framework that combines computational efficiency of a beam description in aeromechanic analysis with a rotor blade structural model that is capable at describing realistic composite rotor blade cross-sections with respect to the structural properties, applied load, stress and strain distributions as well as design constraints [4-6]. - -## Framework: - -Our multidisciplinary rotor blade design framework is named SONATA (Structural Optimization and Aeroelastic Analysis) and is illustrated in the following figure 1. Like most environments it comprises of **three** main components that are wrapped into an optimization framework. - - - -Fig. 1: SONATA: Multidisciplinary Rotor Blade Design Environment for Structural Optimization and Aeroelastic Analysis embedded in OpenMDAO. - -1. As a **first** component, the current state of the art involves an aeromechanical analysis of rotorcraft blades which includes flexible multibody dynamics, nonlinear finite elements and various rotorcraft aerodynamic models. They are often referred to as Comprehensive Analysis. Examples are the widely used Comprehensive Analytical Model of Rotorcraft Aerodynamics and Dynamics II (CAMRAD II) [9] and the software Dymore [10] beyond several others. Both of these codes are presently in use in the rotorcraft industry, academic institutions and government laboratories. The quality of the predictions have been documented in numerous publications. In our SONATA environment Dymore was chosen as aeromechanic tool for both a dynamic analysis in the time domain as well a modal analysis within the frequency domain. In this context classical 1D-beam elements are used to describe the rotor blade due to the much simpler mathematical formulation and reduced computational effort compared to a full three-dimensional finite element model of the composite rotor blade [11]. Typically, this approach decouples the realistic composite blade definition and the manufacturability constraints from the aeromechanic analysis and the predesign of structural blade properties. That way, problems in the blade design cannot be discovered until later in the process where changes are costly and time consuming [12]. - -2. Although the three-dimensional finite element method is the most accurate approach to model realistic rotor blades, it is still not appropriate for the use in rotor blade predesign [11, 13]. The slender characteristic of rotor blades allows the simplification to treat them as one-dimensional body [14] Cesnik and Hodges [15] formulated the Variational Asymptotic Beam Sectional Analysis (VABS) to accurately represent the behavior that is associated with the reduction of two-dimensions. In other words, this method splits the three-dimensional elastic problem into a two-dimensional linear cross-section analysis and a one-dimensional nonlinear beam analysis, which is able to consider initially twisted and curved, anisotropic, non-homogeneous materials to model general composite cross-sectional geometries [13, 15]. VABS is the **second** component of our environment. In the last 20 years, VABS and its variations have become a popular tool in rotor blade predesign and multidisciplinary rotor design optimization and their accuracy and efficiency has been validated in numerous publications [3, 15, 16]. - -3. Consequently, most researches have developed individual parametric mesh generators for the cross-sectional analysis, that reduces their structural model to few design variables in the process. Such a preprocessor for parametric composite rotor blade cross-sections is referred to as *SONATA-CBM* in this framework. It is the **third** component of the SONATA environment. - -Last but not least, the tree components are managed by an environment where design variables and objectives can be defined, constraints to be applied and solvers to be launched. The **SONATA** framework uses [OpenMDAO](http://openmdao.org/) [7, 8, 27], an open-source computing platform for system analysis and multidisciplinary optimization, written in Python. It allows the user to break down the structure of complex optimization tasks into a hierarchic manner while managing the numerical methods. A Python-based wrapper for Dymore has been developed to integrate the dynamic and modal analysis into the OpenMDAO-driven optimizations. Consequently *SONATA-CBM* has been written in Python using the Python wrapper for the CAD-Kernel Opencascade [pythonOCC](http://www.pythonocc.org/ ). - -**Why Python?** - -- Python can be easy to pick up whether you're a first time programmer or you're experienced with other languages. - -- Python is developed under an OSI-approved open source license, making it freely usable and distributable, even for commercial use. - -- The Python Package Index (PyPI) hosts thousands of third-party modules for Python. Both Python's standard library and the community-contributed modules allow for endless possibilities. Two of the most important python modules used in SONATA are the openMDAO and the pythonocc module. - -- [openMDAO](http://openmdao.org/) is an open-source high-performance computing platform for systems analysis and multidisciplinary optimization, written in Python. - -- [pythonOCC](http://www.pythonocc.org/ ) is a python library whose purpose is to provide 3D modeling features. It is intended to developers who aim at developing CAD/PDM/PLM applications. - - -### 1. DYMORE (PYMORE): - -Dymore is a finite element based multibody dynamics code for the comprehensive modeling of flexible multibody systems. Dymore has been implemented in industrial settings, it is presently used by a number of rotorcraft, wind energy, and aerospace companies and research laboratories. Dymore has been developed by Prof. [Bauchau](http://www.dymoresolutions.com/resume/Resume.html) and his students. The documentation can be found in the [user's manual](http://www.dymoresolutions.com/UsersManual/UsersManual.html). - -A Python-based wrapper for Dymore has been developed (called Pymore) to integrate the dynamic and modal analysis into the OpenMDAO-driven optimizations. This module is developed by W. Garre and is hosted @ [https://gitlab.lrz.de/wgarre/Pymore](https://gitlab.lrz.de/wgarre/Pymore) - - - -### 2. VABS (Variational Asymptotic Beam Sectional Analysis): - -"VABS (Variational Asymptotic Beam Sectional Analysis) is a code implementing the various beam theories based on the concept of simplifying the original nonlinear three-dimensional (3D) analysis of slender structures into a two-dimensional (2D) cross-sectional analysis and a one- dimensional (1D) nonlinear beam analysis using a powerful mathematical method, the variational asymptotic method. VABS takes a finite element mesh of the cross section including all the details of geometry and material as inputs to calculate the sectional properties including structural properties and inertial properties. These properties are needed for the 1D beam analysis to predict the global behavior of the slender structure. The 3D pointwise displacement/strain/stress distribution within the structure can also be recovered based on the global behavior of the 1D beam analysis." from [[VABS Manual for Users](docs/man/VABS-Manual.pdf)] - -Thanks to [Analyswift](http://analyswift.com/vabs-cross-sectional-analysis-tool-for-composite-beams/), the distributer of VABS, the Technical University of Munich has received academic licenses free of charge. - - -### 3. SONATA-CBM: -*SONATA-CBM'*s composite topology generation originates from an arbitrary closed curve that can be obtained from various input formats that range from airfoil coordinate tables over a 3D CAD rotor blade surface definition (.step or .iges) with radial station to a parameterized rotor blade with twist, planform, airfoil and chord-line distribution. In the case of the latter two, the 3D surface is intersected at a certain radial station to obtain once again a two-dimensional outer boundary of the cross-section. Figure 2 shows the resulting parameterized 3D surface of the UH-60A rotor blade with a cross-section topology at radial station R = 2000 mm.\\ -While the following methodology is shown with the example of the UH-60A rotor-blade, it should be noted that this procedure can be applied to any closed curve cross-section, and therefore be also used to model rotor blade root sections or any other composite beam cross-sections. - - - -Fig. 2: Parameterized 3D surface of the UH-60A rotor blade created with twist, planform, airfoil and axis information from Davis [29] - -#### Topology Generation - -The process behind the composite topology generation is derived from the manufacturing process, where the layers are placed on top of each other in negative molds in a consecutive manner to avoid complex constraints in the optimization and to keep the solution within proper bounds. Each layer has an assigned material with start and end coordinates, a thickness and fiber orientation (see table 1). Every parameter or groups of them can serve as design variable in the later optimization. After the layup process on top of the outer boundary curve is completed, webs are introduced and subsequently new closed curved geometries are generated where the layup procedure is repeated. Cavities can be filled with core materials and additional trim masses can be inserted. -At first the outer boundary curve, represented as counterclockwise sets of consecutive B-splines, is defined in curve coordinates **s** between zero and one. The origin is typically located at the trailing edge (TE). The curve coordinate system propagates through the layers with an interval tree structure. It allows to efficiently find the intervals/layers that overlap and locate the corresponding coordinate for each layer. -Subsequently, each layer is generated by the following consecutive steps. - -- Determine the relevant set of underlying B-Splines between *start* and *end* coordinate of the layer using an interval tree data structure. -- Discretize the set of B-Splines and perform an parallel offset to return an approximate representation of all points with a given thickness of each layer. -- Generate a new set of B-Splines by interpolation and add smooth layer cutoffs to connect the lower and upper set of B-Splines if necessary. - -In table 1 the layup definition of the cross-section, illustrated in figure 3, is displayed. Note that the shown genetic composite cross-section of the UH-60A serves as demonstration of the modeling capabilities. - - - - - -Table 1: Layup definition of figure 3 - - - -Fig. 3: Topology definition of a generic composite UH-60A rotor blade cross section. - -​ - -Fig. 4: SONATA-CBM discretization of a generic composite UH-60A rotor blade cross-section in reference to [3] to illustrate the modeling capabilities. - -The first set of layers are grouped into Segment 0. The first layer that is generated is a steel erosion protection strip that ranges from coordinate 0.44 to 0.56 with a thickness of 0.82mm. Because of the isotropic material used, the orientation can be neglected for this layer. The material ID represents a reference index of an associated material database. The next 4 layers define the skin of the rotor blade placed in both 0°and +-45° orientation on top of each other. The layers Spar 1 to Spar 7 are unidirectional carbon fiber composite layers that generate a C type spar with ply drops in the leading edge region of the cross-section. -Once the first set of layers (Segment 0) has been created, webs are introduced to the structure. They are defined as straight line between two positions. In this example the first web ranges from coordinate 0.43 to 0.57 while the second is placed behind from 0.30 to 0.70. The three newly generated closed curved geometries are used to repeat the layup procedure. During the manufacturing process this translates to a process of wrapping plies around a core. A core material is assigned to Segment 1 and 3 that fills up the remaining cavity. Segment 2 consists of four carbon fiber layers of different orientation from 0 to 1 to generate a hollow box spar. After the layup is defined a trim mass can be placed on top of the existing layers and will be integrated in the structure during the discretization. - -#### Discretization: - -The discretization follows the topology generation procedure, yet in a reversed direction with respect to the layup definition, starting from the innermost layers and moving outwards. Each layer is meshed by an orthogonal projection with corner style differentiation. -Figure 4 shows the final result of the described procedure. - -Each layer can be described by two sets of B-splines, the inner a_bsplines and outer b_bsplines. The nodes placed on them are called accordingly a_nodes and b_nodes. -The following procedure is applied to each layer, starting at the innermost, and moving outwards. - -- Determine existing a_nodes based on the intervaltree structure of the layup. If sections on the a_bsplines are found with no preexisting nodes, distribute new nodes equidistantly. -- Create an orthogonal projection of each a_node onto the set of b_bsplines. If two or more projections are found determine the angle alpha and the number of potential b_bsplines corners between them. -- Based on a critical angle alpha_crit and the number of exterior corners determine the corner style and as a consequence the meshing procedure. In figure 5-7 the first 6 different corner styles are shown. -- After all nodes are placed on both sets of B-splines, they are connected to form cells with associated material and ply angles. -- In subsequent steps sharp cells, large aspect-ratio cells and cell angles are modified to improve mesh quality. - -As soon as every layer of the segment is meshed, the remaining cavities are triangulated using Shewchuk [30] algorithm with an area constraint. To avoid hanging nodes between two neighboring segments, the cells are consolidated on web interfaces. - - - - - -Fig. 5-7: Corner-style **0**: no exterior corner on bbsplines and α > αcrit ; Corner-style **1**: one exterior corner on bbsplines and α > αcrit. Corner-style **2**: no exterior corner on bbsplines and α < αcrit ; Corner-style **3**: one exterior corner on bbsplines and α < αcrit. Corner-style **4**: two exterior corners on bbsplines and α < αcrit ; Corner-style **5**: three exterior corners on bbsplines and α < αcrit. - -In a final step, the previously defined trim mass is integrated into the described mesh by mapping existing nodes onto the trim mass contour. The corresponding algorithm is schematically illustrated in Figure [8] and described below: - -- Determine the number of inner nodes of the intersected cells. -- Move the inner nodes of the cells marked *1* along the cell edge with shortest distance to the intersecting curve. -- Move the remaining inner nodes of the cells marked *2* along the cell edge with shortest distance to the intersecting curve. -- Move the outer node of the cells marked *3* along the edge direction onto the intersecting curve. -- Delete cells marked *3* and *4.* -- Use the boundary nodes as starting point for the inner triangulation. - - - -Fig. 8-9: (Left) Mapping algorithm to integrate cuves into an existing mesh. (Right) Leading edge region of figure 4 showing the plydrops of the C-Spar and the integration of the trim mass into the existing mesh. - -The final result is displayed in the magnified cutout of the leading edge region in figure 9. Finally, the VABS input files are generated from the mesh together with the material information from an associated database. - -To verify the resulting stiffness properties, simple benchmark testcases for isotropic and anisotropic box-beam cross-sections have been set up and compared to results from [31]. Moreover, the rotor blades of the institute's high altitude synchropter UAV (AREA) [32, 33] have been reengineered with SONATA-CBM and compared to experimental results from Suesse [34]. - - -## Referencencs: -[1] Tarzanin, F. and Young, D., “Boeing rotorcraft experience with rotor design and optimization,” 7th -AIAA/USAF/NASA/ISSMO Symp. Multidiscip. Anal. Optim., American Institute of Aeronautics and Astronautics, Reston, Virigina, sep 1998. [http://arc.aiaa.org/doi/abs/10.2514/6.1998-4733](http://arc.aiaa.org/doi/abs/10.2514/6.1998-4733). - -[2] Adelman, H. M. and Mantay, W. R., “Integrated Multidisciplinary Optimization of Rotorcraft: A Plan for Development,” Tech. rep., NASA, 1989. - -[3] Rohl, P. J., Kumar, D., Dorman, P., Sutton, M., and Cesnik, C. E. S., “A Composite Rotor Blade Structural -Design Environment for Aeromechanical Assessments in Conceptual and Preliminary Design,” American Helicopter Society 68th Annual Forum, American Helicopter Society, 2012. [http://ebooks.cambridge.org/ref/id/CBO9781107415324A009](http://ebooks.cambridge.org/ref/id/CBO9781107415324A009). - -[4] Friedmann, P. P., “Helicopter Vibration Reduction Using Structural Optimization with Aeroelas- -tic/multidisciplinary Constraints - A Survey,” Journal of Aircraft, Vol. 28, No. 1, jan 1991, pp. 8–21. http://dx.doi.org/10.2514/3.45987](http://dx.doi.org/10.2514/3.45987). - -[5] Weller, W. H. and Davis, M. W., “Wind Tunnel Tests of Helicopter Blade Designs Optimized for Minimum Vibration,” American Helicopter Society 44th Annual Forum, 1988. [ - -[9]Johnson, W., “A History of Rotorcraft Comprehensive Analyses,” American Helicopter Society 60th Annual Forum, 2013. - -[10] Bauchau, O., Bottasso, C., and Nikishkov, Y., “Modeling rotorcraft dynamics with finite element multi- -body procedures,” Mathematical and Computer Modelling, Vol. 33, No. 10-11, 2001, pp. 1113–1137. - -[11] Datta, A. and Johnson, W., “Three-Dimensional Finite Element Formulation and Scalable Domain Decomposition for High-Fidelity Rotor Dynamic Analysis,” Journal of the American Helicopter Society, 2011. - -[12] Rohl, P., Dorman, P., Sutton, M., Kumar, D., and Cesnik, C., “A Multidisciplinary Design Environment for Composite Rotor Blades,” 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, No. April, American Institute of Aeronautics and Astronautics (AIAA), Reston, Virigina, apr 2012, pp. 1–15. http://dx.doi.org/10.2514/6.2012-1842](http://dx.doi.org/10.2514/6.2012-1842). - -[13] Li, L., Structural Design of Composite Rotor Blades with Consideration of Manufacturability, Durability, and Manufacturing Uncertainties, Ph.d. thesis, Georgia Institute of Technology, 2008. - -[14] Yeo, H., Truong, K.-V., and Ormiston, R. A., “Asessment of 1D Versus 3D Methods for Modeling Rotor Blade Structural Dynamics,” AIAA, 2010. - -[15] Cesnik, C. E. S. and Hodges, D. H., “VABS: A New Concept for Composite Rotor Blade Cross-Sectional Modeling,” American Helicopter Society 51st Annual Forum, 1995. - -[16] Cesnik, C., Mok, J., Parikh, A., and Shin, S., “Optimum Design Framework for Integrally Twisted Helicopter -Blades,” 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Ameri- -can Institute of Aeronautics and Astronautics (AIAA), apr 2004. - -[17] Kumar, D. and Cesnik, C. E., “Optimization Framework for the Dynamic Analysis and Design of Active Twist Rotors,” American Helicopter Society 68th Annual Forum, 2012. - -[18] Kumar, D. and Cesnik, C. E., “New Hybrid Optimization for Design of Active Twist Rotors,” 54th -AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, American Institute of Aeronautics and Astronautics (AIAA), 2013. - -[19] Silva, C. and Johnson, W., “Multidisciplinary Conceptual Design for Reduced-Emission Rotorcraft,” AHS -Specialists Conference on Aeromechanics Design for Transformative Vertical Flight, AHS, San Francisco, Cali- -fornia, jan 2018. - -[20] Meyn, L., “Rotorcraft Optimization Tools: Incorporating Rotorcraft Design Codes into Multi-Disciplinary Design, Analysis, and Optimization,” . - -[21] Glaz, B., Friedmann, P. P., and Liu, L., “Helicopter Vibration Reduction throughout the Entire Flight Enve- -lope Using Surrogate-Based Optimization,” Journal of the American Helicopter Society, Vol. 54, No. 1, 2009. -[22] Glaz, B., Friedmann, P. P., Liu, L., Kumar, D., and Cesnik, C. E. S., “The AVINOR Aeroelastic Simulation -Code and its Application to Reduced Vibration Composite Rotor Blade Design,” 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Confer- -ence, 2009. - -[23] Yu, W., Volovoi, V., Hodges, D. H., and Hong, X., “Validation of the variational asymptotic beam sectional analysis,” AIAA Journal, Vol. 40, No. 10, jan 2002, -pp. 2105–2112. - -[24] Tian, S., Liu, X., and Yu, W., “PreVABS, -https://cdmhub.org/resources/1597,” Nov 2017. - -[25] Fanjoy, D. and Crossley, W., “Using a Genetic Algorithm for Structural Topology Design of Helicopter Ro- -tor Blades,” 19th AIAA Applied Aerodynamics Conference, American Institute of Aeronautics and Astronautics (AIAA), jun 2001. - -[26] Blasques, J. P., “Multi-material topology optimization of laminated composite beams with eigenfrequency constraints,” Composite Structures, Vol. 111, 2013, pp. 45 – 55. - -[27] Heath, C. M. and Gray, J. S., “OpenMDAO: Framework for Flexible Multidisciplinary Design, Analysis and -Optimization Methods,” 8th AIAA Multidisciplinary Design Optimization Specialist Conference (MDO), Hon- -olulu, Hawaii, 2012, pp. 1–13. - -[28] Paviot, T., “pythonOCC, 3D CAD/CAE/PLM development framework for the Python programming language,http://www.pythonocc.org/,” . - -[29] Davis, S. J., “Predesign Study For a Modern 4-Bladed Rotor for the RSRA,” Tech. Rep. 16155, NASA, 1981. - -[30] Shewchuk, J. R., “Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator,” Applied Computational Geometry: Towards Geometric Engineering, Vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag, 1996, pp. 203–222. - -[31] Popescu, B. and Hodges, D. H., “On asymptotically correct Timoshenko-like anisotropic beam theory,” Inter national Journal of Solids and Structures, Vol. 37, No. 3, 2000, pp. 535 – 558. - -[32] Barth, A., Spiess, C., Kondak, K., and Hajek, M., “Design, Analysis and Flight Testing of a High Altitude -Synchropter UAV,” American Helicopter Society 74th Annual Forum, 2018. - -[33] Pflumm, T., Barth, A., Kondak, K., and Hajek, M., “Auslegung und Konstruktion eines Hauptrotorblattes fuer ein in extremen Flughoehen operierendes Drehfluegel UAV,” Deutscher Luft- und Raumfahrtkongress 2015, Rostock, Germany, 2015. - -[34] Suesse, S. and Hajek, M., “Rotor Blade Displacement and Load Estimation with Fiber-Optical Sensors for a Fu ture Health and Usage Monitoring System,” American Helicopter Society 74th Annual Forum, 2018. - -[35] McColl, C., Palmer, D., Chierichetti, M., Bauchau, O. A., and Ruzzene, M., “Comprehensive UH-60 Loads -Model Validation,” AHS Forum, 2010. [36]Bowen-Davies, G. M., Performance and Loads of Variable Top Speed Rotorcraft at High Advance Ratios, Ph.D. thesis, University of Maryland, 2015. diff --git a/docs/make.bat b/docs/make.bat deleted file mode 100644 index 9534b018..00000000 --- a/docs/make.bat +++ /dev/null @@ -1,35 +0,0 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=sphinx-build -) -set SOURCEDIR=source -set BUILDDIR=build - -if "%1" == "" goto help - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The 'sphinx-build' command was not found. Make sure you have Sphinx - echo.installed, then set the SPHINXBUILD environment variable to point - echo.to the full path of the 'sphinx-build' executable. Alternatively you - echo.may add the Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.http://sphinx-doc.org/ - exit /b 1 -) - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% - -:end -popd diff --git a/docs/man/PatranInputFile.pdf b/docs/man/PatranInputFile.pdf deleted file mode 100644 index f3b12696..00000000 Binary files a/docs/man/PatranInputFile.pdf and /dev/null differ diff --git a/docs/man/VABS-Manual.pdf b/docs/man/VABS-Manual.pdf deleted file mode 100644 index 37557c9e..00000000 Binary files a/docs/man/VABS-Manual.pdf and /dev/null differ diff --git a/docs/source/basic_userguide.rst b/docs/source/basic_userguide.rst deleted file mode 100644 index f8600cbe..00000000 --- a/docs/source/basic_userguide.rst +++ /dev/null @@ -1,26 +0,0 @@ -Basic User Guide -================= - - -SONATA.Blade ------------- - - -SONATA.CBM ----------- - - -VABS ------ - - -ANBAX ------ - - -DYMORE ------- - - -OPENMDAO --------- diff --git a/docs/source/cite.rst b/docs/source/cite.rst deleted file mode 100644 index db48df6d..00000000 --- a/docs/source/cite.rst +++ /dev/null @@ -1,2 +0,0 @@ -How to cite SONATA ------------------- diff --git a/docs/source/conf.py b/docs/source/conf.py deleted file mode 100644 index e1cec728..00000000 --- a/docs/source/conf.py +++ /dev/null @@ -1,76 +0,0 @@ -# Configuration file for the Sphinx documentation builder. -# -# This file only contains a selection of the most common options. For a full -# list see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html - -# -- Path setup -------------------------------------------------------------- - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# -import os -import sys -sys.path.insert(0, os.path.abspath('.')) - - -# -- Project information ----------------------------------------------------- - -project = 'sonata' -copyright = '2020, TPflumm' -author = 'TPflumm' - -# The full version, including alpha/beta/rc tags -release = 'v0.4.0' - - -# -- General configuration --------------------------------------------------- - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. -extensions = [ - "sphinx.ext.autodoc", - "sphinx.ext.autosummary", - "sphinx.ext.coverage", - "sphinx.ext.mathjax", - "sphinx.ext.viewcode", - "sphinx.ext.napoleon", -] -napoleon_google_docstring = False -napoleon_use_param = False -napoleon_use_ivar = True - - - - - -# Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -# This pattern also affects html_static_path and html_extra_path. -exclude_patterns = [] - - -# -- Options for HTML output ------------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# -html_theme = "sphinx_rtd_theme" -html_theme_options = { - 'logo_only': False -} - -html_logo = "../img/logo_wframe.png" -html_show_sourcelink = True - - - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] diff --git a/docs/source/getting_started.rst b/docs/source/getting_started.rst deleted file mode 100644 index 8ff09491..00000000 --- a/docs/source/getting_started.rst +++ /dev/null @@ -1,27 +0,0 @@ -Getting Started -============================ - -Prerequesites -------------- - - - -Installation -------------- - - -Install ``clana`` and execute the example: - -:: - - $ pip install s. - $ python mnist_example.py - -This will generate the clana files. - -Usage ------ - - -Input File Format ------------------ diff --git a/docs/source/index.rst b/docs/source/index.rst deleted file mode 100644 index efd38303..00000000 --- a/docs/source/index.rst +++ /dev/null @@ -1,12 +0,0 @@ -Welcome to SONATA’s documentation! -======================================================= -.. toctree:: - :maxdepth: 2 - :caption: Contents: - :numbered: - - getting_started - theory_manual - basic_userguide - readme - cite diff --git a/docs/source/readme.rst b/docs/source/readme.rst deleted file mode 100644 index 44b9e691..00000000 --- a/docs/source/readme.rst +++ /dev/null @@ -1,4 +0,0 @@ -Readme -=========== - -.. mdinclude:: ../../README.md diff --git a/docs/source/source_docs.rst b/docs/source/source_docs.rst deleted file mode 100644 index da6fa8d0..00000000 --- a/docs/source/source_docs.rst +++ /dev/null @@ -1,9 +0,0 @@ -Source Docs -=========== - -.. automodule:: SONATA.classBlade -Indices and tables ------------------- -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` diff --git a/docs/source/theory_manual.rst b/docs/source/theory_manual.rst deleted file mode 100644 index 5c366a49..00000000 --- a/docs/source/theory_manual.rst +++ /dev/null @@ -1,5 +0,0 @@ -Theory Manual -================= - -Coordinate Systems ------------------- diff --git a/environment.yaml b/environment.yaml index 918ad106..28a5d262 100644 --- a/environment.yaml +++ b/environment.yaml @@ -15,8 +15,12 @@ dependencies: - pre-commit - shapely - windIO + - fenics-basix + - fenics-ufl + - fenics-dolfinx + - mpi4py - pip: - - triangle + - OCC variables: PIP_EXTRA_INDEX_URL: https://pypi.org/simple diff --git a/examples/0_beams/0_SONATA_init_box_beam_HT_antisym_layup_15_6_SI_SmithChopra91.py b/examples/0_beams/0_SONATA_init_box_beam_HT_antisym_layup_15_6_SI_SmithChopra91.py index fd00ce8c..55cb9df3 100644 --- a/examples/0_beams/0_SONATA_init_box_beam_HT_antisym_layup_15_6_SI_SmithChopra91.py +++ b/examples/0_beams/0_SONATA_init_box_beam_HT_antisym_layup_15_6_SI_SmithChopra91.py @@ -27,7 +27,7 @@ Analysis Objectives: ===================== - Generate finite element mesh for the composite cross-section -- Compute sectional properties using ANBAX (Asymptotic Numerical Beam Analysis eXtended) +- Compute sectional properties using b3_secfem (Asymptotic Numerical Beam Analysis eXtended) - Evaluate structural response and material distribution - Validate against reference solutions from literature @@ -107,8 +107,8 @@ print("Generating sections and mesh...") job.blade_gen_section(topo_flag=True, mesh_flag=True) -print("Running ANBAX analysis...") -job.blade_run_anbax() +print("Running b3_secfem analysis...") +job.blade_run_b3_secfem() #----------------------------------------- # Generate plots and summary diff --git a/examples/1_IEA15MW/1_sonata_IEA15.py b/examples/1_IEA15MW/1_sonata_IEA15.py index 083e1711..e6675843 100644 --- a/examples/1_IEA15MW/1_sonata_IEA15.py +++ b/examples/1_IEA15MW/1_sonata_IEA15.py @@ -54,8 +54,8 @@ # ===== User defined radial stations ===== # # Define the radial stations for cross sectional analysis (only used for flag_wt_ontology = True -> otherwise, sections from yaml file are used!) -# radial_stations = [0., 0.01, 0.03, 0.05, 0.075, 0.15, 0.25, 0.3 , 0.4, 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.] -radial_stations = np.linspace(0., 1., 5, endpoint=True) +radial_stations = [0., 0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3 , 0.35, 0.4, 0.45, 0.5 , 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1.] +# radial_stations = np.linspace(0., 1., 5, endpoint=True) # ===== Execute SONATA Blade Component Object ===== # # name - job name of current task # filename - string combining the defined folder directory and the job name diff --git a/examples/5_UH60/5_uh60.py b/examples/5_UH60/5_uh60.py index 9265b006..dda1100e 100644 --- a/examples/5_UH60/5_uh60.py +++ b/examples/5_UH60/5_uh60.py @@ -30,7 +30,7 @@ Analysis Objectives: ===================== - Generate finite element mesh for the composite cross-section -- Compute sectional properties using ANBAX (Asymptotic Numerical Beam Analysis eXtended) +- Compute sectional properties using b3_secfem (Asymptotic Numerical Beam Analysis eXtended) - Evaluate structural response and material distribution - Validate against reference solutions from literature @@ -118,8 +118,8 @@ print("Generating sections and mesh...") job.blade_gen_section(topo_flag=True, mesh_flag=True) -print("Running ANBAX analysis...") -# job.blade_run_anbax() # call through beam_struct_eval instead to save outputs +print("Running b3_secfem analysis...") +# job.blade_run_b3_secfem() # call through beam_struct_eval instead to save outputs flags_dict['flag_csv_export'] = flag_csv_export flags_dict['flag_write_BeamDyn_unit_convert'] = flag_write_BeamDyn_unit_convert diff --git a/pyproject.toml b/pyproject.toml index a9973a32..7be047da 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,7 @@ description = "Multidiciplinary Rotor Blade Design Environment for Structural Op readme = "README.md" requires-python = ">=3.9" license = {text = "MIT"} -keywords = ["Helicopter","Rotor Blade","OpenMdao","Preprocessor","SONATA","VABS","ANBA"] +keywords = ["Helicopter","Rotor Blade","OpenMMDAO","Preprocessor","SONATA", "b3_secfem"] classifiers = [ # Optional # How mature is this project? Common values are # 3 - Alpha @@ -45,13 +45,7 @@ dependencies = [ "openmdao", "windio", "matplotlib", - "fenics-basix", - "fenics-ufl", - "fenics-dolfinx", - "mpi4py", - "OCC" ] -# "dolfinx", [project.urls] # Optional #"Homepage" = "https://gitlab.lrz.de/rfeil/SONATA" diff --git a/tests/regression/0_beams/test_beam_0.py b/tests/regression/0_beams/test_beam_0.py index d7a4057e..7128509f 100644 --- a/tests/regression/0_beams/test_beam_0.py +++ b/tests/regression/0_beams/test_beam_0.py @@ -48,7 +48,7 @@ def test_6x6_beam0(self): job = Blade(name=job_name, filename=filename_str, flags=flags_dict, stations=radial_stations) job.blade_gen_section(topo_flag=True, mesh_flag = True) - job.blade_run_anbax() + job.blade_run_b3_secfem() ########################################################################### ######## Checks on if answer looks consistent with previous runs ########## diff --git a/tests/regression/1_iea15mw/ref_iea15mw_bd_blade.dat b/tests/regression/1_iea15mw/ref_iea15mw_bd_blade.dat index 8072558d..171b9c51 100644 --- a/tests/regression/1_iea15mw/ref_iea15mw_bd_blade.dat +++ b/tests/regression/1_iea15mw/ref_iea15mw_bd_blade.dat @@ -2,7 +2,7 @@ Test Format 1 ---------------------- BLADE PARAMETERS -------------------------------------- 1 station_total - Number of blade input stations (-) - 1 damp_type - Damping type: 0: no damping; 1: damped + 1 damp_type - Damping type (switch) 0: none, 1: stiffness-proportional, 2: modal ---------------------- DAMPING COEFFICIENT------------------------------------ mu1 mu2 mu3 mu4 mu5 mu6 (-) (-) (-) (-) (-) (-) @@ -12,16 +12,16 @@ 0.01 0.01 0.01 0.01 0.01 0.01 zeta - Damping coefficients for mode 1 through n_modes ---------------------- DISTRIBUTED PROPERTIES--------------------------------- 0.700000 - 5.4993413123883523e+07 3.7189702120758761e+06 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.9250458136271890e+06 - 3.7189702120374814e+06 4.2374354894044232e+08 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.6846776106581904e+07 - 0.0000000000000000e+00 0.0000000000000000e+00 1.3971807942476038e+10 7.8978507861388350e+08 -3.9263453066697073e+08 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 7.8978507861332667e+08 2.3839010867776895e+09 -1.4564909497568097e+07 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 -3.9263453066717011e+08 -1.4564909497674715e+07 1.3527237632050970e+09 0.0000000000000000e+00 - 8.9250458135987055e+06 1.6846776106596414e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 6.9494752713175327e+07 + 5.4892682885072835e+07 8.1208098153816937e+06 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.0751844881100357e+06 + 8.1208098153816937e+06 5.0850492968031818e+08 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.0776729856298842e+07 + 0.0000000000000000e+00 0.0000000000000000e+00 1.3971807944309111e+10 7.8978507839228690e+08 -3.9263453051052380e+08 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 7.8978507839228690e+08 2.3839010860614142e+09 -1.4564908963989265e+07 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -3.9263453051052380e+08 -1.4564908963989265e+07 1.3527237629568388e+09 0.0000000000000000e+00 + 9.0751844881100357e+06 2.0776729856298842e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 6.9668125909566522e+07 - 2.2012579405916057e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -5.0210264104882775e+01 - 0.0000000000000000e+00 2.2012579405916057e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 6.8686898749778722e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 2.2012579405916057e+02 5.0210264104882775e+01 -6.8686898749778722e+00 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 5.0210264104882775e+01 1.1480879832955482e+02 -2.4912390503192805e+00 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 -6.8686898749778722e+00 -2.4912390503192805e+00 1.7736253772657896e+01 0.0000000000000000e+00 - -5.0210264104882775e+01 6.8686898749778722e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3254505210221288e+02 + 2.2012579406061363e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -5.0210264103893138e+01 + 0.0000000000000000e+00 2.2012579406061363e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 6.8686898753828114e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 2.2012579406061363e+02 5.0210264103893138e+01 -6.8686898753828114e+00 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 5.0210264103893138e+01 1.1480879833022608e+02 -2.4912390500315311e+00 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -6.8686898753828114e+00 -2.4912390500315311e+00 1.7736253772750903e+01 0.0000000000000000e+00 + -5.0210264103893138e+01 6.8686898753828114e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3254505210297700e+02 diff --git a/tests/regression/1_iea15mw/test_iea15mw.py b/tests/regression/1_iea15mw/test_iea15mw.py index a4c0a5b9..27ad7d23 100644 --- a/tests/regression/1_iea15mw/test_iea15mw.py +++ b/tests/regression/1_iea15mw/test_iea15mw.py @@ -334,9 +334,9 @@ def test_external_mesh_iea15mw(self): job2.blade_custom_mesh(nodes, cells, MatID, split_quads=True, theta_11=theta_11, theta_3=None) - job.blade_run_anbax() + job.blade_run_b3_secfem() - job2.blade_run_anbax() + job2.blade_run_b3_secfem() plt.close('all') matplotlib.use(original_backend) @@ -345,7 +345,7 @@ def test_external_mesh_iea15mw(self): # 6x6 mass matrix npt.assert_allclose(job.beam_properties[i, 1].MM, job2.beam_properties[i, 1].MM, rtol=1e-5) - + # 6x6 timoshenko stiffness matrix npt.assert_allclose(job.beam_properties[i, 1].TS, job2.beam_properties[i, 1].TS, rtol=1e-5) diff --git a/tests/regression/2_iea22mw/ref_iea22mw_bd_blade.dat b/tests/regression/2_iea22mw/ref_iea22mw_bd_blade.dat index 381331d5..9a98e2b6 100644 --- a/tests/regression/2_iea22mw/ref_iea22mw_bd_blade.dat +++ b/tests/regression/2_iea22mw/ref_iea22mw_bd_blade.dat @@ -2,7 +2,7 @@ Test Format 1 ---------------------- BLADE PARAMETERS -------------------------------------- 1 station_total - Number of blade input stations (-) - 1 damp_type - Damping type: 0: no damping; 1: damped + 1 damp_type - Damping type (switch) 0: none, 1: stiffness-proportional, 2: modal ---------------------- DAMPING COEFFICIENT------------------------------------ mu1 mu2 mu3 mu4 mu5 mu6 (-) (-) (-) (-) (-) (-) @@ -12,16 +12,16 @@ 0.01 0.01 0.01 0.01 0.01 0.01 zeta - Damping coefficients for mode 1 through n_modes ---------------------- DISTRIBUTED PROPERTIES--------------------------------- 0.300000 - 3.4726579277802676e+08 3.2331848070463825e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2136780783975720e+07 - 3.2331848070307989e+07 6.7579167322169662e+08 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 5.5148848109211303e+07 - 0.0000000000000000e+00 0.0000000000000000e+00 2.4860932605673111e+10 1.7969593744627364e+09 -6.3438318210125780e+08 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 1.7969593744625189e+09 4.6172668667864792e+10 -2.3734426349607768e+09 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 -6.3438318210093641e+08 -2.3734426349624014e+09 2.1691154207196560e+10 0.0000000000000000e+00 - 9.2136780784259334e+07 5.5148848107970886e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.5595617642793138e+09 + 3.4686917124906659e+08 3.3671107819947802e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2473416285119236e+07 + 3.3671107819947802e+07 6.8569435462077546e+08 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 5.7106137795098320e+07 + 0.0000000000000000e+00 0.0000000000000000e+00 2.4860932605669075e+10 1.7969593745328412e+09 -6.3438318212379301e+08 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 1.7969593745328412e+09 4.6172668667878128e+10 -2.3734426349940510e+09 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -6.3438318212379301e+08 -2.3734426349940510e+09 2.1691154207189724e+10 0.0000000000000000e+00 + 9.2473416285119236e+07 5.7106137795098320e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.5599396727209070e+09 - 6.4658561585132838e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -4.7147642615256109e+02 - 0.0000000000000000e+00 6.4658561585132838e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.4648154690139727e+01 - 0.0000000000000000e+00 0.0000000000000000e+00 6.4658561585132838e+02 4.7147642615256109e+02 -2.4648154690139727e+01 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 4.7147642615256109e+02 2.7220319590339354e+03 -7.8548437132159407e+01 0.0000000000000000e+00 - 0.0000000000000000e+00 0.0000000000000000e+00 -2.4648154690139727e+01 -7.8548437132159407e+01 3.7861641108638673e+02 0.0000000000000000e+00 - -4.7147642615256109e+02 2.4648154690139727e+01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 3.1006483701203342e+03 + 6.4658561585133305e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -4.7147642615259173e+02 + 0.0000000000000000e+00 6.4658561585133305e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.4648154690140316e+01 + 0.0000000000000000e+00 0.0000000000000000e+00 6.4658561585133305e+02 4.7147642615259173e+02 -2.4648154690140316e+01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 4.7147642615259173e+02 2.7220319590341392e+03 -7.8548437132174215e+01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -2.4648154690140316e+01 -7.8548437132174215e+01 3.7861641108638753e+02 0.0000000000000000e+00 + -4.7147642615259173e+02 2.4648154690140316e+01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 3.1006483701205266e+03 diff --git a/tests/regression/6_box_beam/test_stress_recov.py b/tests/regression/6_box_beam/test_stress_recov.py index 94c291ef..b14e77f5 100644 --- a/tests/regression/6_box_beam/test_stress_recov.py +++ b/tests/regression/6_box_beam/test_stress_recov.py @@ -999,7 +999,7 @@ def test_output_maps(self): npt.assert_allclose(stress, elem_stress, atol=1e-18, err_msg='Stresses are different on loaded recovery.') npt.assert_allclose(strain, elem_strain, atol=1e-7, err_msg='Strains are different on loaded recovery.') - + def test_stress_map_zero_twist(self): """ Test that the stress maps are correctly rotated to stay with the GEBT diff --git a/tests/regression/7_viscoelastic/ref_two_term_bd_blade_visc.dat b/tests/regression/7_viscoelastic/ref_two_term_bd_blade_visc.dat index 33c5486d..fc4b08d7 100644 --- a/tests/regression/7_viscoelastic/ref_two_term_bd_blade_visc.dat +++ b/tests/regression/7_viscoelastic/ref_two_term_bd_blade_visc.dat @@ -7,23 +7,23 @@ 0.02299974 0.14994091 inf ---------------------- DISTRIBUTED PROPERTIES--------------------------------- 0.000000 - 1.9665378649788186e+07 6.3055916065999202e+02 -1.0596446510373938e-25 -5.3669926558797432e-17 -2.6736786846400837e-16 -3.4170045171919944e-01 - 6.3055916071755360e+02 5.6354294669792451e+07 3.7987822508967875e-13 9.6698959240565326e-16 -1.8580671993194649e-14 -8.1546799526975047e+00 - 2.7424783318537356e-14 3.6957008766891517e-14 2.1831356082262930e+08 -2.9321663591980022e-08 5.9039144660167171e+00 7.0795168782312469e-15 - -3.5253919275946727e-15 -2.5925105343049129e-15 -2.1829374998076511e-08 9.6951378136403335e+05 -7.6781394842223527e+00 -5.1373376454589595e-16 - -1.3315433037532005e-15 -1.8019369730409791e-15 5.9039144688736513e+00 -7.6781394840564658e+00 2.8937597537821624e+05 -3.4510814436121627e-16 - -3.3776077839723290e-01 -8.1393723459371863e+00 1.8617318932166703e-14 6.5218805541439177e-17 -8.8429149001247509e-16 2.8183512372423348e+05 + 1.9665378682881441e+07 6.3056114727279851e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -3.4242178407083868e-01 + 6.3056114727869146e+02 5.6354294666494809e+07 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -8.1547402643002442e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 2.1831356082261625e+08 -6.2101162109742521e-05 5.9039054767654111e+00 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -6.2049817636476498e-05 9.6951378136405314e+05 -7.6781402586328262e+00 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 5.9039054775349618e+00 -7.6781402585980221e+00 2.8937597537820227e+05 0.0000000000000000e+00 + -3.3848210974078768e-01 -8.1394326327067770e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.8183512372422894e+05 - 1.9877703884593605e+06 6.3736724832657686e+01 -1.5210739576036056e-26 3.3421358933305893e-17 -6.7127917733124537e-18 -3.4538975912377548e-02 - 6.3736724838417707e+01 5.6962746663580090e+06 -8.5803581086749989e-14 -1.6461601493249896e-14 4.2676123407343990e-15 -8.2427252612755020e-01 - 8.0382623699156603e-17 -6.0539779983968755e-15 2.2067067170711074e+07 -2.9923746593763171e-09 5.9676584731645699e-01 -1.1108563242162073e-15 - 5.0942521786175066e-16 -9.5152479761638393e-17 -2.1771827691696946e-09 9.7998153003753119e+04 -7.7610396305948370e-01 -1.2030352847209010e-17 - -5.6884403445063805e-18 2.9694973772141666e-16 5.9676584751106998e-01 -7.7610396304579643e-01 2.9250034043690212e+04 5.4468262809423448e-17 - -3.4140754953416147e-02 -8.2272523797781327e-01 -4.2598158454429705e-15 -8.0610976760451514e-16 2.1116883601880089e-16 2.8487807092025898e+04 + 1.9877703918044183e+06 6.3736925638921548e+01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -3.4611887978849037e-02 + 6.3736925639507866e+01 5.6962746660247017e+06 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -8.2427862242783623e-01 + 0.0000000000000000e+00 0.0000000000000000e+00 2.2067067170709703e+07 -6.2771344747953918e-06 5.9676493861285407e-01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -6.2719090096427066e-06 9.7998153003755637e+04 -7.7610404134191024e-01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 5.9676493877992709e-01 -7.7610404133681987e-01 2.9250034043688756e+04 0.0000000000000000e+00 + -3.4213666886150775e-02 -8.2273133176721114e-01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.8487807092025505e+04 - 1.9454652646447238e+08 6.2380235143867185e+03 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -3.3803893205160702e+00 - 6.2380235149459631e+03 5.5750425530105424e+08 3.2821558390801757e-12 -9.4470036142524460e-13 -1.5638489920510812e-13 -8.0672977998093472e+01 - -2.3187276108286278e-13 -3.8412385045559106e-13 2.1597420367286029e+09 -2.9677403251737724e-07 5.8406505793673830e+01 -7.3015066675758679e-14 - 3.5850638033526972e-14 2.0901212269692965e-14 -2.2433543331337642e-07 9.5912487566487044e+06 -7.5958637408598292e+01 4.2211644389389901e-15 - 1.1237250733948351e-14 1.8749582309920044e-14 5.8406505817779205e+01 -7.5958637407269663e+01 2.8627514300471852e+06 3.5629806261856335e-15 - -3.3414147493328206e+00 -8.0521542237040322e+01 1.5764465703252059e-13 -4.6180516801102872e-14 -7.5078957846561722e-15 2.7881509597489857e+06 + 1.9454652679185775e+08 6.2380431676490389e+03 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -3.3875253496221931e+00 + 6.2380431676995831e+03 5.5750425526842999e+08 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -8.0673574653308364e+01 + 0.0000000000000000e+00 0.0000000000000000e+00 2.1597420367284622e+09 -6.1436038942895403e-04 5.8406416857231456e+01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 -6.1385360150467683e-04 9.5912487566489261e+06 -7.5958645070409617e+01 0.0000000000000000e+00 + 0.0000000000000000e+00 0.0000000000000000e+00 5.8406416868450989e+01 -7.5958645069743199e+01 2.8627514300470375e+06 0.0000000000000000e+00 + -3.3485507643469017e+00 -8.0522138644288617e+01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.7881509597489364e+06 diff --git a/tests/regression/7_viscoelastic/test_viscoelastic.py b/tests/regression/7_viscoelastic/test_viscoelastic.py index 036e596b..058cea35 100644 --- a/tests/regression/7_viscoelastic/test_viscoelastic.py +++ b/tests/regression/7_viscoelastic/test_viscoelastic.py @@ -1,4 +1,3 @@ -import pdb import os import numpy as np import numpy.testing as npt @@ -407,7 +406,7 @@ def test_ortho_viscoelastic(self): err_msg="First 3 rows of viscoelastic do not add to match the elastic.") npt.assert_allclose(tot_prony[3:], stiff[0][3:], - atol=1e-4*tot_prony[3:].max(), rtol=1e-2, + atol=1e-4*tot_prony[3:].max(), rtol=1e-2, err_msg="Last 3 rows of viscoelastic do not add to match the elastic.") if __name__ == "__main__":