pendulumGeomExactBeam2Dsimple.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  Example for GeometricallyExactBeam2D, connected with 2D revolute joint; uses GenerateStraightBeam
 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
20from exudyn.beams import GenerateBeamElementsAlongLine
21import exudyn.graphics as graphics
22
23import numpy as np
24# from math import sin, cos, pi
25
26## setup system container and mbs
27SC = exu.SystemContainer()
28mbs = SC.AddSystem()
29
30## define parameters for beams
31L = 0.5                # length of pendulum 
32E=1e8                  # very soft elastomer
33rho=1000               # elastomer
34h=0.002                # 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
39ks = 10*(1+nu)/(12+11*nu) # shear correction factor
40G = E/(2*(1+nu))          # shear modulus
41
42## create beam template with beam parameters
43beamTemplate = ObjectBeamGeometricallyExact2D(massPerLength=rho*A,
44                                              crossSectionInertia=rho*I,
45                                              bendingStiffness=E*I,
46                                              axialStiffness=E*A,
47                                              shearStiffness=ks*G*A,
48                                              visualization=VObjectBeamGeometricallyExact2D(drawHeight = h), )
49
50## create straight beam with 10 elements, apply gravity and fix (x,y) position of node 0 (rotation left free)
51beamInfo = GenerateBeamElementsAlongLine(mbs, positionStart=[0,0,0], positionEnd=[L,0,0], 
52                                numberOfElements=10, beamTemplate=beamTemplate,
53                                gravity=[0,-9.81,0], groundConstraintsStart=[1,1,0],)
54#beamInfo contains [nodeList, beamList, ...]
55
56## assemble system and define simulation settings
57mbs.Assemble()
58
59simulationSettings = exu.SimulationSettings()
60    
61tEnd = 1
62stepSize = 0.0025
63simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
64simulationSettings.timeIntegration.endTime = tEnd
65simulationSettings.timeIntegration.verboseMode = 1
66simulationSettings.solution.file.writePeriod = 0.005
67simulationSettings.solution.file.write = True
68
69simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
70simulationSettings.timeIntegration.newton.useModifiedNewton = True #for faster simulation
71
72
73## add some visualization settings
74SC.visualizationSettings.nodes.defaultSize = 0.01
75SC.visualizationSettings.nodes.drawNodesAsPoint = False
76SC.visualizationSettings.bodies.beams.crossSectionFilled = True
77
78## run dynamic simulation
79mbs.SolveDynamic(simulationSettings)
80
81## visualize computed solution:
82mbs.SolutionViewer()