ANCFCable2DBendingTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A cantilever of `ANCFCable2D` elements, clamped at one end and pulled down by a tip
  5#           force. Solved dynamically once and statically twice, the second time with the Newton
  6#           tolerances at 1e-14 so that it converges onto the MATLAB reference it is compared with.
  7#           The model compares against a reference value written into it, so its
  8#           result is that difference and its reference solution is 0 (#2632).
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2019-11-01, reworked 2026-09-24
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics
 20
 21testIsActive = exu.sys.get('testIsActive', False)
 22exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28totalError = 0 #total error computed in test
 29
 30L=2                    # length of ANCF element in m
 31Em=2.07e11             # Young's modulus of ANCF element in N/m^2
 32rho=7800               # density of ANCF element in kg/m^3
 33b=0.1                  # width of rectangular ANCF element in m
 34h=0.1                  # height of rectangular ANCF element in m
 35A=b*h                  # cross sectional area of ANCF element in m^2
 36I=b*h**3/12            # second moment of area of ANCF element in m^4
 37f=3*Em*I/L**2           # tip load applied to ANCF element in N
 38
 39EI = Em*I
 40rhoA = rho*A
 41EA = Em*A
 42
 43
 44nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 45
 46nElements = 1
 47lElem = L / nElements
 48
 49for i in range(nElements):
 50    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 51    elem = mbs.AddObject(Cable2D(length=lElem, massPerLength=rhoA, 
 52                                 bendingStiffness=EI, axialStiffness=EA, 
 53                                 nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 54
 55#tip node / force
 56mANCFnode = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast))
 57mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, -f, 0]))
 58
 59
 60#ground node / coordinate
 61nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 62mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action; coordinate number does not matter
 63
 64#constraints
 65mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 66mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 67mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3)) #3
 68mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 69mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 70mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 71#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0],velocityLevel=True))
 72#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1],velocityLevel=True))
 73#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2],velocityLevel=True))
 74
 75#Assemble and Solve:
 76mbs.Assemble()
 77
 78simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 79simulationSettings.solution.file.name = "solution/ANCFCable2D_bending_test.txt"
 80simulationSettings.solution.file.write=False
 81simulationSettings.timeIntegration.numberOfSteps = 1000
 82#simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/1000
 83simulationSettings.timeIntegration.endTime = 0.1
 84simulationSettings.timeIntegration.verboseMode = 0
 85simulationSettings.timeIntegration.newton.useModifiedNewton = False
 86simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6
 87simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
 88simulationSettings.show.statistics = True
 89
 90if not testIsActive: 
 91    SC.renderer.Start()
 92
 93mbs.SolveDynamic(simulationSettings)
 94
 95sol = mbs.systemData.GetODE2Coordinates(); n = len(sol)
 96#tip displacements:
 97u = sol[n-4];v = sol[n-3] #15.12.2019:(-1.040678127615946 -1.444419986874761); 20.10.2019: (-1.0406781266430292 -1.4444199866881322); 17.10.2019: sol= -1.040678126647053 -1.4444199866858678; #28.7.2009: -1.040678126643273, -1.444419986688082
 98if True:
 99    totalError += u+v - (-1.0611305415798655 -1.4396907464589017)  #2021-09-27: new JacobianODE2RHS
100    #totalError += u+v -(-1.0611305415779122 -1.4396907464583975)  #2021-02-04: (-1.0611305415779122 -1.4396907464583975)
101else:
102    totalError += u+v - (-1.0391620828192676 -1.4443521331339881)  #2019-12-26: (-1.0391620828192676 -1.4443521331339881) old (before correct initial accelerations): (-1.040678127615946 -1.444419986874761)
103exu.Print('sol dynamic=',u,v)
104#exu.Print('time integration error =',totalError)
105
106mbs.SolveStatic(simulationSettings)
107
108sol = mbs.systemData.GetODE2Coordinates(); n = len(sol)
109u = sol[n-4]; v = sol[n-3]; #20.10.2019: -0.3622447299987188 -0.9941447593196007; 17.10.2019: sol= -0.3622447299990847 -0.9941447593206921; #28.7.2019: -0.3622447300008477, -0.994144759326213
110exu.Print('sol static (standardTol)=',u,v)
111totalError += u+v - (-0.36224473018839665 -0.9941447595447153) #2021-09-27: new JacobianODE2RHS
112# totalError += u+v-(-0.3622447299987188 -0.9941447593196007)
113
114simulationSettings.staticSolver.newton.relativeTolerance = 1e-14 #in order to converge to MATLAB results
115simulationSettings.staticSolver.newton.absoluteTolerance = 1e-14
116
117mbs.SolveStatic(simulationSettings)
118
119sol = mbs.systemData.GetODE2Coordinates(); n = len(sol)
120#tip displacements: paper GerstmIschrik2008: 1Element: u=-0.362244729891,  v=-0.994144758725; 4 Elements: 0.507428715119 1.205533702233
121u = sol[n-4]; v = sol[n-3];                 #2019-12-17: -0.3622447298904951 -0.9941447587249616
122exu.Print('sol static (tol=1e-14)=',u,v)
123totalError += u+v - (-0.36224472989050654 -0.9941447587249747) #2021-09-27: new JacobianODE2RHS
124# totalError += u+v-(-0.3622447298904951 -0.9941447587249616)
125
126
127#totalError -= -1.3563894893270607-2.4850981133313548 #reference solution with one element and standard settings except: gen-alpha=0.6, useModifiedNewton = False
128
129if not testIsActive: 
130    SC.renderer.DoIdleTasks()
131    SC.renderer.Stop() #safely close rendering window!
132
133testResult = totalError
134exu.Print('solution of ANCFCable2DBendingTest=', testResult)
135exu.sys['testResult'] = testResult