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