ANCFcable2DuserFunction.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for ANCFCable2D test user functions
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2023-12-13
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13
 14import exudyn as exu
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics
 17
 18import numpy as np
 19
 20testIsActive = exu.sys.get('testIsActive', False)
 21
 22from exudyn.beams import *
 23from math import atan
 24
 25#create an environment for mini example
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 30nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 31
 32rhoA = 78.
 33EA = 100000.
 34EI = 2000
 35
 36#example of bending moment user function
 37#limit bending moment with atan function
 38def bendingMomentUserFunction(mbs, t, itemNumber, axialPositionNormalized, curvature, curvature_t, curvatureRef, bendingStiffness, bendingDamping,
 39                                    axialStrain, axialStrain_t, axialStrainRef):
 40    #m = bendingStiffness*(curvature-curvatureRef) + bendingDamping*curvature_t #this is the linear, conventional case
 41    kappa=(curvature-curvatureRef)
 42    kappa = 0.1*atan(10*kappa) #nonlinear behavior, somehow like elasto-plastic
 43    return bendingStiffness*(kappa) + bendingDamping*curvature_t
 44
 45#example of axial force user function
 46#reduce stiffness over time
 47def axialForceUserFunction(mbs, t, itemNumber, axialPositionNormalized, axialStrain, axialStrain_t, axialStrainRef, axialStiffness, axialDamping,
 48            curvature, curvature_t, curvatureRef):
 49    fact = max(0.02,(2-t**0.5)) #make axial stiffness it softer over time
 50    return fact*axialStiffness*(axialStrain-axialStrainRef) + axialDamping*axialStrain_t
 51
 52#create ANCF cable object:
 53cable = ObjectANCFCable2D(massPerLength=rhoA, 
 54                bendingStiffness=EI, 
 55                bendingDamping = EI*0.1,
 56                axialStiffness=EA,
 57                axialDamping=EA*0.05,
 58                bendingMomentUserFunction=bendingMomentUserFunction,
 59                axialForceUserFunction=axialForceUserFunction,
 60                )
 61
 62#create several cable elements
 63ancf=GenerateBeamElementsAlongLine(mbs=mbs,
 64                positionStart=[0,0,0], positionEnd=[2,0,0],
 65                numberOfElements=16, #converged to 4 digits
 66                beamTemplate=cable, #this defines the beam element properties
 67                gravity = [0,-9.81,0],
 68                groundConstraintsStart = [1,1, 0,1],
 69                )
 70
 71#assemble and solve system for default parameters
 72mbs.Assemble()
 73
 74endTime = 0.5
 75stepSize = 5e-3
 76
 77simulationSettings = exu.SimulationSettings()
 78
 79simulationSettings.solution.file.write = False
 80simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
 81simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
 82# simulationSettings.show.computationTime = True
 83# simulationSettings.show.statistics = True
 84
 85simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
 86simulationSettings.timeIntegration.endTime = endTime
 87simulationSettings.timeIntegration.newton.useModifiedNewton = True
 88
 89SC.visualizationSettings.view0.window.renderWindowSize=[1200,1024]
 90#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
 91
 92
 93if not testIsActive:
 94    SC.renderer.Start()              #start graphics visualization
 95    SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
 96
 97mbs.SolveDynamic(simulationSettings)
 98
 99if not testIsActive:
100    SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
101    SC.renderer.Stop()               #safely close rendering window!
102
103#evaluate final (=current) output values
104node = ancf['nodes'][-1]
105p = mbs.GetNodeOutput(node, exu.OutputVariableType.Position)
106exu.Print('ANCFcable2DuserFunction test tip pos=',p)
107
108u=sum(p)
109exu.Print('solution of ANCFcable2DuserFunction test =',u)
110
111exu.sys['testResult'] = u