LShapeGeomExactBeam2D.py

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Test model for GeometricallyExactBeam2D, evaluating behavior of rotated beam
#
# Model:    Planar model of L-shaped beams with tip load; the beam has 4 elements at each part of the L-shape:
#           the element length is 0.25m, the material is steel and height is 0.05m and with 0.1m; 
#           the L-shape is fixed at the left end to ground and a tip load of [1e6,0,0] is applied to the tip.
#
# Author:   Johannes Gerstmayr
# Date:     2021-03-25
#
# 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.
#
# *clean example*
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

## import libaries
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics

import numpy as np
from math import sin, cos, pi

testIsActive = exu.sys.get('testIsActive', False)

## setup system container and mbs
SC = exu.SystemContainer()
mbs = SC.AddSystem()

## define parameters for beams
a = 0.25
lElem = a              # length of one finite element 
E=2.1e11               # Steel; Young's modulus of beam element in N/m^2
rho=7800               # Steel; density of beam element in kg/m^3
b=0.1                  # width of rectangular beam element in m
h=0.05                 # height of rectangular beam element in m
A=b*h                  # cross sectional area of beam element in m^2
I=b*h**3/12            # second moment of area of beam element in m^4
nu = 0.3               # Poisson's ratio for steel
    
EI = E*I
EA = E*A
rhoA = rho*A
rhoI = rho*I
ks = 10*(1+nu)/(12+11*nu) # shear correction factor
G = E/(2*(1+nu))          # shear modulus
GA = ks*G*A               # shear stiffness of beam


## create position list for nodes:
nodeList=[]
pRefList=[[0,0],
          [1*a,0],
          [2*a,0],
          [3*a,0],
          [4*a,0],
          [4*a,1*a],
          [4*a,2*a],
          [4*a,3*a],
          [4*a,4*a],
          ]

## create nodes
for p in pRefList:
    pRef = [p[0],p[1],0] #always angle zero, which is not considered in beam for includeReferenceRotations=False
    ni=mbs.AddNode(NodeRigidBody2D(referenceCoordinates = pRef, 
                                   initialCoordinates = [0,0,0], 
                                   initialVelocities= [0,0,0]))
    nodeList += [ni]

## create beam elements:
for i in range(len(nodeList)-1):
    mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [nodeList[i],nodeList[i+1]], 
                                                 length=lElem,
                                                 massPerLength=rhoA,
                                                 crossSectionInertia=rhoI,
                                                 bendingStiffness=EI,
                                                 axialStiffness=EA,
                                                 shearStiffness=GA,
                                                 includeReferenceRotations=False, #to connect beams at 90° at same node
                                                 visualization=VObjectBeamGeometricallyExact2D(drawHeight = h) ))
    

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## create ground node with marker for coordinate constraints
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))    
mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))

## add markers and constraints for fixed support
n0 = nodeList[0]
mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## add tip load
tipNodeMarker = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodeList[-1]))
mbs.AddLoad(Force(markerNumber = tipNodeMarker, loadVector = [1e6, 0, 0]))
    

## assemble system and define simulation settings
mbs.Assemble()

simulationSettings = exu.SimulationSettings()
    
tEnd = 1
steps = 2000
simulationSettings.timeIntegration.numberOfSteps = steps
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solution.file.writePeriod = tEnd/steps
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.solution.file.write = False

#simulationSettings.solution.file.writePeriod = tEnd/steps
simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
    
simulationSettings.timeIntegration.newton.useModifiedNewton = True


simulationSettings.staticSolver.newton.maxIterations = 50
simulationSettings.staticSolver.numberOfLoadSteps = 10
# simulationSettings.show.computationTime = True
# simulationSettings.show.statistics = True
simulationSettings.staticSolver.newton.relativeTolerance = 1e-6

## add some visualization settings
SC.visualizationSettings.nodes.defaultSize = 0.005
SC.visualizationSettings.bodies.beams.crossSectionFilled = False

## start graphics
if not testIsActive:
    SC.renderer.Start()
    SC.renderer.DoIdleTasks()

## start static solver
mbs.SolveStatic(simulationSettings)

## print some output        
uLast = mbs.GetNodeOutput(nodeList[-1], exu.OutputVariableType.Coordinates)
exu.Print("uTip =", uLast[0:2])

## stop graphics
if not testIsActive:
    SC.renderer.DoIdleTasks()
    SC.renderer.Stop() #safely close rendering window!

exu.Print('solution of LShapeGeomExactBeam2D=',uLast[1]) #use y-coordinate

exu.sys['testResult'] = uLast[1]