# Midterm: due 11/12/2018 at 11:59pm¶

## Problem 01¶

The goal of this midterm is to understand and characterize a simple DeLaval nozzle. These nozzle can generate a supersonic flow from an initial expanding flow.

1. Using the usual definition of the Mach number, given by $$M_A=\frac{u}{c_s}$$ where $u$ is the flow speed and $c_s$ is the ion sound speed given by $c_s=\sqrt{\gamma e T_e/m_i}$, show that the $M_A$ increases along the nozzle axis $r=0$ for T=600K
2. Compare the pressure at the entrance and the exit of the nozzle for T=600 K.
3. Compute the thrust of the nozzle (i.e. the force at the exit) for T=300K, T=600K and T=1200K.

Note: Be sure to answer all these questions at steady state, when the flow is not varying too much in time.

In [1]:
import numpy as np
import math
from math import sqrt
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from numba import jit
from ipywidgets.widgets.interaction import show_inline_matplotlib_plots
from IPython.display import clear_output
from ipywidgets import FloatProgress
from IPython.display import display


## Definitions and constants¶

For practical reason we will define some constants that will be use to remember which component in the array correspond to which quantity, e.g. the temperature can be accessed using tp at the location i,j as Q[tp,i,j]. rh is the number density, mx is the momentum density in the x-direction, my and mz in the y and z directions. en is for the energy density, tp the temperature. nQ is the total number of variables, from rh to tp.

In [2]:
rh=0
mx=1
my=2
mz=3
en=4
tp=5
nQ=6


mu is the atomic number of the fluid, in this case aluminum. aindex is the adiabatic index.

In [3]:
mu=27.
aindex=1.1


Here we define our dimensional parameters. L0 is the geometrical scale length, t0 the time scale, n0 the density scale, v0 is the velocity scale, p0 the pressure scale, te0 the temperature scale. We also define a couple of floor quantities under which no density, temperature and pressure can go.

In [4]:
L0=1.0e-1
t0=1e-7
n0=6.0e22
v0=L0/t0
p0=mu*1.67e-27*n0*v0**2
te0=p0/n0/1.6e-19
rh_floor = 1.0E-8
T_floor = 0.00026/te0
rh_mult = 1.1
P_floor = rh_floor*T_floor


We now define our geometrical functions. They are used quite often so we will compile them inside the code to accelerate computations. To avoid writing two codes, one for slab symmetry and one for cylindrical symmetry, we wrote this code in for slab symmetry and in Cartesian coordinates, i.e. $(x,y,z)$. When using the code for problem with cylindrical symmetry then \begin{align} &x\rightarrow r\\ &y\rightarrow z\\ &z\rightarrow \phi\end{align}

The rc function computes the radius in cylindrical coordinates, which appears in source terms and in the conservation equations. Clearly, this is simply $x$ as mentioned above. That's why rc returns xc. To turn the code into slab symmetry, just have the function rc return 1.

In [5]:
@jit(nopython=True)
def xc(i):
return lxd + (i-ngu+1)/dxi

@jit(nopython=True)
def rc(i):
return xc(i)

@jit(nopython=True)
def yc(j):
return lyd + (j-ngu+1)/dyi


We now add radius and angle function for axi-symmetric distribution in slab geometries, not ot be confused with rc above. Note the the letter c, for center, is now gone.

In [6]:
@jit(nopython=True)
def r(i,j):
return math.sqrt((xc(i)-(lxd+lxu)/2)**2 + (yc(j)-(lyd+lyu)/2)**2)

@jit(nopython=True)
def rz(i,j):
return sqrt(xc(i)**2+yc(j)**2)

@jit(nopython=True)
def theta(i,j):
return math.atan(yc(j)-(lyd+lyu)/2,xc(i)-(lxd+lxu)/2)


We now compute the sources for all our conservation equations.

In [7]:
def get_sources(Qin):
sourcesin=np.zeros((nx,ny,nQ))
for j in range(ngu, ny-ngu):
for i in range(ngu, nx-ngu):

rbp = 0.5*(rc(i+1) + rc(i))
rbm = 0.5*(rc(i) + rc(i-1))

dn = Qin[i,j,rh]
dni=1./dn
vx = Qin[i,j,mx]*dni
vy = Qin[i,j,my]*dni
vz = Qin[i,j,mz]*dni

T = (aindex - 1)*(Qin[i,j,en]*dni - 0.5*(vx**2 + vy**2 + vz**2))
if(T <  T_floor):
T = T_floor
Qin[i,j,tp]=T
sourcesin[i,j,rh] = 0

P = dn*T

sourcesin[i,j,mx] =  (Qin[i,j,mz]*vz + P)
sourcesin[i,j,mz] =  - Qin[i,j,mz]*vx
sourcesin[i,j,en] = 0
return sourcesin


Here do the time advance of the conservation equations. The flux are used here as seen in Eq. $(4.26)$. The flux are computed later.

In [8]:
@jit(nopython=True)
Qout=np.zeros((nx,ny,nQ))
dxt = dxi*dt
dyt = dyi*dt

for j  in range(ngu, ny-ngu):
for i  in range(ngu, nx-ngu):
rbp = 0.5*(rc(i+1) + rc(i))
rbm = 0.5*(rc(i) + rc(i-1))
rci = 1./rc(i)

Qout[i,j,rh] = Qin[i,j,rh]*rc(i) - dxt*(flux_x[i,j,rh]*rbp-flux_x[i-1,j,rh]*rbm) - dyt*(flux_y[i,j,rh]-flux_y[i,j-1,rh])*rc(i)
Qout[i,j,mx] = Qin[i,j,mx]*rc(i) - dxt*(flux_x[i,j,mx]*rbp-flux_x[i-1,j,mx]*rbm) - dyt*(flux_y[i,j,mx]-flux_y[i,j-1,mx])*rc(i) + dt*source[i,j,mx]
Qout[i,j,my] = Qin[i,j,my]*rc(i) - dxt*(flux_x[i,j,my]*rbp-flux_x[i-1,j,my]*rbm) - dyt*(flux_y[i,j,my]-flux_y[i,j-1,my])*rc(i) + dt*source[i,j,my]
Qout[i,j,mz] = Qin[i,j,mz]*rc(i) - dxt*(flux_x[i,j,mz]*rbp-flux_x[i-1,j,mz]*rbm) - dyt*(flux_y[i,j,mz]-flux_y[i,j-1,mz])*rc(i) + dt*source[i,j,mz]
Qout[i,j,en] = Qin[i,j,en]*rc(i) - dxt*(flux_x[i,j,en]*rbp-flux_x[i-1,j,en]*rbm) - dyt*(flux_y[i,j,en]-flux_y[i,j-1,en])*rc(i) + dt*source[i,j,en]

Qout[i,j,rh:en+1] = Qout[i,j,rh:en+1]*rci
return Qout


We now write the function that computes the flux flx and the freezing speeds cfx along the x-direction. These speeds correspond to the largest speed found in the hyperbolic equation at the given location $(x,y)$. Using freezing speed avoids using a Riemann solver.

In [9]:
@jit(nopython=True)
def calc_flux_x(Qx):
cfx=np.zeros((nx,nQ))
flx=np.zeros((nx,nQ))

for i in range(0, nx):

dn = Qx[i,rh]
dni = 1./dn
vx = Qx[i,mx]*dni
vy = Qx[i,my]*dni
vz = Qx[i,mz]*dni

P = (aindex - 1)*(Qx[i,en] - 0.5*dn*(vx**2 + vy**2 + vz**2))
if(P < P_floor):
P = P_floor

flx[i,rh] = Qx[i,mx]
flx[i,mx] = Qx[i,mx]*vx + P
flx[i,my] = Qx[i,my]*vx
flx[i,mz] = Qx[i,mz]*vx
flx[i,en] = (Qx[i,en] + P)*vx

asqr = sqrt(aindex*P*dni)
vf1 = sqrt(vx**2 + aindex*P*dni)

cfx[i,rh] = vf1
cfx[i,mx] = vf1
cfx[i,my] = vf1
cfx[i,mz] = vf1
cfx[i,en] = vf1
return cfx,flx


We now write the function that computes the flux fly and the freezing speeds cfy along the y-direction.

In [10]:
@jit(nopython=True)
def calc_flux_y(Qy):
cfy=np.zeros((ny,nQ))
fly=np.zeros((ny,nQ))

for j in range(0, ny):

dn = Qy[j,rh]
dni = 1./dn
vx = Qy[j,mx]*dni
vy = Qy[j,my]*dni
vz = Qy[j,mz]*dni

P = (aindex - 1)*(Qy[j,en] - 0.5*dn*(vx**2 + vy**2 + vz**2))
if(P < P_floor):
P = P_floor

