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:

$$\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}$$

$$\nabla\cdot\vec{u}=0$$

$$\dfrac{\partial\Theta}{\partial t}+\vec{u}\cdot\nabla\Theta=\dfrac{E}{\rm{Pr}}\nabla^2\Theta$$

and

$$\dfrac{\partial\vec{B}}{\partial t}-\nabla\times(\vec{u}\times\vec{B})=\dfrac{E}{\rm{P_m}}\nabla^2\vec{B}.$$

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 $${\rm{Ra^*}}=\dfrac{\alpha g_o \Delta T}{L\Omega_o^2}=\dfrac{\rm{Ra}}{\rm{Pr}}E^2,$$ 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

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)

print(ri_nd,ro_nd)

[ ]:

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
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)

[ ]:



[ ]: