CHAPTER 21
21.1 Consider the (highly artificial) absorption coefficient of Problem 11.22. Find narrow band averages for the
absorption coefficient and the transmissivity using Monte Carlo integration (use mcint.f90 or write your
own code). Compare with answers from Problem 11.22.
¯
κη=1
κηdη=1
10.5d
Z
A[η+3h(η)]dη
¯
tη=1
∆ηZ
∆η
e−κηLdη=1
10d
10.5d
Z
.5d
e−A[η+3h(η)]dη.
Modifying mcintegral.f90 as shown below (for the transmissivity case only), leads to
a=0.5*d
b=10.5*d
varmax=0.002
.
.
.
FUNCTION F(x)
USE PARAM
478 RADIATIVE HEAT TRANSFER
10000 2.2218E-01 1.6432E-03 0.74
40000 2.2312E-01 4.9700E-04 0.22
AdL= 0.50
no. of bundles integral std dev rel.err(%)
20000 1.1059E-01 1.2477E-02 11.28
80000 8.2596E-02 3.1078E-03 3.76
320000 7.5506E-02 7.8071E-04 1.03
1280000 7.3621E-02 1.8853E-04 0.26
.
.
.
AdL= 3.00
no. of bundles integral std dev rel.err(%)
20000 9.9414E-05 5.6746E-06 5.71
80000 8.8309E-05 1.9382E-06 2.19
320000 8.4443E-05 6.3084E-07 0.75
1280000 8.3445E-05 3.7077E-07 0.44
5120000 8.2644E-05 1.3654E-07 0.17
and in graphical form (including results from Problem 11.22):
AdL
trex, trgr, trmc
10-3
10-2
10-1
trex
trgr
trmc
CHAPTER 21 479
21.8 Consider an absorbing-scattering slab irradiated by a short-pulsed laser, as described in Example 19.3. Prepare
a transient Monte Carlo code to predict the flux exiting the slab as a function of time into either direction
(transmissivity and reflectivity).
Transmissivity, q/q0
0.005
0.01
0.015
0.02
0.025
slab reflectivity
slab transmissivity
τL= 5; ω= 0.9
beta=5. ! 1/m, extinction coefficient
omga=1. ! scattering albedo
tps=0.3 ! nondimensional equivalent pulse duration, beta*c0*t
N=100000
dt=tmax/numt
IF(kap<1.e-8) kap=1.e-8
yi=0.
sx=0.
sz=1.
! time of emission (0<tbm<tps)
tbm=ran1(idum)*tps
! total distance until absorption
labs=-alog(ran1(idum))/kap
! location of next scattering (or the absorption) event
xe=xi+lend*sx
ye=yi+lend*sy
ze=zi+lend*sz
! check whether bundle leaves domain
tbm=tbm+lbot*beta ! nondim time required to do that
it=tbm/dt+1
IF(it>numt) it=numt+1
trmsv(i,it)=trmsv(i,it)+dq
CYCLE
sz=1.-2.*ran1(idum)
sinth=sqrt(1.-sz*sz)
psi=2.*pi*ran1(idum)
sx=sinth*cos(psi)
sy=sinth*sin(psi)
GOTO 10
ENDDO !in
trmsvav(it)=trmsvav(it)/numsmpl
reflcav(it)=reflcav(it)/numsmpl
stddevt(it)=0.
stddevr(it)=0.
DO i=1,numsmpl
DO it=1,numt
tmit=(it-0.5)*dt
write(*,3) tmit,trmsvav(it),stddevt(it),reflcav(it),stddevr(it)
ENDDO
2 format(’no. of bundles =’,i10/’ time trans std dev’&
N=2*N
ns1=numsmplhf+1
GOTO 5
20 write(8,1)
1 format(’variables = “time”, “trans”, “stddevt”, “reflc”, “stddevr”’/’zone’)
END