You can view and download this file on Github: geometricallyExactBeam2Dtest.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test model for GeometricallyExactBeam2D, cantilever beam with tip force and torque
#
# Model: A 2m long shear deformable beam, located between [0,0,0] and [sqrt(2), sqrt(2), 0], which are 45° relative to the x-axis;
# The beam's height is h = 0.5m and the width is b=0.1m; Young's modulus and density are according to a steel material;
# The beam is fixed at [0,0,0], where displacements and rotation are constrained; a force [F,-F,0] with F=5e8 * h**3 * sqrt(0.5) is applied to the tip of the beam.
#
# Author: Johannes Gerstmayr
# Date: 2021-03-25
#
# 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.
#
# *clean example*
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## import libaries
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
import numpy as np
from math import sin, cos, pi
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()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## set up mbs
SC = exu.SystemContainer()
mbs = SC.AddSystem()
## define overall parameters for model
nElements=16 # number of beam finite elements
L=2 # total length of beam
lElem = L/nElements
E=2.07e11 # Young's modulus of beam element in N/m^2
rho=7850 # density of beam element in kg/m^3
b=0.1 # width of rectangular beam element in m
h=0.5 # height of rectangular beam element in m
A=b*h # cross sectional area of beam element in m^2
I=b*h**3/12 # second moment of area of beam element in m^4
nu = 0.3 # Poisson's ratio
EI = E*I
EA = E*A
rhoA = rho*A
rhoI = rho*I
ks = 10*(1+nu)/(12+11*nu)
G = 7.9615e10 #E/(2*(1+nu))
GA = ks*G*A
fEnd=5e8*h**3 # tip load applied to beam element in N
## create nodes with for loop
nodeList=[]
pRefList=[]
elementList=[]
phi = 0.25*pi #angle of beam relative to x-axis
for i in range(nElements+1):
p1Ref = [cos(phi)*lElem*i,sin(phi)*lElem*i,phi]
ni=mbs.AddNode(Rigid2D(referenceCoordinates = p1Ref, initialCoordinates = [0,0,0],
initialVelocities= [0,0,0]))
nodeList += [ni]
pRefList += [p1Ref[0:2]+[0]]
## create elements:
for i in range(nElements):
oBeam = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [nodeList[i],nodeList[i+1]],
physicsLength=lElem,
physicsMassPerLength=rhoA,
physicsCrossSectionInertia=rhoI,
physicsBendingStiffness=EI,
physicsAxialStiffness=EA,
physicsShearStiffness=GA,
visualization=VObjectBeamGeometricallyExact2D(drawHeight = 0.02*h)
))
elementList+=[oBeam]
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## add ground node, markers and constraints for fixed support
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
n0 = nodeList[0]
mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## add tip force and tip torque
tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeList[-1]))
mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [1*fEnd*sin(phi), -1*fEnd*cos(phi), 0]))
mbs.AddLoad(Torque(markerNumber = tipNodeMarker, loadVector = [0, 0, -5e8]))
## assemble system and check some quantities
mbs.Assemble()
n0 = mbs.GetNodeOutput(0, variableType=exu.OutputVariableType.Position, configuration=exu.ConfigurationType.Reference)
exu.Print("n0=",n0)
p = mbs.GetObjectOutputBody(0, variableType=exu.OutputVariableType.Position, localPosition=[0,0,0], configuration=exu.ConfigurationType.Reference)
exu.Print("p=",p)
## set up simulation settings for dynamic and static solution
simulationSettings = exu.SimulationSettings()
tEnd = 1
steps = 2000
simulationSettings.timeIntegration.numberOfSteps = steps
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
#simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.solutionSettings.writeSolutionToFile = False
#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.staticSolver.newton.maxIterations = 50
simulationSettings.staticSolver.numberOfLoadSteps = 10
## change netwon tolerance for larger number of elements
if nElements > 64:
simulationSettings.staticSolver.newton.relativeTolerance = 2e-8
SC.visualizationSettings.nodes.defaultSize = 0.005
## start graphics and solver
if useGraphics:
exu.StartRenderer()
mbs.WaitForUserToContinue()
uTotal = np.zeros(3)
## test two cases: with and without reference rotations
for case in range(2):
for elem in elementList:
#both cases should give the same result for this case!
mbs.SetObjectParameter(elem, 'includeReferenceRotations', case)
mbs.SolveStatic(simulationSettings)
#mbs.SolveDynamic(simulationSettings) #alternative for dynamic solution
uLast = mbs.GetNodeOutput(nodeList[-1], exu.OutputVariableType.Coordinates)
exu.Print("n =",nElements,", uTip =", uLast[0:2])
uTotal += uLast
uTotal = 0.5*uTotal
## stop graphics and print solution
if useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
exu.Print('solution of geometricallyExactBeam2Dtest=',uTotal[1]) #use y-coordinate
exudynTestGlobals.testError = uTotal[1] - (-2.2115028353806547)
exudynTestGlobals.testResult = uTotal[1]