objectGenericODE2Test.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Test for ObjectGenericODE2 with python user function for linear and rotor dynamics
#           This test represents a FEM model of a rotor, which has elastic supports and rotation is locked
#           We compute eigenmodes, compute the linear response as well as the response of the rotor with gyroscopic terms
#
# Author:   Johannes Gerstmayr 
# Date:     2020-05-13
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.FEM import *

useGraphics = True #without test
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
    from modelUnitTests import exudynTestGlobals #for globally storing test results
    useGraphics = exudynTestGlobals.useGraphics
except:
    class ExudynTestGlobals:
        pass
    exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

SC = exu.SystemContainer()
mbs = SC.AddSystem()

import numpy as np
accumulatedError = 0
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Use FEMinterface to import FEM model and create FFRFreducedOrder object
fem = FEMinterface()
inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory

#useGraphics = False

nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
nNodes = len(nodes)
nODE2 = nNodes*3 #total number of ODE2 coordinates in FEM system; size of M and K 

fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
fem.ScaleStiffnessMatrix(1e-2) #for larger deformations, stiffness is reduced to 1%

#+++++++++++ add elastic supports to fem ==> compute correct eigen frequencies
pLeft = [0,0,0]
pRight = [0,0,0.5]
nLeft = fem.GetNodeAtPoint(pLeft)
nRight = fem.GetNodeAtPoint(pRight)
kJoint = 2e8     #joint stiffness
dJoint = kJoint*0.01  #joint damping

fem.AddElasticSupportAtNode(nLeft, springStiffness=[kJoint,kJoint,kJoint])
fem.AddElasticSupportAtNode(nRight, springStiffness=[kJoint,kJoint,kJoint])

#+++++++++++ compute eigenmodes for comparison
nModes = 8
fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = False)
exu.Print("eigen freq.=", fem.GetEigenFrequenciesHz()[0:6+nModes].round(4)) #mode 0 is rigid body mode (free rotation)!
exu.Print("eigen freq. first mode =", fem.GetEigenFrequenciesHz()[1])       #mode1 with supports: 57.6317863976366Hz;  free-free mode6: sparse: 104.63701326020315, dense: 104.63701326063597
accumulatedError += 1e-2*(fem.GetEigenFrequenciesHz()[1]/57.63178639764625 - 1.)   #check mode (with supports); this is subject to small variations between 32 and 64bit! ==>*1e-2

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create generic object for rotor:
forwardWhirl = True    #test: True; switch this flag to turn on rotordynamics effects
backwardWhirl = False   #test: False; switch this flag to turn on rotordynamics effects
excitationSign = 1
if backwardWhirl: excitationSign = -1

forceVector = [0,-1.,0]*nNodes  #static force vector, optional; add some force in y-direction; could also use node mass to compute force due to weight
fUnbalance = 2000               #fictive force, not depending on frequency

nForce = fem.GetNodeAtPoint([0,0,0.15])#position where unbalance (excitation) force acts
exu.Print("excitation node=", nForce)

fRotX = np.array([0]*nODE2)
fRotX[nForce*3+0] = fUnbalance
fRotY = np.array([0]*nODE2)
fRotY[nForce*3+1] = fUnbalance

#sweep parameters:
t1 = 5 #end time of sweep (t0=0)
f0 = 50 #start frequency
f1 = 250 #terminal frequency
omega1=np.pi*2.*f1

G = fem.GetGyroscopicMatrix(rotationAxis=2, sparse=False) #gyroscopic matrix for rotordynamics effects

useSparseG=True #speed up computation, using sparse matrix in user function
if useSparseG:
    from scipy.sparse import csr_matrix
    G=csr_matrix(G) #convert to sparse matrix

def UFforce(mbs, t, itemIndex, q, q_t):
    #print("UFforce")
    omega = 2.*np.pi*FrequencySweep(t, t1, f0,f1)
    fact = omega/omega1
    force = (fact*SweepCos(t, t1, f0, f1))*fRotY + (excitationSign*fact*SweepSin(t, t1, f0, f1))*fRotX #excitationSign: '+' = forward whirl, '-' = backward whirl

    if forwardWhirl or backwardWhirl:
        #omegaQ_t = omega * np.array(q_t)
        force -= G @ (omega * np.array(q_t)) #negative sign: because term belongs to left-hand-side!!!
        #force -= omega*(G @q_t) #negative sign: because term belongs to left-hand-side!!!
    return force

#add nodes:
nodeList = [] #create node list
for node in fem.GetNodePositionsAsArray():
    nodeList += [mbs.AddNode(NodePoint(referenceCoordinates=node))]

#now add generic body built from FEM model with mass and stiffness matrix (optional damping could be added):
oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers = nodeList, 
                                                massMatrix=fem.GetMassMatrix(sparse=False), 
                                                stiffnessMatrix=fem.GetStiffnessMatrix(sparse=False), 
                                                forceVector=forceVector, forceUserFunction=UFforce,
                                                visualization=VObjectGenericODE2(triangleMesh = fem.GetSurfaceTriangles(), 
                                                                                 color=graphics.color.lightred)))

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add markers and joints
nodeDrawSize = 0.0025 #for joint drawing

