pendulumGeomExactBeam2D.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example for GeometricallyExactBeam2D, connected with 2D revolute joint; pendulum is modeled with 10 element
  5#
  6# Model:    Planar model of a highly flexible pendulum of length 0.5m with h=0.002m, b=0.01m, E=1e8 and density rho=1000kg/m^3;
  7#           The pendulum is released from the horizontal position under gravity acting 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 libaries
 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
 23# from math import sin, cos, pi
 24
 25## setup system container and mbs
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29## define parameters for beams
 30nElements = 10
 31L = 0.5                # length of pendulum
 32lElem = L/nElements    # length of one finite element
 33E=1e8                  # very soft elastomer
 34rho=1000               # elastomer
 35h=0.002                # height of rectangular beam element in m
 36b=0.01                 # width 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
 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
 49g = [0,-9.81,0]           # gravity load
 50
 51## create nodes (one more than number of elements)
 52for i in range(nElements+1):
 53    pRef = [i*lElem,0,0]
 54    n = mbs.AddNode(NodeRigidBody2D(referenceCoordinates = pRef))
 55    if i==0: firstNode = n
 56
 57## create beam elements:
 58listBeams = []
 59for i in range(nElements):
 60    oBeam = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [firstNode+i,firstNode+i+1],
 61                                                            physicsLength=lElem,
 62                                                            physicsMassPerLength=rhoA,
 63                                                            physicsCrossSectionInertia=rhoI,
 64                                                            physicsBendingStiffness=EI,
 65                                                            physicsAxialStiffness=EA,
 66                                                            physicsShearStiffness=GA,
 67                                                            visualization=VObjectBeamGeometricallyExact2D(drawHeight = h)
 68                                                ))
 69    listBeams += [oBeam]
 70
 71## create ground node with marker for coordinate constraints
 72oGround = mbs.CreateGround(referencePosition=[0,0,0])
 73mGround = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0]))
 74
 75## add revolute joint between ground and first node
 76mNode0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=firstNode))
 77mbs.AddObject(ObjectJointRevolute2D(markerNumbers=[mGround, mNode0]))
 78
 79#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 80## add gravity loading for beams
 81for beam in listBeams:
 82    marker = mbs.AddMarker(MarkerBodyMass(bodyNumber=beam))
 83    mbs.AddLoad(LoadMassProportional(markerNumber=marker, loadVector=g))
 84
 85
 86## assemble system and define simulation settings
 87mbs.Assemble()
 88
 89simulationSettings = exu.SimulationSettings()
 90
 91tEnd = 5
 92stepSize = 0.0025
 93simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
 94simulationSettings.timeIntegration.endTime = tEnd
 95simulationSettings.timeIntegration.verboseMode = 1
 96simulationSettings.solutionSettings.solutionWritePeriod = 0.005
 97simulationSettings.solutionSettings.writeSolutionToFile = True
 98
 99simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
100simulationSettings.timeIntegration.newton.useModifiedNewton = True #for faster simulation
101
102
103## add some visualization settings
104SC.visualizationSettings.nodes.defaultSize = 0.01
105SC.visualizationSettings.nodes.drawNodesAsPoint = False
106SC.visualizationSettings.bodies.beams.crossSectionFilled = True
107
108## run dynamic simulation
109mbs.SolveDynamic(simulationSettings)
110
111## visualize computed solution:
112mbs.SolutionViewer()