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]