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 * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics #only import if it does not conflict
 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 system and define parameters
 38SC = exu.SystemContainer()
 39mbs = SC.AddSystem()
 40
 41lElem = 0.5            # length of one finite element
 42lElemY = lElem*1
 43E=2.1e11*0.1               # Steel; Young's modulus of beam element in N/m^2
 44rho=7800               # Steel; density of beam element in kg/m^3
 45b=0.1                  # width of rectangular beam element in m
 46h=0.1                 # height of rectangular beam element in m
 47A=b*h                  # cross sectional area of beam element in m^2
 48I=b*h**3/12            # second moment of area of beam element in m^4
 49nu = 0.3               # Poisson's ratio for steel
 50
 51EI = E*I
 52EA = E*A
 53rhoA = rho*A
 54rhoI = rho*I
 55ks = 10*(1+nu)/(12+11*nu) # shear correction factor
 56G = E/(2*(1+nu))          # shear modulus
 57GA = ks*G*A               # shear stiffness of beam
 58
 59nX = 16             #grid size x
 60nY = 4              #grid size y
 61
 62## create nodes
 63nodeIndices = np.zeros((nX, nY), dtype=int)
 64
 65for j in range(nY):
 66    for i in range(nX):
 67        pRef = [i*lElem,j*lElemY, 0] #zero angle; reference rotation not used!
 68
 69        ni=mbs.AddNode(NodeRigidBody2D(referenceCoordinates = pRef,
 70                                       initialCoordinates = [0,0,0],
 71                                       initialVelocities= [0,0,0]))
 72        nodeIndices[i,j]=int(ni)
 73
 74## create elements in x-direction
 75for j in range(nY):
 76    for i in range(nX-1):
 77        n0 = nodeIndices[i, j]
 78        n1 = nodeIndices[i+1, j]
 79        oGeneric = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [n0,n1],
 80                                                                physicsLength=lElem,
 81                                                                physicsMassPerLength=rhoA,
 82                                                                physicsCrossSectionInertia=rhoI,
 83                                                                physicsBendingStiffness=EI,
 84                                                                physicsAxialStiffness=EA,
 85                                                                physicsShearStiffness=GA,
 86                                                                includeReferenceRotations=False,
 87                                                                visualization=VObjectBeamGeometricallyExact2D(drawHeight = h)
 88                                                    ))
 89
 90## create elements in y-direction
 91for j in range(nY-1):
 92    for i in range(int(nX/4)):
 93        n0 = nodeIndices[(i+1)*4-1, j]
 94        n1 = nodeIndices[(i+1)*4-1, j+1]
 95        mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [n0,n1],
 96                                                     physicsLength=lElemY,
 97                                                     physicsMassPerLength=rhoA,
 98                                                     physicsCrossSectionInertia=rhoI,
 99                                                     physicsBendingStiffness=EI,
100                                                     physicsAxialStiffness=EA,
101                                                     physicsShearStiffness=GA,
102                                                     includeReferenceRotations=False,
103                                                     visualization=VObjectBeamGeometricallyExact2D(drawHeight = h) ))
104
105
106
107#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
108## create ground node and marker coordinate attached to ground node
109nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
110mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
111
112## add markers and constraints for fixed supports
113for j in range(nY):
114    n0 = nodeIndices[0,j]
115    mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
116    mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
117    mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
118    mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
119    mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
120    mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))
121
122    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123    #add tip load
124    tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeIndices[-1,j]))
125    mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [0, -1e5, 0]))
126
127
128## assemble
129mbs.Assemble()
130
131## set up simulation settings
132simulationSettings = exu.SimulationSettings()
133
134tEnd = 0.1
135steps = 100
136simulationSettings.timeIntegration.numberOfSteps = steps
137simulationSettings.timeIntegration.endTime = tEnd
138simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
139simulationSettings.timeIntegration.verboseMode = 1
140simulationSettings.solutionSettings.writeSolutionToFile = False
141#simulationSettings.timeIntegration.simulateInRealtime = True
142#simulationSettings.timeIntegration.realtimeFactor = 0.1
143
144#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
145simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
146
147simulationSettings.timeIntegration.newton.useModifiedNewton = True
148
149
150simulationSettings.staticSolver.newton.maxIterations = 50
151simulationSettings.staticSolver.numberOfLoadSteps = 10
152# simulationSettings.displayComputationTime = True
153# simulationSettings.displayStatistics = True
154
155
156SC.visualizationSettings.nodes.defaultSize = 0.005
157# SC.visualizationSettings.bodies.beams.crossSectionFilled = False
158SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.ForceLocal
159SC.visualizationSettings.contour.outputVariableComponent = 0
160
161## start graphics
162if useGraphics:
163    exu.StartRenderer()
164    mbs.WaitForUserToContinue()
165
166## start dynamic solver
167mbs.SolveDynamic(simulationSettings)
168
169## stop graphics
170if useGraphics:
171    SC.WaitForRenderEngineStopFlag()
172    exu.StopRenderer() #safely close rendering window!
173
174## read and print solution
175uLast = mbs.GetNodeOutput(nodeIndices[-1,-1], exu.OutputVariableType.Coordinates)
176exu.Print("grid =",nodeIndices.shape,", uTip =", uLast[0:2])
177
178exu.Print('solution of gridGeomExactBeam2D=',uLast[1]) #use y-coordinate
179
180exudynTestGlobals.testError = uLast[1] - (-2.2115028353806547)
181exudynTestGlobals.testResult = uLast[1]