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