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