You can view and download this file on Github: pendulumVerify.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
#
# Author: Johannes Gerstmayr
# Date: 2021-04-20
#
# 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.itemInterface import *
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.FEM import *
from exudyn.graphicsDataUtilities import *
SC = exu.SystemContainer()
mbs = SC.AddSystem()
import numpy as np
import time
import exudyn.basicUtilities as eb
import exudyn.rigidBodyUtilities as rb
import exudyn.utilities as eu
import numpy as np
useGraphics = True
fileName = 'testData/pendulumSimple' #for load/save of FEM data
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#netgen/meshing part:
femInterface = FEMinterface()
#geometrical parameters:
L = 400 #Length of plate (X)
ff=1 #factor for higher thickness
h = 8*ff #height of plate (Y)
w = 4*ff #width of plate (Z)
nModes = 10
meshH = 2*ff
print("mesh h=", meshH)
#steel:
rho = 2.7e-9 #tons/mm^3
Emodulus=70e3 #N/mm^2
g = 9810 #[mm/s^2]
nu=0.35
gForce=[0,-g,0]
V = L*h*w
mass = V*rho
print("V=", V, " (mm^3), mass=", mass*1000, "(kg)")
useGravity = True
A=w*h
q=rho*A*g
EI = Emodulus*w*h**3/12.
F=1 #N
if useGravity:
tipDisp = q*L**4/(8*EI)
print("tip displ (weight)=", tipDisp)
else:
tipDisp = F*L**3/(3*EI)
print("tip displ (tip F )=", tipDisp)
#h*w=8*4, nu=0.35, E=70e3:
#F=1
#meshH=2: w = 1.5989535
#meshH=1: w = 1.73735541
#meshH=0.5: w = 1.77126479
#analytical:w = 1.78571428
#weight:
#meshH=2: w = 0.20311787
#meshH=1: w = 0.220752717
#analytical:w = 0.2270314285714286
meshCreated = False
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
import ngsolve as ngs
import netgen
from netgen.meshing import *
from netgen.geom2d import unit_square
#import netgen.libngpy as libng
from netgen.csg import *
geo = CSGeometry()
#plate
block = OrthoBrick(Pnt(0,-0.5*h, -0.5*w),Pnt(L, 0.5*h, 0.5*w))
geo.Add(block)
mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH))
mesh.Curve(1)
if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
# import netgen
import netgen.gui
ngs.Draw(mesh)
netgen.Redraw()
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#Use femInterface to import femInterface model and create FFRFreducedOrder object
eigenModesComputed = False
femInterface.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
computeEigenmodes=eigenModesComputed, verbose=False, excludeRigidBodyModes = 6,
numberOfModes = nModes, maxEigensolveIterations=20)
P1=[0,0,0]
nodesLeft = femInterface.GetNodesInPlane(P1, [1,0,0])
#print("boundary nodes bolt=", nodesLeft)
nodesLeftLen = len(nodesLeft)
nodesLeftWeights = np.array((1./nodesLeftLen)*np.ones(nodesLeftLen))
P2=[L,0,0]
nodesRight = femInterface.GetNodesInPlane(P2, [1,0,0])
#print("boundary nodes bolt=", nodesRight)
nodesRightLen = len(nodesRight)
nodesRightWeights = np.array((1./nodesRightLen)*np.ones(nodesRightLen))
print("nNodes=",femInterface.NumberOfNodes())
strMode = ''
boundaryList = [nodesLeft, nodesRight]
print("compute HCB modes... ")
start_time = time.time()
femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
nEigenModes=nModes,
useSparseSolver=True,
computationMode = HCBstaticModeSelection.RBE2)
print("eigen freq.=", femInterface.GetEigenFrequenciesHz())
print("HCB modes needed %.3f seconds" % (time.time() - start_time))
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
if False:
mat = KirchhoffMaterial(Emodulus, nu, rho)
varType = exu.OutputVariableType.StressLocal
#varType = exu.OutputVariableType.StrainLocal
print("ComputePostProcessingModes ... (may take a while)")
start_time = time.time()
femInterface.ComputePostProcessingModes(material=mat,
outputVariableType=varType)
print(" ... needed %.3f seconds" % (time.time() - start_time))
SC.visualizationSettings.contour.reduceRange=False
SC.visualizationSettings.contour.outputVariable = varType
SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
else:
SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
SC.visualizationSettings.contour.outputVariableComponent = 1
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
print("create CMS element ...")
cms = ObjectFFRFreducedOrderInterface(femInterface)
objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
initialVelocity=[0,0,0],
initialAngularVelocity=[0,0,0],
color=[0.9,0.9,0.9,1.],
)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#add markers and joints
nodeDrawSize = 1 #for joint drawing
#mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
if True:
oGround = mbs.AddObject(ObjectGround(referencePosition= P1))
#compute offset of nodes at plane (because the average does not give [0,0,0]):
pOff = np.zeros(3)
nodes = femInterface.nodes['Position']
for i in nodesLeft:
pOff += nodes[i]
pOff *= 1/len(nodesLeft)
#create marker:
altApproach = True
mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
meshNodeNumbers=np.array(nodesLeft), #these are the meshNodeNumbers
offset=-pOff,
useAlternativeApproach=altApproach,
weightingFactors=nodesLeftWeights))
lockedAxes=[1,1,1,1,1,1]
if True:
mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
localPosition=P1,
visualization=VMarkerBodyRigid(show=True)))
mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
constrainedAxes = lockedAxes,
visualization=VGenericJoint(show=False, axesRadius=0.1*w, axesLength=0.1*h)))
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#animate modes
SC.visualizationSettings.markers.show = True
SC.visualizationSettings.markers.defaultSize=1
SC.visualizationSettings.markers.drawSimplified = False
SC.visualizationSettings.loads.show = False
SC.visualizationSettings.loads.drawSimplified = False
SC.visualizationSettings.loads.defaultSize=10
SC.visualizationSettings.loads.defaultRadius = 0.1
SC.visualizationSettings.openGL.multiSampling=4
SC.visualizationSettings.openGL.lineWidth=2
if False: #activate to animate modes
from exudyn.interactive import AnimateModes
mbs.Assemble()
SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.openGL.showFaceEdges = True
SC.visualizationSettings.openGL.multiSampling=4
SC.visualizationSettings.openGL.lineWidth=2
SC.visualizationSettings.window.renderWindowSize = [1600,1080]
SC.visualizationSettings.contour.showColorBar = False
SC.visualizationSettings.general.textSize = 16
#%%+++++++++++++++++++++++++++++++++++++++
#animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n')
import sys
sys.exit()
oFFRF = objFFRF['oFFRFreducedOrder']
if useGravity:
#add gravity (not necessary if user functions used)
mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= gForce))
else:
mRight = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
meshNodeNumbers=np.array(nodesRight), #these are the meshNodeNumbers
useAlternativeApproach=altApproach,
weightingFactors=nodesRightWeights))
mbs.AddLoad(LoadForceVector(markerNumber=mRight, loadVector=[0,-F,0]))
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
fileDir = 'solution/'
# sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
# fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
# outputVariableType = exu.OutputVariableType.Position))
# sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
# fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
# outputVariableType = exu.OutputVariableType.Position))
nTip = femInterface.GetNodeAtPoint([L,0.5*h,0.5*w])
sensTip = mbs.AddSensor(SensorSuperElement(bodyNumber=oFFRF,
meshNodeNumber=nTip,
fileName=fileDir+'displacementTip.txt',
outputVariableType = exu.OutputVariableType.DisplacementLocal))
# print("tip0=",mbs.GetSensorValues(sensTip))
mbs.Assemble()
simulationSettings = exu.SimulationSettings()
SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
SC.visualizationSettings.nodes.basisSize = 0.12
SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
SC.visualizationSettings.openGL.showFaceEdges = True
SC.visualizationSettings.openGL.showFaces = True
SC.visualizationSettings.sensors.show = True
SC.visualizationSettings.sensors.drawSimplified = False
SC.visualizationSettings.sensors.defaultSize = 2
SC.visualizationSettings.loads.defaultSize = 10
simulationSettings.solutionSettings.solutionInformation = "pendulum verification"
h=0.25e-3*4
tEnd = 0.25*8
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.writeSolutionToFile = False
simulationSettings.timeIntegration.verboseMode = 1
#simulationSettings.timeIntegration.verboseModeFile = 3
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.solutionSettings.sensorsWritePeriod = h
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
#simulationSettings.displayStatistics = True
#simulationSettings.displayComputationTime = True
#create animation:
# simulationSettings.solutionSettings.recordImagesInterval = 0.005
# SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
SC.visualizationSettings.window.renderWindowSize=[1920,1080]
SC.visualizationSettings.openGL.multiSampling = 4
useGraphics=False
if useGraphics:
if useGraphics:
SC.visualizationSettings.general.autoFitScene=False
exu.StartRenderer()
if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
mbs.WaitForUserToContinue() #press space to continue
#SC.RedrawAndSaveImage()
if False:
mbs.SolveDynamic(simulationSettings=simulationSettings)
else:
mbs.SolveStatic(simulationSettings=simulationSettings)
# print("tip1=",mbs.GetSensorValues(sensTip))
if useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
data = np.loadtxt('solution/displacementTip.txt', comments='#', delimiter=',')
print("tip disp=", data[-1,1:])