search

jgerstmayr / EXUDYN

Multibody Dynamics Simulation: Rigid and flexible multibody systems
https://www.youtube.com/channel/UCsQD2bIPBXB_4J23WtqKkVw/playlists
Other
173 stars 23 forks source link

How to make ANCFCable2D fixed to a rigid body? #1

Closed MingLocke closed 4 years ago

MingLocke commented 4 years ago

Problem Description

First of all, thank you for the development of EXUDYN program! I appreciate your work very much.

Recently, I tried to simulate a case, which can be described as a very flexible beam and its base rotating about the base's center.

But I don't know how to make a beam fixed to its rigid base, could you give me some advice?

rotating beam and base

Python Script

import exudyn as exu
from exudyn.itemInterface import *
from exudyn.utilities import *
import numpy as np

SC = exu.SystemContainer()
mbs = SC.AddSystem()

# cable
L = 5           # length of ANCF element in m
h = 2e-2        # height of rectangular ANCF element in m
b = 2e-2        # width of rectangular ANCF element in m
S = h*b         # cross sectional area of ANCF element in m^2
Iz = b*h**3/12  # second moment of area of ANCF element in m^4
rho = 2766.7    # density of ANCF element in kg/m^3
E = 68.952e9    # Young's modulus of ANCF element in N/m^2

# base
a = 1           # radius of the base
Joh = 5         # second moment of area of base in m^4

# torque
tau0 = 250      # N*m

# graphices for base
graphicsBase = GraphicsDataCylinder([0,0,-0.1], [0,0,0.2], a, color=[0.1,0.1,0.8,1])

iBase = InertiaCylinder(10/np.pi, 1, a, 2)

# base
[n0,b0]=AddRigidBody(mainSys = mbs, 
                     inertia = iBase,
                     nodeType = str(exu.NodeType.RotationEulerParameters), 
                     position = [0, 0, 0], 
                     rotationMatrix = np.diag([1,1,1]),
                     angularVelocity = [0,0,0],
                     graphicsDataList = [graphicsBase])

nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[a,0,1,0]))
nElements = 10
lElem = L / nElements
for i in range(nElements):
    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1)+a,0,1,0]))
    mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*S, 
                          physicsBendingStiffness=E*Iz, 
                          physicsAxialStiffness=E*S, 
                          nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))

# ground body and marker
oGround = mbs.AddObject(ObjectGround())
markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))

# markers for base:
markerBase = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[0,0,0]))
marker0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[a,0,0]))
# bmarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=0))
# bmarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=1))
# bmarker2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=5))

# markers for beam:
marker1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nc0))
# mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
# mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
# mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))

# revolute joint (free z-axis)
mbs.AddObject(GenericJoint(markerNumbers=[markerGround, markerBase], 
                           constrainedAxes=[1,1,1,1,1,0],
                           visualization=VObjectJointGeneric(axesRadius=0.01, axesLength=0.1)))

# # fixed constraint
mbs.AddObject(GenericJoint(markerNumbers=[marker0, marker1],
                            constrainedAxes=[1,1,1,1,1,1],
                            visualization=VObjectJointGeneric(axesRadius=0.01, axesLength=0.1)))
# mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker0,mANCF0]))
# mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker1,mANCF1]))
# mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker2,mANCF2]))

# add torque
def userTorque(t, loadVector):
    if t <= 2:
        return [0,0,loadVector[2]*np.sin(2*np.pi*t/2)]
    else:
        return [0,0,0]

mbs.AddLoad(Torque(markerNumber=markerBase,
                   loadVector=[0,0,tau0],
                   loadVectorUserFunction=userTorque))

mbs.Assemble()

simulationSettings = exu.SimulationSettings() #takes currently set values or default values

simulationSettings.timeIntegration.numberOfSteps = 100000
simulationSettings.timeIntegration.endTime = 4 #0.2 for testing

simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True

exu.StartRenderer()
mbs.WaitForUserToContinue()

SC.TimeIntegrationSolve(mbs, 'GeneralizedAlpha', simulationSettings)

SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
jgerstmayr commented 4 years ago 
Dear Luo,
nice example.
The reason, why the constraint does not work is that the beam is
  2D. For this reason, there may be only three constraints active in
  the generic joint:
[1,1,0, 0,0,1]
which is the x/y position and the rotation around the z-axis.
I modified this in the attached example. Works nicely.
REMARK: instead of the 3D body, you could also use a 2D body,
  which would fit better to your example.

Kind regards,
Johannes

