Cylindrical radiation

  • Notes follow.
  • This work was supported by the National Science Foundation under grant number CBET-1403403.
In [5]:
from IPython.display import HTML
HTML('<iframe src=https://ignite.byu.edu/nsf/FVDOM_cylindrical.pdf width=800 height=600></iframe>')
Out[5]:
In [1]:
import numpy as np
import scipy.constants   as const
import matplotlib.pyplot as plt
from   scipy.integrate import quad
import sympy as sp
%matplotlib inline

Setup

In [2]:
R  = 1
Tw = 300
σ  = const.Stefan_Boltzmann
I0 = σ*Tw**4/np.pi

kConst  = 1
TgConst = 1500
IbConst = σ*TgConst**4/np.pi

FVDOM

  • Solve both halves of a circle.
  • For use with ODT, which is peicewise symmetric
  • $-R <= r <= R$
In [3]:
nr = 401   #1001
 = 40    #100
 = 80    #200

######################################

rf = np.linspace(-R,R,nr+1)
r  = 0.5*(rf[1:]+rf[:-1])
κ  = np.ones(nr) * kConst
θ  = np.linspace(0,np.pi/2,)[::-1]
ψ  = np.linspace(0,np.pi,)
V  = np.ones(nr)
for k in range(nr):
    if rf[k+1]*rf[k] > 0:
        V[k] = 0.5*np.abs(rf[k+1]**2 - rf[k]**2)
    else:
        V[k] = 0.5*np.abs(rf[k+1]**2 + rf[k]**2)
Asn = np.ones(nr)
for k in range(nr):
    Asn[k] = (rf[k+1] - rf[k])


 = θ[1]-θ[0]
 = ψ[1]-ψ[0]
dr = r[1]-r[0]

T  = np.ones(nr)*TgConst
Ib = σ*T**4/np.pi
Iinf = σ*Tw**4/np.pi

I  = np.zeros((nr,,))
α  = np.zeros((,))
β  = np.zeros((,))
q  = np.zeros(nr)
Q  = np.zeros(nr)

for i in range():
    for j in range():
        α[i,j] = np.sin(θ[i]) * np.cos(ψ[j])
        β[i,j] = np.sin(θ[i]) * np.sin(ψ[j])

        
######################################

