solutionViewerMultipleSimulations.py
You can view and download this file on Github: solutionViewerMultipleSimulations.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for multiple static solutions merged into one solution file
5#
6# Author: Johannes Gerstmayr
7# Date: 2022-12-20
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
13import exudyn as exu
14from exudyn.itemInterface import *
15from exudyn.utilities import *
16
17SC = exu.SystemContainer()
18mbs = SC.AddSystem()
19
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#set up simple ANCF model
22#background
23rect = [-0.5,-2,2.5,0.5] #xmin,ymin,xmax,ymax
24background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
25oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
26
27#cable:
28
29L=2 # length of ANCF element in m
30E=2.07e11 # Young's modulus of ANCF element in N/m^2
31rho=7800 # density of ANCF element in kg/m^3
32b=0.1 # width of rectangular ANCF element in m
33h=0.1 # height of rectangular ANCF element in m
34A=b*h # cross sectional area of ANCF element in m^2
35I=b*h**3/12 # second moment of area of ANCF element in m^4
36f=2*3*E*I/L**2 # tip load applied to ANCF element in N
37
38nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
39mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
40
41#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
42#generate ANCF beams with utilities function
43cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
44 physicsMassPerLength = rho*A,
45 physicsBendingStiffness = E*I,
46 physicsAxialStiffness = E*A,
47 useReducedOrderIntegration = 1,
48 #nodeNumbers = [0, 0], #will be filled in GenerateStraightLineANCFCable2D(...)
49 )
50
51positionOfNode0 = [0, 0, 0] # starting point of line
52positionOfNode1 = [L, 0, 0] # end point of line
53numberOfElements = 16
54
55#alternative to mbs.AddObject(Cable2D(...)) with nodes:
56ancf=GenerateStraightLineANCFCable2D(mbs,
57 positionOfNode0, positionOfNode1,
58 numberOfElements,
59 cableTemplate, #this defines the beam element properties
60 massProportionalLoad = [0,-9.81*0,0], #optionally add gravity
61 fixedConstraintsNode0 = [1,1,0,1], #add constraints for pos and rot (r'_y)
62 fixedConstraintsNode1 = [0,0,0,0])
63mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=ancf[0][-1])) #ancf[0][-1] = last node
64nLoad = mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [f*0, -f, 0])) #will be changed in load steps
65
66
67#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
68mbs.Assemble()
69# print(mbs)
70simulationSettings = exu.SimulationSettings() #takes currently set values or default values
71
72simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
73simulationSettings.solutionSettings.writeSolutionToFile = True
74simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/1000
75simulationSettings.displayComputationTime = False
76#simulationSettings.displayStatistics = True
77#simulationSettings.displayComputationTime = True
78
79SC.visualizationSettings.nodes.defaultSize = 0.01
80
81simulationSettings.solutionSettings.solutionInformation = "Cantilever"
82simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
83
84simulationSettings.staticSolver.verboseMode = 0
85simulationSettings.staticSolver.newton.newtonResidualMode = 1
86
87#adapt these settings for better solution file with multiple simulations:
88#**************************************************
89simulationSettings.solutionSettings.appendToFile = False
90simulationSettings.solutionSettings.writeFileFooter = False #never write footer as it would be seen between the solution steps
91#**************************************************
92
93useGraphics=False
94if useGraphics:
95 exu.StartRenderer()
96 mbs.WaitForUserToContinue()
97
98nLoadSteps = 25 #this is the number of individual computations; could also be done with staticSolver.numberOfLoadSteps
99 # but here, we want to show how to do multiple steps merged into one solution file
100for loadSteps in range(nLoadSteps):
101 #loadValue = f**((loadSteps+1)/nLoadSteps) #geometric increment of loads
102 loadValue = 2*f*(loadSteps+1)/(nLoadSteps)
103
104 mbs.SetLoadParameter(nLoad, 'loadVector', [0, -loadValue,0])
105 #print('load vector=' + str(mbs.GetLoadParameter(nLoad, 'loadVector')) )
106
107 simulationSettings.staticSolver.loadStepStart = loadSteps
108 # simulationSettings.staticSolver.numberOfLoadSteps = 5
109 mbs.SolveStatic(simulationSettings, updateInitialValues=True)
110
111 #**************************************************
112 #after first STEP, add this:
113 simulationSettings.solutionSettings.writeInitialValues = False #to avoid duplication of output times (start/end)
114 simulationSettings.solutionSettings.writeFileHeader = False
115 simulationSettings.solutionSettings.appendToFile = True
116 #**************************************************
117
118 sol = mbs.systemData.GetODE2Coordinates()
119
120 n = len(sol)
121 print('load=',loadValue, ', tip: x='+str(sol[n-4])+', y='+str(sol[n-3]))
122
123if useGraphics:
124 SC.WaitForRenderEngineStopFlag()
125 exu.StopRenderer() #safely close rendering window!
126
127if True:
128 #%%
129
130 t=LoadSolutionFile('solution/coordinatesSolution.txt', verbose=False, safeMode=True)
131 mbs.SolutionViewer(solution=t)