Custom reference state for a non-dimensional Boussinesq MHD setup in a convective spherical shell
Custom reference state for a non-dimensional Boussinesq MHD setup in a convective spherical shell#
Custom reference state for a Boussinesq convective setup. We non-dimensionalize the MHD Boussinesq equations using the rotation period such that \([t]=1/\Omega_o\) is the timescale, the shell depth \([r]=r_o-r_i=L\) is the lengthscale (where \(r_o\) is the outer radius and \(r_i\) is the inner radius), \([u]=L\Omega_o\) is the velocity scale, and \([T]=\Delta T\) is the temperature scale*. Assuming that \(\Theta\) are the temperature perturbations, the non-dimensional Boussinesq equations can be written as:
\begin{equation} \dfrac{\partial\vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}+2\hat{z}\times\vec{u}=-\dfrac{\nabla P}{\rho_m}+{\rm{Ra^*}}\left(\dfrac{r_o}{r}\right)^n\Theta \hat{e}_r+\dfrac{1}{4\pi}\dfrac{E}{\rm{P_m}}(\nabla\times\vec{B})\times{\vec{B}}+E\nabla^2\vec{u} \end{equation}
\begin{equation} \nabla\cdot\vec{u}=0 \end{equation}
\begin{equation} \dfrac{\partial\Theta}{\partial t}+\vec{u}\cdot\nabla\Theta=\dfrac{E}{\rm{Pr}}\nabla^2\Theta \end{equation}
and
\begin{equation} \dfrac{\partial\vec{B}}{\partial t}-\nabla\times(\vec{u}\times\vec{B})=\dfrac{E}{\rm{P_m}}\nabla^2\vec{B}. \end{equation}
We have assumed that \(\nu\), \(\kappa\), and \(\eta\) are constants, \(\rho_m\) is the mean denisty and \(n\) is the gravity power (hence e.g. \(n=0\) for constant gravity). We also have the modified Rayleigh number Ra\(^*\) given by \begin{equation} {\rm{Ra^*}}=\dfrac{\alpha g_o \Delta T}{L\Omega_o^2}=\dfrac{\rm{Ra}}{\rm{Pr}}E^2, \end{equation} where \(\vec{g}(r)=g_o(r_o/r)^n\) , Pr=\(\nu/\kappa\) is the Prandtl number, \(E=\dfrac{\nu}{L^2\Omega_o}\) is the Ekman number and Ra\(=\alpha g_o\Delta T L^3/(\kappa\nu)\). Finally, the magnetic Prandtl number is P\(_m=\nu/\eta\). Then the corresponding fuctions \(f\) used here are:
\(f_1(r)\rightarrow 1\),
\(f_2(r)\rightarrow (r_o/r)^n\),
\(f_3(r)\rightarrow 1\),
\(f_4(r)\rightarrow 1\),
\(f_5(r)\rightarrow 1\),
\(f_6(r)\rightarrow 0\),
\(f_7(r)\rightarrow 1\),
\(f_8(r)\rightarrow 0\),
\(f_9(r)\rightarrow 0\),
\(f_{10}(r)\rightarrow 0\),
\(f_{14}(r)\rightarrow 0\),
\(f_{15}(r)\rightarrow 0\),
\(f_{16}(r)\rightarrow 1\)
and the constants \(c\) are:
\(c_1\rightarrow 2\),
\(c_2\rightarrow Ra E^2/Pr\),
\(c_3\rightarrow 1\),
\(c_4\rightarrow E/(4\pi{\rm{P_m})}\),
\(c_5\rightarrow E\),
\(c_6\rightarrow E/{\rm{Pr}}\),
\(c_7\rightarrow E/{\rm{P_m}}\),
\(c_8\rightarrow 0\),
\(c_9\rightarrow 0\),
\(c_{10}\rightarrow 0\).
*This is the relevant temperature scale for isothermal BCs – for fixed flux, fixed temperature BCs, the temperature scale should be something like \([T]=L dT_o/dr\).
[ ]:
#######################
import numpy
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot, draw, show
import os, sys
sys.path.insert(0, os.path.abspath('../../'))
import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook
[ ]:
# Grid Parameters
nr = 500 # Number of radial points
ri = 0.7e0 # Inner boundary of radial domain
ro = 1.0e0 # Outer boundary of radial domain
# radial grid
r=numpy.linspace(ri,ro,nr)
#aspect ratio
beta=ri/ro
# shell depth depending on the non-dimensionalization
d=1.e0
#non-dimensional r_i
ri_nd=beta*d/(1-beta)
#non-dimensional r_o
ro_nd=d/(1-beta)
#non-dimensional radial grid
radius1=numpy.linspace(ri_nd,ro_nd,nr)
print(ri_nd,ro_nd)
#print(radius1[0],radius1[nr-1])
[ ]:
ones = numpy.ones(nr,dtype='float64')
zeros = numpy.zeros(nr,dtype='float64')
# Here we define the reference state i.e. density, temperature,etc.
# For a classic RBC setup, this is the reference state to be used.
density = ones # density rho
dlnrho= zeros # dlnrho/dr
d2lnrho= zeros # d^2lnrho/dr^2
temperature=ones # temperature T
dlnt=zeros # dlnT/dr
pressure=ones # pressure P
entropy = zeros # entropy S, not used in Boussinesq -- set it = 0
gravity=zeros # gravity -- it is part of the buoyancy term in the non-dimensional momentum equation (see notes above)
hprofile=zeros # heating function (if we want one)
dsdr=zeros # dS/dr -- not useful in Boussinesq -- set it = 0
[ ]:
my_ref = rt.equation_coefficients(radius1)
[ ]:
## Here we define all the functions and constants that will be written in our data file and
## read by Rayleigh if we choose the custom reference state (=4). For more info, check notes above!
## Also, chech main_input_Boussinesq to see how to run a simulation with Rayleigh and this custom
## reference state ( "Boussinesq.dat" input file generated below!)
unity = numpy.ones(nr,dtype='float64')
gravity_power=0.0
buoy = (radius1[nr-1]/radius1)**gravity_power # buoyancy term calculation
my_ref.set_function(density,1) # denisty rho
my_ref.set_function(buoy,2) # buoyancy term
my_ref.set_function(unity,3) # nu(r) -- can be overwritten via nu_type in Rayleigh
my_ref.set_function(temperature,4) # temperature T
my_ref.set_function(unity,5) # kappa(r) -- works like nu
my_ref.set_function(hprofile,6) # heating function
my_ref.set_function(dlnrho,8) # dlnrho/dr
my_ref.set_function(d2lnrho,9) # d^2lnrho/dr^2
my_ref.set_function(dlnt,10) # dlnT/dr
my_ref.set_function(unity,7) # eta -- works like nu and kappa
my_ref.set_function(dsdr,14) # This is not used in Boussinesq -- set it = 0
# The constants can all be set/overridden in the input file
# NOTE that they default to ZERO, but we want
# most of them to be UNITY. These constants will explicitly depend on the non-dimensionalization chosen.
# The comments corresponding to each one of the constants are generic but here, we also specify
# what they are exactly in our example which is based on the non-dimensionalization used in this notebook!
## This is a non-magnetic example with Pr=1, E=0.001, and Ra=10^6, such that Ra*=1 !
my_ref.set_constant(2.0,1) # multiplies the Coriolis term, here it is: 2
my_ref.set_constant(1.0,2) # multiplies the buoyancy, here it is Ra*=Ra.E^2/Pr as defined above
my_ref.set_constant(1.0,3) # multiplies the pressure gradient
my_ref.set_constant(0.0 , 4) # multiplies the lorentz force, here it is: E/(4*pi*Pm)
my_ref.set_constant(0.001e0,5) # multiplies the viscosity, here it is: E
my_ref.set_constant(0.001e0,6) # multiplies the entropy diffusion (kappa), here it is: E/Pr
my_ref.set_constant(0.0,7) # multiplies eta in induction equation, here it is E/Pm
my_ref.set_constant(0.0,8) # multiplies viscous heating, here it is always 0, since we assume the Boussinesq approximation
my_ref.set_constant(1.0,9) # multiplies ohmic heating, here it is always 0, since we assume the Boussinesq approximation
my_ref.set_constant(1.0,10) # multiplies the heating, here it is 0, since we have assumed that there is no heating function!
my_ref.write('Boussinesq.dat') # Here we write our data file to be used to run our simulation with Rayleigh!
print(my_ref.fset)
print(my_ref.cset)
[ ]:
[ ]: