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)

Time advance algorithms

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

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

            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

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

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:
            MMask[i,j]=1
            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
    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")
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 [ ]: