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