fly[j,rh] = Qy[j,my]
fly[j,mx] = Qy[j,mx]*vy
fly[j,my] = Qy[j,my]*vy + P
fly[j,mz] = Qy[j,mz]*vy
fly[j,en] = (Qy[j,en] + P)*vy

asqr = sqrt(aindex*P*dni)
vf1 = sqrt(vy**2 + aindex*P*dni)

cfy[j,rh] = vf1
cfy[j,mx] = vf1
cfy[j,my] = vf1
cfy[j,mz] = vf1
cfy[j,en] = vf1
return cfy,fly


We now smooth out the different fluxes using a total diminishing variation scheme (or TVD]. Basically, this schemes applies numerical viscosity to smooth out sharp transitions caused by shocks or sharp interfaces. Typically, sharp transitions cause numerical oscillations (called Gibbs phenomenon). Note that this scheme is conservative, in the sense that it does conserve all global quantities $n, n\vec u,\epsilon$.

In [11]:
def tvd2(Qin,n,ff,cfr):
sl =1 #use 0.75 for first order
flux2=np.zeros((n,nQ))
wr = cfr*Qin + ff
wl = cfr*Qin - ff

fr = wr
fl = np.roll(wl,-1,axis=0)

dfrp = np.roll(fr,-1,axis=0) - fr
dfrm = fr - np.roll(fr,+1,axis=0)
dflp = fl - np.roll(fl,-1,axis=0)
dflm = np.roll(fl,+1,axis=0) - fl
dfr=np.zeros((n,nQ))
dfl=np.zeros((n,nQ))

for l in range(nQ):
for i in range(n):
if(dfrp[i,l]*dfrm[i,l] > 0) :
dfr[i,l] = dfrp[i,l]*dfrm[i,l]/(dfrp[i,l] + dfrm[i,l])
else:
dfr[i,l] = 0
if(dflp[i,l]*dflm[i,l] > 0) :
dfl[i,l] = dflp[i,l]*dflm[i,l]/(dflp[i,l] + dflm[i,l])
else:
dfl[i,l] = 0
flux2 = 0.5*(fr - fl + sl*(dfr - dfl))
return flux2


We truly compute all the fluxes here.

In [12]:
def get_flux(Qin):
flux_x=np.zeros((nx,ny,nQ))
flux_y=np.zeros((nx,ny,nQ))
fx=np.zeros((nx,nQ))
cfx=np.zeros((nx,nQ))
ffx=np.zeros((nx,nQ))
fy=np.zeros((ny,nQ))
cfy=np.zeros((ny,nQ))
ffy=np.zeros((ny,nQ))
for j  in range(0, ny):
cfx,ffx=calc_flux_x(Qin[:,j,:])
flux_x[:,j,:]=tvd2(Qin[:,j,:],nx,ffx,cfx)
for i  in range(0, nx):
cfy,ffy=calc_flux_y( Qin[i,:,:])
flux_y[i,:,:]=tvd2(Qin[i,:,:],ny,ffy,cfy)
return flux_x,flux_y


We now limit flows to avoid minor instabilities caused by round errors.

In [13]:
def limit_flow(Qin):

en_floor = rh_floor*T_floor/(aindex-1)

for j in range(ngu-1, ny-ngu+1):
for i in range(ngu-1, nx-ngu+1):

if(Qin[i,j,rh]  <=  rh_floor):
Qin[i,j,rh] = rh_floor
Qin[i,j,mx] = 0.0
Qin[i,j,my] = 0.0
Qin[i,j,mz] = 0.0
Qin[i,j,en] = en_floor

Qin[i,j,mx] = 0.0
Qin[i,j,my] = 0.0
Qin[i,j,mz] = 0.0
Qin[i,j,en] = en_floor


### Grid definition¶

In [14]:
ngu=2
nx=30
ny=nx*4

lxu=5e-1/L0 #note that 15e-3 is in m while 15e-3/L0 is dimensionless
lxd=0

lyd=0
lyu = lyd + (ny-2*ngu)*(lxu-lxd)/(nx-2*ngu)

dxi = (nx-2*ngu)/(lxu-lxd)
dyi = (ny-2*ngu)/(lyu-lyd)


### Time step computation¶

In [15]:
def get_min_dt(Q):
cfl=0.025 #this is where the CFL condition is changed. Keep the rest identical
vmax =sqrt(aindex*1.6e-19*T_floor*te0/(mu*1.6e-27))/v0
for j in range(ny):
for i in range(nx):
dn = Q[i,j,rh]
dni=1./dn
vx = Q[i,j,mx]*dni
vy = Q[i,j,my]*dni
vz = Q[i,j,mz]*dni

T = (aindex - 1)*(Q[i,j,en]*dni - 0.5*(vx**2 + vy**2 + vz**2))
if(T <  T_floor):
T = T_floor
cs=sqrt(aindex*1.6e-19*T*te0/(mu*1.6e-27))/v0
if(vmax < cs  and T > rh_mult*T_floor):
vmax = cs
v = sqrt(vx**2+vy**2+vz**2)
if(vmax < v  and  Q[i,j,rh] > rh_mult*rh_floor):
vmax = v
return min(1./(vmax*dxi),1./(vmax*dyi))*cfl


Below we defined the initial ti and final tf times. n_out corresponds to to the total number of outputs. The rest are variable initialization. Note the most important one, the time t.

In [16]:
ti = 00.0E-2
tf = 100e-6/t0
n_out=20

n_steps=0
t_count=0
dt=0.
t = ti


And now all our arrays are defined here. Note that we use a mask MMask for parts of the simulation space were we want to preclude motion. That will come handy later in the course.

In [17]:
Q=np.zeros((nx,ny,nQ))
Q1=np.copy(Q)
Q2=np.copy(Q)
sources=np.zeros((nx,ny,nQ))
flux_x=np.zeros((nx,ny,nQ))
flux_y=np.zeros((nx,ny,nQ))


### Initial conditions¶

That's for an empty domain.

In [18]:
Q[:,:,rh]=rh_floor
sigma=(lyu-lyd)/16
for j in range(ny):
rx=yc(j)-(lyu-lyd)/4
r_nozzle=1-math.exp(-rx**2/sigma**2)
for i in range(nx):
if rc(i)>(lxu-lxd)/3+((lxu-lxd)/4)*r_nozzle and yc(j)<(lyu-lyd)/2:
Q[i,j,rh]=rh_floor
Q[:,:,en]=0.026/te0*Q[:,:,rh]/(aindex-1)


### Boundary conditions¶

Below we define our boundary conditions. In this example, we do not use the time t and focus only on steady state solutions. The first block define the default boundary conditions for the whole domain. Note that the conditions for the y-direction are pretty simple. These boundary conditions correspond to an open boundary. everything can go through it.

Boundary conditions are a bit more complex for the x-direction. We are using a cylindrical coordinate system with open boundary conditions on the right but reflecting boundary conditions on the left. That's because nothing can pass through the axis of symmetry $r=0$.

At the end we have the boundary conditions that make up our problem, e.g. material free falling toward the axis of our coordinate system. Note that we are using very large density because the author works in high energy density plasmas. the initial parameters L0, t0 and n0 can be changed to accommodate other regimes. But, chances are that the CFL scaler cflm may have to be adjusted.

In [19]:
def set_bc(Qin,t):
if ((rc(0)==rc(1))or(lxd!=0)):
#left open boundary conditions at x=lxu
Qin[1,:,:] = Qin[2,:,:]*rc(2)/rc(1)
Qin[0,:,:] = Qin[1,:,:]*rc(1)/rc(0)
else:
#left reflecting boundary conditions at x=lxd when the symmetry axis is located there
for l in range(ngu):
for k in range(nQ):
if(k == mx):
Qin[l,:,k] = -Qin[2*ngu-l-1,:,k]
if (k == mz) :
Qin[l,:,k] = -Qin[2*ngu-l-1,:,k]
if (k == my or k == rh or k == en) :
Qin[l,:,k] = Qin[2*ngu-l-1,:,k]
#right open boundary conditions at x=lxu
Qin[nx-ngu,:,:] = Qin[nx-ngu-1,:,:]*rc(nx-ngu-1)/rc(nx-ngu)
Qin[nx-ngu+1,:,:] = Qin[nx-ngu,:,:]*rc(nx-ngu)/rc(nx-ngu+1)
#bottom open boundary conditions at y=lyd
Qin[:,1,:] = Qin[:,2,:]
Qin[:,0,:] = Qin[:,1,:]
#top open boundary conditions at ly=lyu
Qin[:,ny-ngu,:] = Qin[:,ny-ngu-1,:]
Qin[:,ny-ngu+1,:] = Qin[:,ny-ngu,:]

T_bc=600/11604.525 #conversion from kelvin to volt
for i in range (nx):
if (rc(i)<lxu/2):
for k in range(ngu):
Q[i,k,rh]=6e22/n0
Q[i,k,en]=T_bc/te0*Q[i,k,rh]/(aindex-1)