Am 27.09.2020 um 15:54 schrieb Luo
  Ming:

    import exudyn as exu

from exudyn.itemInterface import from exudyn.utilities import import numpy as np

SC = exu.SystemContainer() mbs = SC.AddSystem()

cable

L = 5 # length of ANCF element in m h = 2e-2 # height of rectangular ANCF element in m b = 2e-2 # width of rectangular ANCF element in m S = hb # cross sectional area of ANCF element in m^2 Iz = bh**3/12 # second moment of area of ANCF element in m^4 rho = 2766.7 # density of ANCF element in kg/m^3 E = 68.952e9 # Young's modulus of ANCF element in N/m^2

base

a = 1 # radius of the base Joh = 5 # second moment of area of base in m^4

torque

tau0 = 250 # N*m

graphices for base

graphicsBase = GraphicsDataCylinder([0,0,-0.1], [0,0,0.2], a, color=[0.1,0.1,0.8,1])

iBase = InertiaCylinder(10/np.pi, 1, a, 2)

base

[n0,b0]=AddRigidBody(mainSys = mbs, inertia = iBase, nodeType = str(exu.NodeType.RotationEulerParameters), position = [0, 0, 0], rotationMatrix = np.diag([1,1,1]), angularVelocity = [0,0,0], graphicsDataList = [graphicsBase])

nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[a,0,1,0])) nElements = 10 lElem = L / nElements for i in range(nElements): nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem(i+1)+a,0,1,0])) mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rhoS, physicsBendingStiffness=EIz, physicsAxialStiffness=ES, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))

ground body and marker

oGround = mbs.AddObject(ObjectGround()) markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))

markers for base:

markerBase = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[0,0,0])) marker0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[a,0,0]))

bmarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=0))

bmarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=1))

bmarker2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=5))

markers for beam:

marker1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nc0))

mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))

mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))

mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))

revolute joint (free z-axis)

mbs.AddObject(GenericJoint(markerNumbers=[markerGround, markerBase], constrainedAxes=[1,1,1,1,1,0], visualization=VObjectJointGeneric(axesRadius=0.01, axesLength=0.1)))

fixed constraint

mbs.AddObject(GenericJoint(markerNumbers=[marker0, marker1], constrainedAxes=[1,1,1,1,1,1], visualization=VObjectJointGeneric(axesRadius=0.01, axesLength=0.1)))

mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker0,mANCF0]))

mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker1,mANCF1]))

mbs.AddObject(CoordinateConstraint(markerNumbers=[bmarker2,mANCF2]))

add torque

def userTorque(t, loadVector): if t <= 2: return [0,0,loadVector[2]np.sin(2np.pi*t/2)] else: return [0,0,0]

mbs.AddLoad(Torque(markerNumber=markerBase, loadVector=[0,0,tau0], loadVectorUserFunction=userTorque))

mbs.Assemble()

simulationSettings = exu.SimulationSettings() #takes currently set values or default values

simulationSettings.timeIntegration.numberOfSteps = 100000 simulationSettings.timeIntegration.endTime = 4 #0.2 for testing

simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True

exu.StartRenderer() mbs.WaitForUserToContinue()

SC.TimeIntegrationSolve(mbs, 'GeneralizedAlpha', simulationSettings)

SC.WaitForRenderEngineStopFlag() exu.StopRenderer() #safely close rendering window!

--

Univ.-Prof. DI Dr. Johannes Gerstmayr Department of Mechatronics, Group of Machine Elements and Design, Leopold-Franzens-Universität Innsbruck, Technikerstraße 13 (3. Stock, Raum 322), 6020 Innsbruck, Österreich Tel.: +43 (0) 512 507 62750 Tel. (mobil): +43 (0) 676 8725 62750 E-Mail: johannes.gerstmayr@uibk.ac.at

EXUDYN - Flexible Multibody Dynamics Simulation with Python and C++ https://github.com/jgerstmayr/EXUDYN

MingLocke commented 4 years ago

It worked! Thank you. I will take your advice to make a 2D body later.

oridong commented 3 years ago

I try to run this , but I get this

Exception has occurred: AttributeError (note: full exception trace is shown but execution is paused at: _run_module_as_main) 'exudyn.exudynCPP.SystemContainer' object has no attribute 'TimeIntegrationSolve' File "D:\vsc_dyn\EXUDYN\myexp\cable_rotbase.py", line 105, in SC.TimeIntegrationSolve(mbs, 'GeneralizedAlpha', simulationSettings)

version 1.0.276