Solar Wind Acceleration

We consider the solar wind acceleration using 1D spherical MHD equations, $$ \frac{\partial}{\partial t}(r^2\rho)+\frac{\partial}{\partial t}(r^2\rho v_r)=0;\\ \frac{\partial}{\partial t}(r^2\rho v_r)+\frac{\partial}{\partial t}[r^2(\rho v_r^2+p)]=2rp-GM_\odot\rho;\\ p=\frac{2\rho T_0}{m_p} $$ Here, we assume that the temperature is constant and $T_0=2\times10^6 K$.

The analytical and steady state solution of this system is described as $$\left(\frac{v_r}{c_s}\right)^2-\ln\left(\frac{v_r}{c_s}\right)^2=4\ln\left(\frac{r}{r_c}\right)+4\frac{r_c}{r}+C.$$ When the constant $C=-3$, the solution passes the critical point ($v_r=c_s$ and $r=r_c$).
Here, $c_s$ is the sonic speed (constant) and $r_c=\frac{GM_\odot m_p}{4R_sT_0}$ where $R_s$ is the solar radius.

In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import math as mt
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
import matplotlib.animation as animation
import matplotlib as mpl

plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
mpl.rcParams['animation.html'   ]='jshtml'
In [2]:
nx=2000
nt=20000
isav=100
dx=0.01
dt=0.001*isav
In [3]:
f1 =open('data/rohst.d'  ,'rb')
f2 =open('data/uxhst.d'  ,'rb')

rr=1+dx*np.arange(nx)
ro =np.fromfile(f1, dtype='float64').reshape(-1,nx)
ux =np.fromfile(f2, dtype='float64').reshape(-1,nx)
In [4]:
plt.imshow(ux,origin='lower',aspect='auto',extent=(rr[0],rr[-1],0,nt*dt))
plt.colorbar(label='$v_r$')
plt.xlabel('$r/R_s$')
plt.ylabel('$time$')
plt.show()
In [5]:
u=np.linspace(1e-4,5,1000)
R,U=np.meshgrid(rr,u)
rc=5.75/2
fSW=U**2-2*np.log(U)-4*np.log(R/rc)-4*(rc/R)+3

def update_anim(it):
    
    fig.clf() #clear the figure
    
    ax=fig.subplots(1, 1) #add subplots

    fig.tight_layout() #reduce spacing around the figure

    ax.plot(rr,ux[it,:],'.-',markevery=20)
    ax.contour(R,U,fSW,[0])
    ax.set_ylim(-1,5)
    ax.set_xlabel('$x/R_s$')
    ax.set_ylabel('$v_r/c_s$')

fig=plt.figure(figsize=(6,3))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
#plt.subplots_adjust(wspace=0.1,hspace=0.01, bottom=0,top=1)

plt.close()
anim
Out[5]:
In [6]:
plt.plot(rr,ux[-1,:],'.',markevery=40)
plt.contour(R,U,fSW,[0])
plt.xlim(1,21)
plt.xlabel(r'$x/R_s$')
plt.ylabel(r'$v_r/c_s$')
plt.show()
In [ ]:
!jupyter nbconvert --to html --template tmp.tpl plotting_tips.ipynb