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