beamTutorial.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Tutorial for GeometricallyExactBeam2D and ANCFCable2D
  5#
  6# Model:    Planar model of two highly flexible beams, modeled once with a geometrically exact beam and once with an ANCF cable element;
  7#           the beam has length 2m with h=0.005m, b=0.01m, E=1e9 and density rho=2000kg/m^3;
  8#           the shear deformable beam is rigidly attached to ground and the cable is rigidly attached to a moving ground.
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2024-02-13
 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
 24
 25## setup system container and mbs
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29## define parameters for beams
 30numberOfElements = 16
 31L = 2                  # length of pendulum
 32E=2e11                 # steel
 33rho=7800               # elastomer
 34h=0.005                # height of rectangular beam element in m
 35b=0.01                 # width of rectangular beam element in m
 36A=b*h                  # cross sectional area of beam element in m^2
 37I=b*h**3/12            # second moment of area of beam element in m^4
 38nu = 0.3               # Poisson's ratio
 39
 40EI = E*I
 41EA = E*A
 42rhoA = rho*A
 43rhoI = rho*I
 44ks = 10*(1+nu)/(12+11*nu) # shear correction factor
 45G = E/(2*(1+nu))          # shear modulus
 46GA = ks*G*A               # shear stiffness of beam
 47
 48g = [0,-9.81,0]           # gravity load
 49
 50positionOfNode0 = [0,0,0] # 3D vector
 51positionOfNode1 = [0+L,0,0] # 3D vector
 52
 53#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 54## build geometrically exact 2D beam template (Timoshenko-Reissner), which includes all parameters
 55beamTemplate = Beam2D(nodeNumbers = [-1,-1],
 56                      physicsMassPerLength=rhoA,
 57                      physicsCrossSectionInertia=rhoI,
 58                      physicsBendingStiffness=EI,
 59                      physicsAxialStiffness=EA,
 60                      physicsShearStiffness=GA,
 61                      physicsBendingDamping=0.02*EI,
 62                      visualization=VObjectBeamGeometricallyExact2D(drawHeight = h))
 63
 64beamData = GenerateStraightBeam(mbs, positionOfNode0, positionOfNode1,
 65                                numberOfElements, beamTemplate, gravity= g,
 66                                fixedConstraintsNode0=[1,1,1],
 67                                fixedConstraintsNode1=None)
 68
 69#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 70## build ANCF cable elemente (Bernoulli-Euler)
 71beamTemplate = Cable2D(nodeNumbers = [-1,-1],
 72                       physicsMassPerLength=rhoA,
 73                       physicsBendingStiffness=EI,
 74                       physicsAxialStiffness=EA,
 75                       physicsBendingDamping=0.02*EI,
 76                       visualization=VCable2D(drawHeight = h))
 77
 78cableData = GenerateStraightBeam(mbs, positionOfNode0, positionOfNode1,
 79                                 numberOfElements, beamTemplate, gravity= g,
 80                                 fixedConstraintsNode0=None,
 81                                 fixedConstraintsNode1=None)
 82
 83#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 84## create ground object to attach cable with generic joint
 85oGround = mbs.CreateGround(referencePosition=[0,0,0])
 86mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
 87
 88mCable = mbs.AddMarker(MarkerNodeRigid(nodeNumber=cableData[0][0]))
 89
 90## user function which represents translation and rotation in joint
 91def UFoffset(mbs, t, itemNumber, offsetUserFunctionParameters):
 92    x = SmoothStep(t, 2, 4, 0, 0.5)   #translate in local joint coordinates
 93    phi = SmoothStep(t, 5, 10, 0, pi) #rotates frame of mGround
 94    return [x, 0,0,0,0,phi]
 95
 96## add rigid joint (2D displacements and rotation around Z fixed)
 97mbs.AddObject(GenericJoint(markerNumbers=[mGround, mCable],
 98                           constrainedAxes=[1,1,0, 0,0,1],
 99                           offsetUserFunction=UFoffset,
100                           visualization=VGenericJoint(axesRadius=0.01,
101                                                       axesLength=0.02)))
102
103## assemble system and define simulation settings
104mbs.Assemble()
105
106simulationSettings = exu.SimulationSettings()
107
108tEnd = 10
109stepSize = 0.005
110simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
111simulationSettings.timeIntegration.endTime = tEnd
112simulationSettings.timeIntegration.verboseMode = 1
113simulationSettings.solutionSettings.solutionWritePeriod = 0.005
114simulationSettings.solutionSettings.writeSolutionToFile = True
115
116simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
117simulationSettings.timeIntegration.newton.useModifiedNewton = True #for faster simulation
118
119
120## add some visualization settings
121SC.visualizationSettings.nodes.defaultSize = 0.01
122SC.visualizationSettings.nodes.drawNodesAsPoint = False
123SC.visualizationSettings.bodies.beams.crossSectionFilled = True
124
125exu.StartRenderer()
126## run dynamic simulation
127mbs.SolveDynamic(simulationSettings)
128exu.StopRenderer()
129
130## visualize computed solution:
131mbs.SolutionViewer()