nMid = fem.GetNodeAtPoint([0,0,0.25])
nTopMid = fem.GetNodeAtPoint([0., 0.05, 0.25]) #lock rotation of body
exu.Print("nMid=",nMid)
exu.Print("nTopMid=",nTopMid)

nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0], visualization = VNodePointGround(show=False))) #ground node for coordinate constraint
mGroundCoordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action

#add constraint to avoid rotation of body
mTopRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nTopMid, coordinate=0)) #x-coordinate of node at y-max
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundCoordinate, mTopRight]))

oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))

mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))

#++++++++++++++++++++++++++++++++++++++++++
#find nodes at left and right surface:
nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])


lenLeft = len(nodeListLeft)
lenRight = len(nodeListRight)
weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
weightsRight = np.array((1./lenRight)*np.ones(lenRight))

addDampers = True
if addDampers:
    for i in range(3):
        mLeft = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLeft, coordinate=i))
        mRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRight, coordinate=i))
    
        mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGroundCoordinate,mLeft], 
                                             stiffness=0, damping=dJoint))
        if i != 2: #exclude double constraint in z-direction (axis)
            mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGroundCoordinate,mRight],
                                                 stiffness=0, damping=dJoint))
    

addJoint = False #set this flag, if not adding supports to stiffness matrix in fem
if addJoint:

    useSpringDamper = True

    mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2, 
                                                    meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
                                                    weightingFactors=weightsLeft))
    mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2, 
                                                    meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers 
                                                    weightingFactors=weightsRight))

    oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
                                        stiffness=[kJoint,kJoint,kJoint], damping=[dJoint,dJoint,dJoint]))
    oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
                                        stiffness=[kJoint,kJoint,0], damping=[dJoint,dJoint,0]))

                                                    

fileDir = 'solution/'
sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=oGenericODE2, meshNodeNumber=nMid, #meshnode number!
                         storeInternal=True,#fileName=fileDir+'nMidDisplacementLinearTest.txt', 
                         outputVariableType = exu.OutputVariableType.Displacement))

mbs.Assemble()

simulationSettings = exu.SimulationSettings()

SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize

SC.visualizationSettings.bodies.show = True
#SC.visualizationSettings.connectors.show = False

SC.visualizationSettings.bodies.deformationScaleFactor = 10 #use this factor to scale the deformation of modes
if SC.visualizationSettings.bodies.deformationScaleFactor !=1:
    SC.visualizationSettings.nodes.show = False #nodes are not scaled

SC.visualizationSettings.openGL.showFaceEdges = True
SC.visualizationSettings.openGL.showFaces = True

#SC.visualizationSettings.sensors.show = True
#SC.visualizationSettings.sensors.drawSimplified = False
SC.visualizationSettings.sensors.defaultSize = 0.01
#SC.visualizationSettings.markers.drawSimplified = False
#SC.visualizationSettings.markers.show = True
SC.visualizationSettings.markers.defaultSize = 0.01

SC.visualizationSettings.loads.drawSimplified = False

SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component

simulationSettings.solutionSettings.solutionInformation = "ObjectGenericODE2 test"
simulationSettings.solutionSettings.writeSolutionToFile=False

h=1e-3
tEnd = 0.05
#if useGraphics:
#    tEnd = 0.1

simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.solutionWritePeriod = h
simulationSettings.timeIntegration.verboseMode = 1
#simulationSettings.timeIntegration.verboseModeFile = 3
simulationSettings.timeIntegration.newton.useModifiedNewton = True

simulationSettings.solutionSettings.sensorsWritePeriod = h

simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
simulationSettings.displayStatistics = False
simulationSettings.displayComputationTime = False

#create animation:
#simulationSettings.solutionSettings.recordImagesInterval = 0.0002
#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"

#useGraphics = True
if useGraphics:
    exu.StartRenderer()
    if 'lastRenderState' in vars():
        SC.SetRenderState(lastRenderState) #load last model view

    #mbs.WaitForUserToContinue() #press space to continue


mbs.SolveDynamic(simulationSettings)
    
if useGraphics:
    SC.WaitForRenderEngineStopFlag()
    exu.StopRenderer() #safely close rendering window!
    lastRenderState = SC.GetRenderState() #store model view for next simulation

accumulatedError += mbs.GetNodeOutput(nMid,exu.OutputVariableType.Position)[0] #take x-coordinate of position

exu.Print('solution of ObjectGenericODE2=',accumulatedError)

exudynTestGlobals.testError = accumulatedError - (-2.2737401292182432e-05) #2020-05-18: -2.2737401292182432e-05 
exudynTestGlobals.testResult = accumulatedError

##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
#plot results
if useGraphics:
    
    
    mbs.PlotSensor(sDisp, components=1, closeAll=True, labels=['uMid,linear'])