gridGeomExactBeam2D.py
You can view and download this file on Github: gridGeomExactBeam2D.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test model for GeometricallyExactBeam2D, evaluating a grid of beams
5#
6# Model: Planar model of beams arranged at grid (horizontal and vertical lines), rigidly connected;
7# The grid has a size of 16 elements in x-direction and 4 elements in y-direction.
8#
9# Author: Johannes Gerstmayr
10# Date: 2021-03-25
11#
12# 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.
13#
14# *clean example*
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17## import libraries
18import exudyn as exu
19from exudyn.utilities import *
20
21import numpy as np
22from math import sin, cos, pi
23
24useGraphics = True #without test
25#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
27try: #only if called from test suite
28 from modelUnitTests import exudynTestGlobals #for globally storing test results
29 useGraphics = exudynTestGlobals.useGraphics
30except:
31 class ExudynTestGlobals:
32 pass
33 exudynTestGlobals = ExudynTestGlobals()
34#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35
36## set up system and define parameters
37SC = exu.SystemContainer()
38mbs = SC.AddSystem()
39
40lElem = 0.5 # length of one finite element
41lElemY = lElem*1
42E=2.1e11*0.1 # Steel; Young's modulus of beam element in N/m^2
43rho=7800 # Steel; density of beam element in kg/m^3
44b=0.1 # width of rectangular beam element in m
45h=0.1 # height of rectangular beam element in m
46A=b*h # cross sectional area of beam element in m^2
47I=b*h**3/12 # second moment of area of beam element in m^4
48nu = 0.3 # Poisson's ratio for steel
49
50EI = E*I
51EA = E*A
52rhoA = rho*A
53rhoI = rho*I
54ks = 10*(1+nu)/(12+11*nu) # shear correction factor
55G = E/(2*(1+nu)) # shear modulus
56GA = ks*G*A # shear stiffness of beam
57
58nX = 16 #grid size x
59nY = 4 #grid size y
60
61## create nodes
62nodeIndices = np.zeros((nX, nY), dtype=int)
63
64for j in range(nY):
65 for i in range(nX):
66 pRef = [i*lElem,j*lElemY, 0] #zero angle; reference rotation not used!
67
68 ni=mbs.AddNode(NodeRigidBody2D(referenceCoordinates = pRef,
69 initialCoordinates = [0,0,0],
70 initialVelocities= [0,0,0]))
71 nodeIndices[i,j]=int(ni)
72
73## create elements in x-direction
74for j in range(nY):
75 for i in range(nX-1):
76 n0 = nodeIndices[i, j]
77 n1 = nodeIndices[i+1, j]
78 oGeneric = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [n0,n1],
79 physicsLength=lElem,
80 physicsMassPerLength=rhoA,
81 physicsCrossSectionInertia=rhoI,
82 physicsBendingStiffness=EI,
83 physicsAxialStiffness=EA,
84 physicsShearStiffness=GA,
85 includeReferenceRotations=False,
86 visualization=VObjectBeamGeometricallyExact2D(drawHeight = h)
87 ))
88
89## create elements in y-direction
90for j in range(nY-1):
91 for i in range(int(nX/4)):
92 n0 = nodeIndices[(i+1)*4-1, j]
93 n1 = nodeIndices[(i+1)*4-1, j+1]
94 mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [n0,n1],
95 physicsLength=lElemY,
96 physicsMassPerLength=rhoA,
97 physicsCrossSectionInertia=rhoI,
98 physicsBendingStiffness=EI,
99 physicsAxialStiffness=EA,
100 physicsShearStiffness=GA,
101 includeReferenceRotations=False,
102 visualization=VObjectBeamGeometricallyExact2D(drawHeight = h) ))
103
104
105
106#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107## create ground node and marker coordinate attached to ground node
108nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
109mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
110
111## add markers and constraints for fixed supports
112for j in range(nY):
113 n0 = nodeIndices[0,j]
114 mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
115 mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
116 mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
117 mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
118 mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
119 mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))
120
121 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122 #add tip load
123 tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeIndices[-1,j]))
124 mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [0, -1e5, 0]))
125
126
127## assemble
128mbs.Assemble()
129
130## set up simulation settings
131simulationSettings = exu.SimulationSettings()
132
133tEnd = 0.1
134steps = 100
135simulationSettings.timeIntegration.numberOfSteps = steps
136simulationSettings.timeIntegration.endTime = tEnd
137simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
138simulationSettings.timeIntegration.verboseMode = 1
139simulationSettings.solutionSettings.writeSolutionToFile = False
140#simulationSettings.timeIntegration.simulateInRealtime = True
141#simulationSettings.timeIntegration.realtimeFactor = 0.1
142
143#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
144simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
145
146simulationSettings.timeIntegration.newton.useModifiedNewton = True
147
148
149simulationSettings.staticSolver.newton.maxIterations = 50
150simulationSettings.staticSolver.numberOfLoadSteps = 10
151# simulationSettings.displayComputationTime = True
152# simulationSettings.displayStatistics = True
153
154
155SC.visualizationSettings.nodes.defaultSize = 0.005
156# SC.visualizationSettings.bodies.beams.crossSectionFilled = False
157SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.ForceLocal
158SC.visualizationSettings.contour.outputVariableComponent = 0
159
160## start graphics
161if useGraphics:
162 exu.StartRenderer()
163 mbs.WaitForUserToContinue()
164
165## start dynamic solver
166mbs.SolveDynamic(simulationSettings)
167
168## stop graphics
169if useGraphics:
170 SC.WaitForRenderEngineStopFlag()
171 exu.StopRenderer() #safely close rendering window!
172
173## read and print solution
174uLast = mbs.GetNodeOutput(nodeIndices[-1,-1], exu.OutputVariableType.Coordinates)
175exu.Print("grid =",nodeIndices.shape,", uTip =", uLast[0:2])
176
177exu.Print('solution of gridGeomExactBeam2D=',uLast[1]) #use y-coordinate
178
179exudynTestGlobals.testError = uLast[1] - (-2.2115028353806547)
180exudynTestGlobals.testResult = uLast[1]