connectorGravityTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for ObjectConnectorGravity, realizing gravitational forces between 
  5#           two masses (used for astrospace or small-scale astronomical investigations, e.g., satellites);
  6#           in this case we simulate a small planetary system, initializing planets with orbital velocity
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2022-01-30
 10#
 11# 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.
 12#
 13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 14import sys
 15sys.path.append('../TestModels')
 16
 17import exudyn as exu
 18from exudyn.itemInterface import *
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics
 21import numpy as np
 22
 23testIsActive = exu.sys.get('testIsActive', False)
 24
 25#create an environment for mini example
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29#some drawing properties, to see the objects ...
 30massStar = 1e20
 31massSatellite = 1e3
 32mass = [massSatellite, massSatellite, massSatellite, massSatellite]
 33sizeMass0 = 1e5 #just for drawing
 34sizeMass = [2e4,2e4,2e4,2e4] #just for drawing
 35rOrbit = [2e5, 4e5, 8e5, 10e5]
 36vOrbitEps = [1.,1.,1.,1.] #factors to make non-circular orbits...
 37
 38background = graphics.CheckerBoard(point=[0,0,-2*sizeMass0], size=2.1*max(rOrbit), 
 39                                      color=[0,0,0,1], alternatingColor=[0.05,0,0,1])
 40
 41oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 42                                   visualization=VObjectGround(graphicsData=[background])))
 43# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 44node0 = mbs.AddNode(NodePoint(referenceCoordinates = [0,0,0])) #planet
 45gMass0 = graphics.Sphere(radius=1e5, color=graphics.color.blue, nTiles=32)
 46oMassPoint0 = mbs.AddObject(MassPoint(nodeNumber = node0, mass=massStar,
 47                                      visualization=VMassPoint(graphicsData=[gMass0])))
 48m0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node0))
 49
 50#create satellites:
 51for i,r in enumerate(rOrbit):
 52    G = 6.6743e-11
 53    vOrbit = vOrbitEps[i]*np.sqrt(G*massStar/r) #orbital velocity, assumption of heavy star
 54    #exu.Print('vOrbit'+str(i)+'=',vOrbit)
 55    node1 = mbs.AddNode(NodePoint(referenceCoordinates = [r,0,0], 
 56                                  initialVelocities=[0,vOrbit,0])) #satellite
 57    
 58    gMass1 = graphics.Sphere(radius=sizeMass[i], color=graphics.colorList[i], nTiles=24)
 59    
 60    oMassPoint1 = mbs.AddObject(MassPoint(nodeNumber = node1, mass=mass[i],
 61                                          visualization=VMassPoint(graphicsData=[gMass1])))
 62    
 63    m1 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node1))
 64    
 65    mbs.AddObject(ObjectConnectorGravity(markerNumbers=[m0,m1],
 66                                         mass0 = massStar, mass1=mass[i]))
 67
 68#assemble and solve system for default parameters
 69mbs.Assemble()
 70simulationSettings = exu.SimulationSettings()
 71
 72tEnd = 1e6
 73h = 1000
 74simulationSettings.solution.file.write = False
 75simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 76simulationSettings.timeIntegration.endTime = tEnd
 77simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
 78simulationSettings.timeIntegration.newton.useModifiedNewton = True
 79
 80simulationSettings.show.statistics = True
 81simulationSettings.timeIntegration.verboseMode = 1
 82
 83# SC.visualizationSettings.nodes.drawNodesAsPoint = False
 84
 85if not testIsActive:
 86    SC.renderer.Start()              #start graphics visualization
 87    SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
 88
 89#start solver:
 90# mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.TrapezoidalIndex2)
 91#gives 7 digits of accuracy for tEnd=1e6, h=1e3:
 92mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.RK67)
 93
 94if not testIsActive:
 95    SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
 96    SC.renderer.Stop()               #safely close rendering window!
 97
 98#check result at default integration time
 99#node1 is last node
100pos = mbs.GetNodeOutput(node1, exu.OutputVariableType.Position)
101
102exu.sys['testResult'] = pos[0] + pos[1] + pos[2]
103
104exu.Print("result for ObjectConnectorGravity =", exu.sys['testResult'])
105#exu.sys['testResult'] = 1014867.2330320379