for i in range():
    
    #AAAAAAAAAAAAAAAAAA RHS, march <---- from BC, ψ goes clockwise from ----> dir
    
    j = -1
    
    k = nr-1
    I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*np.abs(rf[k+1])*Iinf) / \
               (-α[i,j]*np.abs(rf[k]) - α[i,j]*Asn[k] + V[k]*κ[k])
        
    for k in range(nr-2, int(nr/2), -1):
        I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*np.abs(rf[k+1])*I[k+1,i,j]) / \
                   (-α[i,j]*np.abs(rf[k]) - α[i,j]*Asn[k] + V[k]*κ[k])
            
    # Other j 
            
    for j in range(-2, int(/2)-1, -1):
        
        k = nr-1
        I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*np.abs(rf[k+1])*Iinf             + β[i,j]*Asn[k]/*I[k,i,j+1]) / \
                   (-α[i,j]*np.abs(rf[k]) - α[i,j]*Asn[k] + V[k]*κ[k]         + β[i,j]*Asn[k]/)
         
        for k in range(nr-2, int(nr/2), -1):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*np.abs(rf[k+1])*I[k+1,i,j]   + β[i,j]*Asn[k]*I[k,i,j+1]/) / \
                       (-α[i,j]*np.abs(rf[k]) - α[i,j]*Asn[k] + V[k]*κ[k]     + β[i,j]*Asn[k]/)
                
    #BBBBBBBBBBBBBBBBBBB LHS, march ----> from BC, ψ goes counterclockwise from ----> dir
    
    j = 0
    
    k = 0
    I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*np.abs(rf[k])*Iinf) / \
               ( α[i,j]*np.abs(rf[k+1]) + α[i,j]*Asn[k] + V[k]*κ[k])
        
    for k in range(1, int(nr/2)):
        I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*np.abs(rf[k])*I[k-1,i,j]) / \
                   ( α[i,j]*np.abs(rf[k+1]) + α[i,j]*Asn[k] + V[k]*κ[k])
            
    # Other j 
            
    for j in range(1, int(/2)):
        
        k = 0
        I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*np.abs(rf[k])*Iinf             + β[i,j]*Asn[k]/*I[k,i,j-1]) / \
                   ( α[i,j]*np.abs(rf[k+1]) + α[i,j]*Asn[k] + V[k]*κ[k]     + β[i,j]*Asn[k]/)
         
        for k in range(1, int(nr/2)):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*np.abs(rf[k])*I[k-1,i,j]       + β[i,j]*Asn[k]/*I[k,i,j-1]) / \
                       ( α[i,j]*np.abs(rf[k+1]) + α[i,j]*Asn[k] + V[k]*κ[k]     + β[i,j]*Asn[k]/)
                
    #                                                                         ^
    #                                                                        /                                                                           
    #CCCCCCCCCCCCCCCCCC RHS, march ----> from center, ψ goes clockwise from / dir
    
    for j in range(int(/2)-1,-1,-1):
        
        k = int(nr/2)
        I[k,i,j] = I[k-1,i,j] #+ (r[k]-r[k-1])*(I[k-1,i,j]-I[k-2,i,j])/(r[k-1]-r[k-2]) # extrapolation
        
        for k in range(int(nr/2)+1, nr):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*np.abs(rf[k])*I[k-1,i,j]       + β[i,j]*Asn[k]/*I[k,i,j+1]) / \
                       ( α[i,j]*np.abs(rf[k+1]) - α[i,j]*Asn[k] + V[k]*κ[k]     + β[i,j]*Asn[k]/)
                
    #                                                                             ^
    #                                                                              \                                                                           
    #DDDDDDDDDDDDDDDDDD LHS, march <---- from center, ψ goes counterclockwise from \ dir
                
    for j in range(int(/2), ):
        
        k = int(nr/2)
        I[k,i,j] = I[k+1,i,j] #+ (r[k]-r[k+1])*(I[k+1,i,j]-I[k+2,i,j])/(r[k+1]-r[k+2]) # extrapolation
        
        for k in range(int(nr/2)-1, -1, -1):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*np.abs(rf[k+1])*I[k+1,i,j]   + β[i,j]*Asn[k]/*I[k,i,j-1]) / \
                       (-α[i,j]*np.abs(rf[k]) + α[i,j]*Asn[k] + V[k]*κ[k]     + β[i,j]*Asn[k]/)
    
######################################

q = np.zeros(nr)
for k in range(nr):
    for i in range():
        for j in range():
            q[k] -= 4***I[k,i,j]*α[i,j]*np.sin(θ[i])

Q = np.zeros(nr)
Q[0]    = -(q[1] -q[0])  /(r[1] -r[0])
Q[1:-1] = -(q[2:]-q[:-2])/(r[2:]-r[:-2])
Q[-1]   = -(q[-1]-q[-2]) /(r[-1]-r[-2])

######################################

j1 = 3
j2 = -4 
plt.rc("font", size=14)
plt.plot(r/R,I[:,0,j1] /1000,'-')
plt.plot(r/R,I[:,0,j2]/1000,'-')
plt.plot(r/R,np.ones(len(r))*Ib[0]/1000);
plt.xlabel('r/R')
plt.ylabel(r'I (kw/m$^2$*st)')
plt.legend([r'$\psi=$'+'{:3.1f}'.format(ψ[j1]*180/np.pi)+r'$^o$', \
            r'$\psi=$'+'{:3.1f}'.format(ψ[j2]*180/np.pi)+r'$^o$', \
            r'I$_b$'], frameon=False)

plt.figure()
plt.plot(r/R,q/1000)
plt.xlabel('r/R')
plt.ylabel(r'q (kW/m$^2$)')

plt.figure()
plt.plot(r/R,Q/1000)
plt.xlabel('r/R')
plt.ylabel(r'Q (kW/m$^3$)');

Qdom2 = Q.copy()
rdom2 = r.copy()

FVDOM

  • Solve symmetric problem
  • $0 <= r <= R$
In [4]:
nr = 401
 = 40
 = 80

######################################

rf = np.linspace(0,R,nr+1)
r  = 0.5*(rf[1:]+rf[:-1])
κ  = np.ones(nr) * kConst
θ  = np.linspace(0,np.pi/2,)[::-1]
ψ  = np.linspace(0,np.pi,)
V  = np.ones(nr)
for k in range(nr):
    if rf[k+1]*rf[k] > 0:
        V[k] = 0.5*np.abs(rf[k+1]**2 - rf[k]**2)
    else:
        V[k] = 0.5*np.abs(rf[k+1]**2 + rf[k]**2)
Asn = np.ones(nr)
for k in range(nr):
    Asn[k] = (rf[k+1] - rf[k])

 = θ[1]-θ[0]
 = ψ[1]-ψ[0]
dr = r[1]-r[0]

T  = np.ones(nr)*TgConst
Ib = σ*T**4/np.pi
Iinf = σ*Tw**4/np.pi

I  = np.zeros((nr,,))
α  = np.zeros((,))
β  = np.zeros((,))
q  = np.zeros(nr)
Q  = np.zeros(nr)

for i in range():
    for j in range():
        α[i,j] = np.sin(θ[i]) * np.cos(ψ[j])
        β[i,j] = np.sin(θ[i]) * np.sin(ψ[j])
        
######################################

for i in range():
    
    #---------------------
    
    j = -1
    
    k = nr-1
    I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*rf[k+1]*Iinf) / \
               (-α[i,j]*rf[k] - α[i,j]*Asn[k] + V[k]*κ[k])
    
    for k in range(nr-2,-1,-1):
        I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*rf[k+1]*I[k+1,i,j]) / \
                   (-α[i,j]*rf[k] - α[i,j]*Asn[k] + V[k]*κ[k])
            
    #---------------------
    
    for j in range(-2,int(/2)-1,-1):
        
        k = nr-1
        I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*rf[k+1]*Iinf + β[i,j]*Asn[k]/*I[k,i,j+1]) / \
                   (-α[i,j]*rf[k] - α[i,j]*Asn[k] + V[k]*κ[k] + β[i,j]*Asn[k]/)
            
        for k in range(nr-2,-1,-1):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] - α[i,j]*rf[k+1]*I[k+1,i,j] + β[i,j]*Asn[k]/*I[k,i,j+1]) / \
                       (-α[i,j]*rf[k] - α[i,j]*Asn[k] + V[k]*κ[k] + β[i,j]*Asn[k]/)
                
    #---------------------
                
    for j in range(int(/2)-1,-1,-1):
        
        k = 0
        
        I[k,i,j] = I[0,i,-1-j]
            
        for k in range(1,nr):
            I[k,i,j] = (V[k]*κ[k]*Ib[k] + α[i,j]*rf[k]*I[k-1,i,j] + β[i,j]*Asn[k]*I[k,i,j+1]/) / \
                       ( α[i,j]*rf[k+1] - α[i,j]*Asn[k] + V[k]*κ[k] + β[i,j]*Asn[k]/)
            
    #break

######################################

q = np.zeros(nr)
for k in range(nr):
    for i in range():
        for j in range():
            q[k] -= 4***I[k,i,j]*α[i,j]*np.sin(θ[i])
            
r = np.hstack((-r[::-1],r))
q = np.hstack((-q[::-1], q))

Q = np.zeros(len(r))
Q[0]    = -(q[1] -q[0])  /(r[1] -r[0])
Q[1:-1] = -(q[2:]-q[:-2])/(r[2:]-r[:-2])
Q[-1]   = -(q[-1]-q[-2]) /(r[-1]-r[-2])


######################################

j1 = 3
j2 = -4 
II = np.hstack(( I[::-1,0,j1], I[:,0,j2] ))
plt.rc("font", size=14)
plt.plot(r/R,II/1000,'-')
plt.plot(r/R,II[::-1]/1000,'-')
plt.plot(r/R,np.ones(len(r))*Ib[0]/1000);
plt.xlabel('r/R')
plt.ylabel(r'I (kw/m$^2$*st)')
plt.legend([r'$\psi=$'+'{:3.1f}'.format(ψ[j1]*180/np.pi)+r'$^o$', \
            r'$\psi=$'+'{:3.1f}'.format(ψ[j2]*180/np.pi)+r'$^o$', \
            r'I$_b$'], frameon=False)

plt.figure()
plt.plot(r/R,q/1000)
plt.xlabel('r/R')
plt.ylabel(r'q (kW/m$^2$)')

plt.figure()
plt.plot(r/R,Q/1000)
plt.xlabel('r/R')
plt.ylabel(r'Q (kW/m$^3$)');

rdom1 = r.copy()
Qdom1 = Q.copy()

Ray tracing numerical integration

  • Exact solution, but implemented numerically.
  • I is computed using implicit Euler method.
In [5]:
def getκ(r):
    return kConst
    #return 1 + 0.2*np.cos(8*r)
def Tg(r):
    return TgConst
    #return 1200 + 900* np.cos(8*r)
def getIb(r):
    return σ*Tg(r)**4/np.pi

def I_IE(R,r,θ,ψ):
    
    ns = 200
    s = ( r*np.cos(np.pi-ψ) + np.sqrt(R**2 - r**2 + r**2*np.cos(np.pi-ψ)**2) ) 
    ns = int(ns*s/2/R) + 1
    s /= np.abs(np.sin(θ))
    ds = s/ns
    
    I  = I0
    sp = 0
    for i in range(ns):
        sp += ds
        #rl = np.sqrt( r**2 + (s-sp)**2*np.cos(ψ)**2 - 2*r*(s-sp)*np.cos(ψ)*np.cos(np.pi-θ) )
        rl = np.sqrt( r**2 + (s-sp)**2*np.sin(θ)**2 - 2*r*(s-sp)*np.sin(θ)*np.cos(np.pi-ψ) )
        κrl = getκ(rl)
        I = (I + ds*κrl*getIb(rl))/(1+κrl*ds)
    
    return I
In [6]:
nr = 200
 = 5
 = 10

 = np.pi/2/
 = np.pi/
θ  = np.linspace(/2, np.pi/2-/2, )
ψ  = np.linspace(/2, np.pi-/2, )
rq = np.linspace(0,R,nr+1)
r  = 0.5*(rq[1:]+rq[:-1])
 
q  = np.zeros(nr+1)

for ir in range(nr+1):
    for  in range():
        for  in range():
            s = ( rq[ir]*np.cos(np.pi-ψ[]) + np.sqrt(R**2 - rq[ir]**2 + rq[ir]**2*np.cos(np.pi-ψ[])**2) ) /np.abs(np.sin(θ[]))
            II = I_IE(R,rq[ir],θ[],ψ[])
            q[ir] -= 4*II*np.cos(ψ[])*np.sin(θ[])**2**  # 4* since ψ goes from 0:π/2 instead of 0:2π 

Q = -(q[1:] - q[:-1])/(rq[1]-rq[0])

QrayN = np.hstack((Q[::-1], Q))
rrayN = np.hstack((-r[::-1], r))

Ray tracing, analytic solution

  • Assumes constant $I_b$ and constant $k$
In [7]:
nr = 200
 = 5 
 = 10

 = np.pi/2/
 = np.pi/
θ  = np.linspace(/2, np.pi/2-/2, )
ψ  = np.linspace(/2, np.pi-/2, )
rq = np.linspace(0,R,nr+1)
r  = 0.5*(rq[1:]+rq[:-1])
 
q  = np.zeros(nr+1)

for ir in range(nr+1):
    for  in range():
        for  in range():
            s = ( rq[ir]*np.cos(np.pi-ψ[]) + np.sqrt(R**2 - rq[ir]**2 + rq[ir]**2*np.cos(np.pi-ψ[])**2) ) /np.abs(np.sin(θ[]))
            I = IbConst - (IbConst-I0)*np.exp(-kConst*s) 
            q[ir] -= 4*I*np.cos(ψ[])*np.sin(θ[])**2**  # 4* since ψ goes from 0:π/2 instead of 0:2π 

Q = -(q[1:] - q[:-1])/(rq[1]-rq[0])

QrayA = np.hstack((Q[::-1], Q))
rrayA = np.hstack((-r[::-1],r))

Exact solution

  • Assumes constant $I_b$ and constant $k$
  • Integrate with quad
In [8]:
def get_Q(r):
    
    def fI_1(θ,ψ):
        s = ( r*np.cos(np.pi-ψ) + \
              np.sqrt(R**2 - r**2 + r**2*np.cos(np.pi-ψ)**2) ) / \
            np.abs(np.sin(θ))
        dsdr = ( np.cos(np.pi-ψ) + (r*np.cos(np.pi-ψ)**2-r)/
               np.sqrt(R**2-r**2+r**2*np.cos(np.pi-ψ)**2) )/ \
               np.abs(np.sin(θ))
        return np.cos(ψ)*np.sin(θ)**2*kConst*(IbConst-I0)*np.exp(-kConst*s)*dsdr
               
    def fI_2(ψ):
        return quad(fI_1, 1E-6, np.pi/2, args=(ψ,))[0]
    return 4*quad(fI_2, 0, np.pi)[0] #  4* for symmetry in θ and ψ

#---------------------------------

r = np.linspace(0,R,100)
Q    = np.zeros(len(r))
for i in range(len(r)):
    Q[i] = get_Q(r[i])

rexct = np.hstack((-r[::-1], r))
Qexct = np.hstack((Q[::-1], Q))

Plot

In [9]:
plt.rc('font', size=14)
plt.plot(rexct/R, Qexct/1000, color='gray',  linewidth=12, label='exact') 
#plt.plot(rrayA/R, QrayA/1000, color='red', linewidth=4, label='RT Analytic') 
plt.plot(rrayN/R, QrayN/1000, color='orange',  linewidth=3, label='RT Numerical') 
plt.plot(rdom2/R, Qdom2/1000, color='black',  linewidth=2, label='FVDOM full')
#plt.plot(rdom1/R, Qdom1/1000, color='blue', linewidth=2, label='FVDOM symm')
plt.xlabel("r/R")
plt.ylabel(r"Q (kW/m$^3$)")
plt.xlim([-1,1])
plt.legend(frameon=False)
plt.grid();
In [ ]: