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.
Note: Be sure to answer all these questions at steady state, when the flow is not varying too much in time.
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
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
.
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.
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.
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
.
@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.
@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.
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.
@jit(nopython=True)
def advance_time_level(Qin,flux_x,flux_y,source):
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.
@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.
@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$.
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.
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.
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
if(MMask[i,j]==1):
Qin[i,j,mx] = 0.0
Qin[i,j,my] = 0.0
Qin[i,j,mz] = 0.0
Qin[i,j,en] = en_floor
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)
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
.
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.
Q=np.zeros((nx,ny,nQ))
MMask=np.zeros((nx,ny))
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))
That's for an empty domain.
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:
MMask[i,j]=1
Q[i,j,rh]=rh_floor
Q[:,:,en]=0.026/te0*Q[:,:,rh]/(aindex-1)
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.
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)
Let's define some plotting parameters
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!
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
Q1=advance_time_level(Q,flux_x,flux_y,sources) #advance the solution
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)
Q2=advance_time_level(Q1,flux_x,flux_y,sources)
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
fig.add_subplot(rows, columns, k)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
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 = fig.colorbar(im,fraction=0.046, pad=0.04)
cb.ax.set_ylabel('$T(K)$', rotation=-90)
im.set_interpolation('quadric')
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")
FloatProgress(value=975.686503855416, description='Progress', max=1000.0000000000001)
FloatProgress(value=0.0, description='Next output', max=50.00000000000001)
DONE
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.