Hi,
I am currently working on a program to compute the evolution of planetary rings over large time scale. For the moment I’m trying to replicate the results of an algorithm described p6 of this article :
https://hal.science/hal-00676215/document
While the article is using a Runge-Kutta method I choose to use an explicit Euler method.
The problems I encounter are :
- the evolution is too slow
- As the time increases, my Gaussian function shifts on the left while its center should not move.
Here is the code :
import numpy as np
import matplotlib.pyplot as plt
G = 6.67430e-11 #gravitational constant
Ms = 5.6834e26 #Saturn's mass in kg
Mm = 3.84e19 #Mima's mass en kg
rhoP = 900.0 #particle density
rp = 1.0 #particle radius in m
mp = rhoP*4/3*np.pi*rp**3 #particle mass in kg
nu = 0.1 #constant viscosity in m2.s-1
dispR = np.sqrt(G*mp/rp) #radial dispertion en m.s-1
Rmin, Rmax = 60e6, 140e6
dr = 100e3 #element size in m
Rmid = np.arange(Rmin, Rmax+dr, dr) #position of elements mid points
RInf = np.arange(Rmin-dr/2, Rmax, dr) #position of elements lower limit
RSup = np.arange(Rmin+dr/2, Rmax+1+dr/2, dr) #position of elements upper limit
Nr = int((Rmax-Rmin)/dr) #number of elements for the finite element method
Tmin, Tmax = 0, 1e5*365.25*24*3600
dt = 0.5*(dr**2) #time step
T = np.arange(Tmin, Tmax, dt)
Nt = int((Tmax-Tmin)/dt) #number of time steps
#omega = np.sqrt(G*Ms/(R**3)) #orbital frequency in Hz
#Initial state
def repartSigma(centre, width):
return 6*np.exp(-(Rdemi-centre)**2 / (2*(largeur**2)))
#evolution
def f(nu, sigma) :
return nu*sigma*np.sqrt(Rdemi)
def calcDerivF(nu, sigma) :
fct = f(nu, sigma)
derivF = np.zeros_like(Rdemi)
for i in range (Nr) :
derivF[i] = (fct[i+1]-fct[i])/(dr)
return derivF
def calcFlux(R, nu, sigma) :
derivF = calcDerivF(nu, sigma)
flux = 6*np.pi*np.sqrt(R)*derivF
return flux
def calcDerivSigma(nu, sigma):
fluxInf = calcFlux(RInf, nu, sigma)
fluxSup = calcFlux(RSup, nu, sigma)
fluxInf[0] = fluxSup[0]
fluxSup[-1] = fluxInf[-1]
deriveeSigma = (fluxSup - fluxInf)/(2*np.pi*Rmid*dr)
return derivSigma
def euler(nu, centre, witdh) :
sigma = repartSigma(centre, width)
for i in range (Nt) :
derivSigma = calcDerivSigma(nu, sigma)
sigma += derivSigma*dt
print (i/Nt) #progress check
return sigma
def calcMass(R, sigma) :
mass = 0
for i in range (Nr) :
mass += sigma[i]*2*np.pi*R[i]*dr
return masse
sigmaIni = repartSigma(110e6, 3e5)
sigmaFin= euler(0.1, 110e6, 3e5)
plt.figure()
plt.plot(Rdemi, sigmaIni, "r", label='0 years')
plt.plot(Rdemi, sigmaFin, "b", label=f'{Tmax/(365.25*24*3600)} ')
plt.title(f'{Tmax/(365.25*24*3600)} years')
plt.xlabel('Distance à Saturn (100000 km)')
plt.ylabel(r'$\Sigma$ (${10}^{4}kg.m^{-2}$)')
plt.xlim(100e6, 120e6)
plt.ylim(0, 8)
plt.legend()
plt.show()