### Computation and outputs¶

Let's define some plotting parameters

In [20]:
matplotlib.rcParams.update({'font.size': 22})
xi = np.linspace(lxd*L0, lxu*L0, nx-2*ngu-1)
yi = np.linspace(lyd*L0, lyu*L0, ny-2*ngu-1)
progress_bar = FloatProgress(min=ti, max=tf,description='Progress')
output_bar = FloatProgress(min=0, max=(tf-ti)/n_out,description='Next output')
columns = 3
rows = 2
box=np.array([lxd,lxu,lyd,lyu])*L0


And here we go!

In [21]:
while t<tf:
dt=get_min_dt(Q) #compute the optimal dt

set_bc(Q,t) #compute the boundary cnditions
sources=get_sources(Q) #compute the sources
flux_x,flux_y=get_flux(Q) #and ten the fluxes
limit_flow(Q1) #limit flows that are too fast

set_bc(Q1,t+dt) #let's do it again
sources=get_sources(Q1)
flux_x,flux_y=get_flux(Q1)
limit_flow(Q2)

Q=0.5*(Q+Q2) #here we do a simple second average
limit_flow(Q)

n_steps+=1 #time step increased by 1
t_count+=dt #and the time by dt

progress_bar.value=t #we have defined some progress bars cause python is slow
output_bar.value+=dt # we need to monitor our progress and reduce resolution if it takes too long

if ((t==ti) or (t_count>(tf-ti)/n_out)): #everything below is plotting....
fig=plt.figure(figsize=(20, 19))
for i in range(1, rows+1):
for j in range(1, columns+1):
k=j+(i-1)*columns
if (k==1):
data= np.flip(np.transpose(Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]*n0),0)
im = plt.imshow(data, cmap=plt.cm.nipy_spectral, norm=colors.LogNorm(vmin=data.min(), vmax=data.max()),extent=box)
cb.ax.set_ylabel('$\log_{10}n_i(m^{-3})$', rotation=-90)
if (k==2):
data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mx]**2
data+=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,my]**2
data+=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mz]**2
data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2
data=np.sqrt(data)*L0/t0
data= np.flip(np.transpose(data),0)
im = plt.imshow(data, extent=box)
cb.ax.set_ylabel('$u(m/s)$', rotation=-90)
if (k==4):
data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mx]**2
data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2
data=np.sqrt(data)*L0/t0
data= np.flip(np.transpose(data),0)
im = plt.imshow(data, extent=box)
cb.ax.set_ylabel('$u_r(m/s)$', rotation=-90)
if (k==5):
data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,my]**2
data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2
data=np.sqrt(data)*L0/t0
data= np.flip(np.transpose(data),0)
im = plt.imshow(data, extent=box)
cb.ax.set_ylabel('$u_z(m/s)$', rotation=-90)
if (k==6):
data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mz]**2
data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2
data=np.sqrt(data)*L0/t0
data= np.flip(np.transpose(data),0)
im = plt.imshow(data, extent=box)
cb.ax.set_ylabel('$u_\phi(m/s)$', rotation=-90)
if (k==3):
data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,tp]*te0*11604.525
data= np.flip(np.transpose(data),0)
im = plt.imshow(data, cmap=plt.cm.jet,  extent=box)
cb.ax.set_ylabel('$T(K)$', rotation=-90)
cb.ax.yaxis.set_label_coords(6.5, 0.5)
plt.gca().axes.get_yaxis().set_visible(False)
plt.gca().axes.get_xaxis().set_visible(False)
if (j==1):
plt.ylabel('z(m)', rotation=90)
plt.gca().axes.get_yaxis().set_visible(True)
if (i==rows):
plt.xlabel('r(m)', rotation=0)
plt.gca().axes.get_xaxis().set_visible(True)
plt.tight_layout()
clear_output()
output_bar.value=0
display(progress_bar) # display the bar
display(output_bar) # display the bar
t_count=0
show_inline_matplotlib_plots()
t=t+dt

#we are now done. Let's turn off the progress bars.
output_bar.close()
del output_bar
progress_bar.close()
del progress_bar

print("DONE")

DONE


## Problem 2¶

Using the binary collision function of HW05, plot the distribution function of 10 positrons and 10 electrons locked inside a box of size 2mmx2mm, randomly distributed, for each of the following initial parameters: $v_p=v_e=1000\,m/s$, all with random directions

Note that we will suppose that there is no displacement along the Z-direction and that all particles reflect back inside the box when they hit the walls of the box.

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: