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 #only import if it does not conflict
 21import numpy as np
 22
 23useGraphics = True #without test
 24#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 25#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 26try: #only if called from test suite
 27    from modelUnitTests import exudynTestGlobals #for globally storing test results
 28    useGraphics = exudynTestGlobals.useGraphics
 29except:
 30    class ExudynTestGlobals:
 31        pass
 32    exudynTestGlobals = ExudynTestGlobals()
 33#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 34
 35#create an environment for mini example
 36SC = exu.SystemContainer()
 37mbs = SC.AddSystem()
 38
 39#some drawing properties, to see the objects ...
 40massStar = 1e20
 41massSatellite = 1e3
 42mass = [massSatellite, massSatellite, massSatellite, massSatellite]
 43sizeMass0 = 1e5 #just for drawing
 44sizeMass = [2e4,2e4,2e4,2e4] #just for drawing
 45rOrbit = [2e5, 4e5, 8e5, 10e5]
 46vOrbitEps = [1.,1.,1.,1.] #factors to make non-circular orbits...
 47
 48background = graphics.CheckerBoard(point=[0,0,-2*sizeMass0], size=2.1*max(rOrbit),
 49                                      color=[0,0,0,1], alternatingColor=[0.05,0,0,1])
 50
 51oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 52                                   visualization=VObjectGround(graphicsData=[background])))
 53# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 54node0 = mbs.AddNode(NodePoint(referenceCoordinates = [0,0,0])) #planet
 55gMass0 = graphics.Sphere(radius=1e5, color=graphics.color.blue, nTiles=64)
 56oMassPoint0 = mbs.AddObject(MassPoint(nodeNumber = node0, physicsMass=massStar,
 57                                      visualization=VMassPoint(graphicsData=[gMass0])))
 58m0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node0))
 59
 60#create satellites:
 61for i,r in enumerate(rOrbit):
 62    G = 6.6743e-11
 63    vOrbit = vOrbitEps[i]*np.sqrt(G*massStar/r) #orbital velocity, assumption of heavy star
 64    #exu.Print('vOrbit'+str(i)+'=',vOrbit)
 65    node1 = mbs.AddNode(NodePoint(referenceCoordinates = [r,0,0],
 66                                  initialVelocities=[0,vOrbit,0])) #satellite
 67
 68    gMass1 = graphics.Sphere(radius=sizeMass[i], color=graphics.colorList[i], nTiles=24)
 69
 70    oMassPoint1 = mbs.AddObject(MassPoint(nodeNumber = node1, physicsMass=mass[i],
 71                                          visualization=VMassPoint(graphicsData=[gMass1])))
 72
 73    m1 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node1))
 74
 75    mbs.AddObject(ObjectConnectorGravity(markerNumbers=[m0,m1],
 76                                         mass0 = massStar, mass1=mass[i]))
 77
 78#assemble and solve system for default parameters
 79mbs.Assemble()
 80simulationSettings = exu.SimulationSettings()
 81
 82tEnd = 1e6
 83h = 1000
 84simulationSettings.solutionSettings.writeSolutionToFile = False
 85simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 86simulationSettings.timeIntegration.endTime = tEnd
 87simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
 88simulationSettings.timeIntegration.newton.useModifiedNewton = True
 89
 90simulationSettings.displayStatistics = True
 91simulationSettings.timeIntegration.verboseMode = 1
 92
 93# SC.visualizationSettings.nodes.drawNodesAsPoint = False
 94
 95if useGraphics:
 96    exu.StartRenderer()              #start graphics visualization
 97    mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
 98
 99#start solver:
100# mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.TrapezoidalIndex2)
101#gives 7 digits of accuracy for tEnd=1e6, h=1e3:
102mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.RK67)
103
104if useGraphics:
105    SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
106    exu.StopRenderer()               #safely close rendering window!
107
108#check result at default integration time
109#node1 is last node
110pos = mbs.GetNodeOutput(node1, exu.OutputVariableType.Position)
111
112exudynTestGlobals.testResult = pos[0] + pos[1] + pos[2]
113
114exu.Print("result for ObjectConnectorGravity =", exudynTestGlobals.testResult)
115#exudynTestGlobals.testResult = 1014867.2330320379