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