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 #only import if it does not conflict
 17
 18import numpy as np
 19
 20useGraphics = True #without test
 21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 23try: #only if called from test suite
 24    from modelUnitTests import exudynTestGlobals #for globally storing test results
 25    useGraphics = exudynTestGlobals.useGraphics
 26except:
 27    class ExudynTestGlobals:
 28        pass
 29    exudynTestGlobals = ExudynTestGlobals()
 30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 31
 32from exudyn.beams import *
 33from math import atan
 34
 35#create an environment for mini example
 36SC = exu.SystemContainer()
 37mbs = SC.AddSystem()
 38
 39oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 40nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 41
 42rhoA = 78.
 43EA = 100000.
 44EI = 2000
 45
 46#example of bending moment user function
 47#limit bending moment with atan function
 48def bendingMomentUserFunction(mbs, t, itemNumber, axialPositionNormalized, curvature, curvature_t, curvatureRef, physicsBendingStiffness, physicsBendingDamping,
 49                                    axialStrain, axialStrain_t, axialStrainRef):
 50    #m = physicsBendingStiffness*(curvature-curvatureRef) + physicsBendingDamping*curvature_t #this is the linear, conventional case
 51    kappa=(curvature-curvatureRef)
 52    kappa = 0.1*atan(10*kappa) #nonlinear behavior, somehow like elasto-plastic
 53    return physicsBendingStiffness*(kappa) + physicsBendingDamping*curvature_t
 54
 55#example of axial force user function
 56#reduce stiffness over time
 57def axialForceUserFunction(mbs, t, itemNumber, axialPositionNormalized, axialStrain, axialStrain_t, axialStrainRef, physicsAxialStiffness, physicsAxialDamping,
 58            curvature, curvature_t, curvatureRef):
 59    fact = max(0.02,(2-t**0.5)) #make axial stiffness it softer over time
 60    return fact*physicsAxialStiffness*(axialStrain-axialStrainRef) + physicsAxialDamping*axialStrain_t
 61
 62#create ANCF cable object:
 63cable = ObjectANCFCable2D(physicsMassPerLength=rhoA,
 64                physicsBendingStiffness=EI,
 65                physicsBendingDamping = EI*0.1,
 66                physicsAxialStiffness=EA,
 67                physicsAxialDamping=EA*0.05,
 68                bendingMomentUserFunction=bendingMomentUserFunction,
 69                axialForceUserFunction=axialForceUserFunction,
 70                )
 71
 72#create several cable elements
 73ancf=GenerateStraightLineANCFCable(mbs=mbs,
 74                positionOfNode0=[0,0,0], positionOfNode1=[2,0,0],
 75                numberOfElements=16, #converged to 4 digits
 76                cableTemplate=cable, #this defines the beam element properties
 77                massProportionalLoad = [0,-9.81,0],
 78                fixedConstraintsNode0 = [1,1, 0,1],
 79                )
 80
 81#assemble and solve system for default parameters
 82mbs.Assemble()
 83
 84endTime = 0.5
 85stepSize = 5e-3
 86
 87simulationSettings = exu.SimulationSettings()
 88
 89simulationSettings.solutionSettings.writeSolutionToFile = False
 90simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
 91simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
 92# simulationSettings.displayComputationTime = True
 93# simulationSettings.displayStatistics = True
 94
 95simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
 96simulationSettings.timeIntegration.endTime = endTime
 97simulationSettings.timeIntegration.newton.useModifiedNewton = True
 98
 99SC.visualizationSettings.window.renderWindowSize=[1200,1024]
100#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
101
102
103if useGraphics:
104    exu.StartRenderer()              #start graphics visualization
105    mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
106
107mbs.SolveDynamic(simulationSettings)
108
109if useGraphics:
110    SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
111    exu.StopRenderer()               #safely close rendering window!
112
113#evaluate final (=current) output values
114node = ancf[0][-1]
115p = mbs.GetNodeOutput(node, exu.OutputVariableType.Position)
116exu.Print('ANCFcable2DuserFunction test tip pos=',p)
117
118u=sum(p)
119exu.Print('solution of ANCFcable2DuserFunction test =',u)
120
121exudynTestGlobals.testError = u - (0.6015588367721232)
122exudynTestGlobals.testResult = u