compareAbaqusAnsysRotorEigenfrequencies.py

You can view and download this file on Github: compareAbaqusAnsysRotorEigenfrequencies.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example to compute eigenfrequencies of a rotor
  5#           NOTE that this example requires files from a subdirectory testData as provided on github
  6#
  7# Author:   Stefan Holzinver, Johannes Gerstmayr
  8# Date:     2020-05-18
  9#
 10# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import numpy as np
 15from scipy.sparse import linalg
 16import scipy as sp
 17
 18import exudyn as exu
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics
 21from exudyn.FEM import *
 22
 23testIsActive = exu.sys.get('testIsActive', False)
 24
 25numberOfModes = 18
 26useSparseSolverRoutine = False
 27
 28errorResult = 0
 29
 30testDataDir = "testData/"
 31#if not testIsActive:
 32#    testDataDir = "testData/"
 33
 34###############################################################################
 35# Ansys - lumped mass matrix formulation - Sparse Matrix - MMF format
 36###############################################################################
 37
 38#read finite element model
 39exu.sys['testResult'] = 0
 40for testNumber in range(2):
 41    if testNumber == 1:
 42        useSparseSolverRoutine = True
 43    
 44    fem = FEMinterface()
 45    #exu.Print("read matrices ...")
 46    fem.ReadMassMatrixFromAnsys(fileName=testDataDir + 'rotorAnsysMassMatrixSparse.txt', 
 47                                dofMappingVectorFile=testDataDir + 'rotorAnsysMatrixDofMappingVector.txt')
 48    fem.ReadStiffnessMatrixFromAnsys(fileName=testDataDir + 'rotorAnsysStiffnessMatrixSparse.txt',
 49                                dofMappingVectorFile=testDataDir + 'rotorAnsysMatrixDofMappingVector.txt')
 50    
 51    #exu.Print("compute eigenmodes ...")
 52    fem.ComputeEigenmodes(numberOfModes, useSparseSolver = useSparseSolverRoutine)
 53    
 54    if not testIsActive:
 55        exu.Print('natural frequencies from Ansys model (Lumped Mass Matrix, MMF-Format)', fem.GetEigenFrequenciesHz()[0:numberOfModes])
 56    
 57    if not useSparseSolverRoutine:
 58        f6 = fem.GetEigenFrequenciesHz()[6] - 104.63701055079783 #reference solution
 59    else:
 60        #compare with dense matrix solution; 
 61        #due to random initialization, the results change with every computation ==> approx. 1e-11 repeatability
 62        f6 = fem.GetEigenFrequenciesHz()[6] - 104.6370105507865 
 63    f6*=1e-6 #use offset also for abaqus, as it gives non-reproducible results in dense case (32/64bit?)
 64    exu.Print('natural frequencies from Ansys model, sparse=',str(useSparseSolverRoutine),":", fem.GetEigenFrequenciesHz()[6] )
 65    errorResult += f6
 66
 67    exu.sys['testResult'] += 1e-6*fem.GetEigenFrequenciesHz()[6]
 68
 69    ###############################################################################
 70    # Abaqus
 71    ###############################################################################
 72    
 73    #read finite element model
 74    fem = FEMinterface()
 75    #exu.Print("read matrices ...")
 76    fem.ReadMassMatrixFromAbaqus(fileName=testDataDir + 'rotorDiscTestMASS1.mtx')
 77    fem.ReadStiffnessMatrixFromAbaqus(fileName=testDataDir + 'rotorDiscTestSTIF1.mtx')
 78    fem.ScaleStiffnessMatrix(1e-2) #in Ansys, the stiffness matrix is already scaled!
 79        
 80    #exu.Print("compute eigenmodes ...")
 81    fem.ComputeEigenmodes(numberOfModes, useSparseSolver = useSparseSolverRoutine)
 82
 83    if not testIsActive:
 84        exu.Print('natural frequencies from Abaqus model (Lumped Mass Matrix)',fem.GetEigenFrequenciesHz()[0:numberOfModes])
 85    
 86    if not useSparseSolverRoutine:
 87        f6 = fem.GetEigenFrequenciesHz()[6] - 104.6370132606291 #reference solution
 88    else:
 89        #compare with dense matrix solution; 
 90        #due to random initialization, the results change with every computation ==> approx. 1e-11 repeatability
 91        f6 = fem.GetEigenFrequenciesHz()[6] - 104.6370132606291 
 92    f6*=1e-6 #use offset also for abaqus, as it gives non-reproducible results in dense case (32/64bit?)
 93    exu.Print('natural frequencies from Abaqus model, sparse=',str(useSparseSolverRoutine),":", fem.GetEigenFrequenciesHz()[6] )
 94    errorResult += f6
 95    exu.sys['testResult'] += 1e-6*fem.GetEigenFrequenciesHz()[6]
 96
 97exu.Print('error of compareAbaqusAnsysRotorEigenfrequencies (due to sparse solver)=',errorResult)
 98if abs(errorResult) < 1e-15: #usually of size 1e-17
 99    errorResult = 0 #due to randomized sparse solver results, take this treshold!
100    
101exu.Print('solution of compareAbaqusAnsysRotorEigenfrequencies (with treshold)=',errorResult)
102#exu.sys['testResult'] computed above