LShapeGeomExactBeam2D.py

You can view and download this file on Github: LShapeGeomExactBeam2D.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for GeometricallyExactBeam2D, evaluating behavior of rotated beam
  5#
  6# Model:    Planar model of L-shaped beams with tip load; the beam has 4 elements at each part of the L-shape:
  7#           the element length is 0.25m, the material is steel and height is 0.05m and with 0.1m;
  8#           the L-shape is fixed at the left end to ground and a tip load of [1e6,0,0] is applied to the tip.
  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 #only import if it does not conflict
 22
 23import numpy as np
 24from math import sin, cos, pi
 25
 26useGraphics = True #without test
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 29try: #only if called from test suite
 30    from modelUnitTests import exudynTestGlobals #for globally storing test results
 31    useGraphics = exudynTestGlobals.useGraphics
 32except:
 33    class ExudynTestGlobals:
 34        pass
 35    exudynTestGlobals = ExudynTestGlobals()
 36#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 37
 38## setup system container and mbs
 39SC = exu.SystemContainer()
 40mbs = SC.AddSystem()
 41
 42## define parameters for beams
 43a = 0.25
 44lElem = a              # length of one finite element
 45E=2.1e11               # Steel; Young's modulus of beam element in N/m^2
 46rho=7800               # Steel; density of beam element in kg/m^3
 47b=0.1                  # width of rectangular beam element in m
 48h=0.05                 # height of rectangular beam element in m
 49A=b*h                  # cross sectional area of beam element in m^2
 50I=b*h**3/12            # second moment of area of beam element in m^4
 51nu = 0.3               # Poisson's ratio for steel
 52
 53EI = E*I
 54EA = E*A
 55rhoA = rho*A
 56rhoI = rho*I
 57ks = 10*(1+nu)/(12+11*nu) # shear correction factor
 58G = E/(2*(1+nu))          # shear modulus
 59GA = ks*G*A               # shear stiffness of beam
 60
 61
 62## create position list for nodes:
 63nodeList=[]
 64pRefList=[[0,0],
 65          [1*a,0],
 66          [2*a,0],
 67          [3*a,0],
 68          [4*a,0],
 69          [4*a,1*a],
 70          [4*a,2*a],
 71          [4*a,3*a],
 72          [4*a,4*a],
 73          ]
 74
 75## create nodes
 76for p in pRefList:
 77    pRef = [p[0],p[1],0] #always angle zero, which is not considered in beam for includeReferenceRotations=False
 78    ni=mbs.AddNode(NodeRigidBody2D(referenceCoordinates = pRef,
 79                                   initialCoordinates = [0,0,0],
 80                                   initialVelocities= [0,0,0]))
 81    nodeList += [ni]
 82
 83## create beam elements:
 84for i in range(len(nodeList)-1):
 85    mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [nodeList[i],nodeList[i+1]],
 86                                                 physicsLength=lElem,
 87                                                 physicsMassPerLength=rhoA,
 88                                                 physicsCrossSectionInertia=rhoI,
 89                                                 physicsBendingStiffness=EI,
 90                                                 physicsAxialStiffness=EA,
 91                                                 physicsShearStiffness=GA,
 92                                                 includeReferenceRotations=False, #to connect beams at 90° at same node
 93                                                 visualization=VObjectBeamGeometricallyExact2D(drawHeight = h) ))
 94
 95
 96#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 97## create ground node with marker for coordinate constraints
 98nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 99mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
100
101## add markers and constraints for fixed support
102n0 = nodeList[0]
103mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
104mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
105mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
106mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
107mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
108mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))
109
110#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111## add tip load
112tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeList[-1]))
113mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [1e6, 0, 0]))
114
115
116## assemble system and define simulation settings
117mbs.Assemble()
118
119simulationSettings = exu.SimulationSettings()
120
121tEnd = 1
122steps = 2000
123simulationSettings.timeIntegration.numberOfSteps = steps
124simulationSettings.timeIntegration.endTime = tEnd
125simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
126simulationSettings.timeIntegration.verboseMode = 1
127simulationSettings.solutionSettings.writeSolutionToFile = False
128
129#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
130simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
131
132simulationSettings.timeIntegration.newton.useModifiedNewton = True
133
134
135simulationSettings.staticSolver.newton.maxIterations = 50
136simulationSettings.staticSolver.numberOfLoadSteps = 10
137# simulationSettings.displayComputationTime = True
138# simulationSettings.displayStatistics = True
139simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
140
141## add some visualization settings
142SC.visualizationSettings.nodes.defaultSize = 0.005
143SC.visualizationSettings.bodies.beams.crossSectionFilled = False
144
145## start graphics
146if useGraphics:
147    exu.StartRenderer()
148    mbs.WaitForUserToContinue()
149
150## start static solver
151mbs.SolveStatic(simulationSettings)
152
153## print some output
154uLast = mbs.GetNodeOutput(nodeList[-1], exu.OutputVariableType.Coordinates)
155exu.Print("uTip =", uLast[0:2])
156
157## stop graphics
158if useGraphics:
159    SC.WaitForRenderEngineStopFlag()
160    exu.StopRenderer() #safely close rendering window!
161
162exu.Print('solution of LShapeGeomExactBeam2D=',uLast[1]) #use y-coordinate
163
164exudynTestGlobals.testError = uLast[1] - (-2.2115028353806547)
165exudynTestGlobals.testResult = uLast[1]