geometricallyExactBeam2Dtest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for GeometricallyExactBeam2D, cantilever beam with tip force and torque
  5#
  6# 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;
  7#           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;
  8#           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.
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2021-03-25
 12#
 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.
 14#
 15# *clean example*
 16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 17
 18## import libaries
 19import exudyn as exu
 20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 21import exudyn.graphics as graphics #only import if it does not conflict
 22
 23import numpy as np
 24from math import sin, cos, pi
 25
 26useGraphics = True #without test
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 29try: #only if called from test suite
 30    from modelUnitTests import exudynTestGlobals #for globally storing test results
 31    useGraphics = exudynTestGlobals.useGraphics
 32except:
 33    class ExudynTestGlobals:
 34        pass
 35    exudynTestGlobals = ExudynTestGlobals()
 36#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 37
 38## set up mbs
 39SC = exu.SystemContainer()
 40mbs = SC.AddSystem()
 41
 42
 43## define overall parameters for model
 44nElements=16           # number of beam finite elements
 45L=2                    # total length of beam
 46lElem = L/nElements
 47E=2.07e11              # Young's modulus of beam element in N/m^2
 48rho=7850               # density of beam element in kg/m^3
 49b=0.1                  # width of rectangular beam element in m
 50h=0.5                  # height of rectangular beam element in m
 51A=b*h                  # cross sectional area of beam element in m^2
 52I=b*h**3/12            # second moment of area of beam element in m^4
 53nu = 0.3               # Poisson's ratio
 54
 55EI = E*I
 56EA = E*A
 57rhoA = rho*A
 58rhoI = rho*I
 59ks = 10*(1+nu)/(12+11*nu)
 60G = 7.9615e10           #E/(2*(1+nu))
 61GA = ks*G*A
 62fEnd=5e8*h**3        # tip load applied to beam element in N
 63
 64
 65## create nodes with for loop
 66nodeList=[]
 67pRefList=[]
 68elementList=[]
 69phi = 0.25*pi #angle of beam relative to x-axis
 70
 71for i in range(nElements+1):
 72    p1Ref = [cos(phi)*lElem*i,sin(phi)*lElem*i,phi]
 73    ni=mbs.AddNode(Rigid2D(referenceCoordinates = p1Ref, initialCoordinates = [0,0,0],
 74                         initialVelocities= [0,0,0]))
 75    nodeList += [ni]
 76    pRefList += [p1Ref[0:2]+[0]]
 77
 78
 79## create elements:
 80for i in range(nElements):
 81    oBeam = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [nodeList[i],nodeList[i+1]],
 82                                                            physicsLength=lElem,
 83                                                            physicsMassPerLength=rhoA,
 84                                                            physicsCrossSectionInertia=rhoI,
 85                                                            physicsBendingStiffness=EI,
 86                                                            physicsAxialStiffness=EA,
 87                                                            physicsShearStiffness=GA,
 88                                                            visualization=VObjectBeamGeometricallyExact2D(drawHeight = 0.02*h)
 89                                                ))
 90    elementList+=[oBeam]
 91
 92
 93#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 94## add ground node, markers and constraints for fixed support
 95nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 96mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
 97
 98n0 = nodeList[0]
 99mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
100mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
101mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
102mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
103mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
104mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))
105
106#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107## add tip force and tip torque
108tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeList[-1]))
109mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [1*fEnd*sin(phi), -1*fEnd*cos(phi), 0]))
110mbs.AddLoad(Torque(markerNumber = tipNodeMarker, loadVector = [0, 0, -5e8]))
111
112
113## assemble system and check some quantities
114mbs.Assemble()
115n0 = mbs.GetNodeOutput(0, variableType=exu.OutputVariableType.Position, configuration=exu.ConfigurationType.Reference)
116exu.Print("n0=",n0)
117p = mbs.GetObjectOutputBody(0, variableType=exu.OutputVariableType.Position, localPosition=[0,0,0], configuration=exu.ConfigurationType.Reference)
118exu.Print("p=",p)
119
120## set up simulation settings for dynamic and static solution
121simulationSettings = exu.SimulationSettings()
122
123tEnd = 1
124steps = 2000
125simulationSettings.timeIntegration.numberOfSteps = steps
126simulationSettings.timeIntegration.endTime = tEnd
127simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
128#simulationSettings.timeIntegration.verboseMode = 1
129simulationSettings.solutionSettings.writeSolutionToFile = False
130
131#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
132simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
133
134simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
135simulationSettings.timeIntegration.newton.useModifiedNewton = True
136
137
138simulationSettings.staticSolver.newton.maxIterations = 50
139simulationSettings.staticSolver.numberOfLoadSteps = 10
140
141## change netwon tolerance for larger number of elements
142if nElements > 64:
143    simulationSettings.staticSolver.newton.relativeTolerance = 2e-8
144
145SC.visualizationSettings.nodes.defaultSize = 0.005
146
147## start graphics and solver
148if useGraphics:
149    exu.StartRenderer()
150    mbs.WaitForUserToContinue()
151
152uTotal = np.zeros(3)
153
154## test two cases: with and without reference rotations
155for case in range(2):
156    for elem in elementList:
157        #both cases should give the same result for this case!
158        mbs.SetObjectParameter(elem, 'includeReferenceRotations', case)
159
160    mbs.SolveStatic(simulationSettings)
161    #mbs.SolveDynamic(simulationSettings) #alternative for dynamic solution
162
163    uLast = mbs.GetNodeOutput(nodeList[-1], exu.OutputVariableType.Coordinates)
164    exu.Print("n =",nElements,", uTip =", uLast[0:2])
165    uTotal += uLast
166
167uTotal = 0.5*uTotal
168
169## stop graphics and print solution
170if useGraphics:
171    SC.WaitForRenderEngineStopFlag()
172    exu.StopRenderer() #safely close rendering window!
173
174exu.Print('solution of geometricallyExactBeam2Dtest=',uTotal[1]) #use y-coordinate
175
176exudynTestGlobals.testError = uTotal[1] - (-2.2115028353806547)
177exudynTestGlobals.testResult = uTotal[1]