ccc
ccc Title: A validated model for the dispersion of exhaust emissions from aircraft in their take-off roll
ccc Authors: A. Graham1, M. Bennett*2, S. Christie2 and D.W. Raper2
ccc
ccc * To whom correspondence should be addressed: m.bennett@mmu.ac.uk
ccc
ccc Affliation 1: Norwegian Centre of Offshore Wind Energy, Postboks 6031, N-5892 Bergen, Norway angus.graham@uni.no
ccc Affiliation 2 Centre for Aviation, Transport and the Environment, Manchester Metropolitan University, Manchester M1 5GD, UK
ccc
ccc Copyright. Manchester Metropolitan University 2018.
ccc
ccc
ccc The following is a FORTRAN subroutine which codes the mathematical model described in the above paper
ccc Specimen output files can be found in: TOdispResScenA.txt; and TOdispResScenB.txt.
ccc
ccc
ccc
ccc The source code is F77 including the use of mathematical intrinsic functions such as hyperbolic, which most code libraries should support.
ccc It is fully commented so as to be reasonably easily adaptable by future users.
ccc
ccc
real*8 fricv,wdir,invMO,F0max,MTOM,wingar,dum
print *,'friction velocity (m s-1) ?'
read *,fricv
print *,'wind dir rel to aircraft hding (+- 180 degs) ?'
read *,wdir
print *,'inverse Monin Obukhov length (m-1) ?'
read *,invMO
print *,'boundary-layer height (m) ?'
read *,dum
iblht = idint(dum)
if (iblht.lt.1) then
print *,'value too low: prog ends'
stop
endif
print *,'top static thrust from pooled engines (kN) ?'
read *,dum
F0max = 1000*dum
print *,'maximum takeoff mass (tonnes) ?'
read *,dum
MTOM = 1000*dum
print *,'basic wing planform area (m2) ?'
read *,wingar
call getTOd(fricv,wdir,invMO,iblht,F0max,MTOM,wingar)
end
************************************************************************
subroutine getTOd(fricv,wdir,invMO,iblht,F0max,MTOM,wingar)
real*8 zero,sig,sig1,dtor,sixth,quart,third,twothd,quarpi,one,
+ cthrhf
parameter (zero=0,sig=1e-9,sig1=1e-3,quart=.25,one=1,
+ sixth=one/6,third=2*sixth,twothd=2*third,cthrhf=1.5)
real*8 kmF,kmM,rollco,gacc,corot,CLmax,rho,therc1,therc2,therc3,
+ therc4,tau,xres,coke,cokF,corat,therc5,therc6,therc7,
+ therc8
parameter (kmF=.85,kmM=.90,rollco=.025,gacc=9.8,corot=1.2,
+ CLmax=2.25,rho=1.25,therc1=11.5,therc2=59,therc3=75,
+ therc4=5,tau=80,xres=.5,coke=.3,cokF=.7,therc5=34,
+ therc6=third,therc7=1.5,therc8=1.4)
ccc kmF .. fraction of max thrust characterisrtically
ccc used during takeoff
ccc kmM .. fraction of max takeoff weight
ccc characteristically holding for departing
ccc aircraft
ccc rollco .. coefficient of friction of aircraft tyres on
ccc runway
ccc gacc .. free acceleration due to gravity (m s-2)
ccc corot .. airspeed of aircraft when upward rotation of
ccc aircraft is begun, expressed as a factor of
ccc stall speed
ccc CLmax .. maximum attainable lift coefficient of aircraft
ccc at takeoff flap setting
ccc rho .. density of ambient air at ground level (kg m-3)
ccc therc1 .. log_e (z_ref / z_r) / k_vK, where z_ref is
ccc reference height of 5 m, z_r is representative
ccc ambient roughness length of 5 cm and k_vK is
ccc von Karman's constant, 0.4
ccc therc2 .. alpha_s * z_ref / k_vK (m), where alpha_s is
ccc 4.7 (in stable ambient stratification, shear at
ccc height z_ref, as normalised on
ccc fricv / (k_vK * z_ref), is equal to
ccc 1 + alpha_s * z_ref / L_M, L_M being Monin-
ccc Obukhov length)
ccc therc3 .. alpha_c * z_ref (m), where alpha_c is 15
ccc (in ambient convection, shear at height z_ref,
ccc as normalised on fricv / (k_vK * zref), is equal
ccc to 1 / (1 - alpha_c * z_ref / L_M) ** [1 / 4])
ccc therc4 .. 2 / k_vK
ccc tau .. F0 / B_F (s), where B_F is buoyant force
ccc imparted to exhaust pool in unit release time
ccc xres .. resolution of grid in direction of release over
ccc which exhaust ages are solved for during ground
ccc run, prior to upward rotation of aircraft (m)
ccc coke .. rate of change (as total derivative) of square
ccc root of cross-sectional area of exhausts, as
ccc normalised on their excess speed in direction
ccc of release
ccc cokF .. factor by which excess momentum of exhausts in
ccc in direction of release is reduced as result of
ccc drag losses at ground
ccc therc5 .. value in stable ambient stratification of k_rho,
ccc equal to fractional rate of change with height
ccc of potential density at height, L_M, as
ccc normalised on -(fricv / L_M) ** 2 / gacc (a
ccc thermal layer is identified between ground and
ccc height L_M in which density gradient is treated
ccc as constant, and assigned a value as in reality
ccc holding at height L_M)
ccc therc6 .. in stable ambient stratification, potential-
ccc density gradient in upper thermal layer holding
ccc between heights L_M and iblht (where L_M',
+ iTfir,' to continue) ?'
read *,iTlas
if (iTlas.gt.iTfir) then
iTold = iTfir
ccc when some of ground run prior to rotation
ccc remains to be processed
120 if (iTold.lt.iTrot) then
ccc set start grid time to be processed
mTolgs = iTfir*ngsps + 1
ccc set end grid time to be processed
mTgs = min(iTrot,iTlas)*ngsps
ccc step through grid times
do 500 iTgs=mTolgs,mTgs
ccc update cX
cX = (mxnew - mx)*xresnd - delXnd
if (cX.lt.-xresnd-sig .or. cX.gt.sig)
+ then
print *,'algorithm fails: prog ends'
stop
endif
ccc update furthest grid location in direction
ccc of release from source to have been
ccc processed
mx = mxnew
m = mx + 1
ccc initialise youngest appropriate gridded age
iolagt = 1
Tnd = Tnd + Tresnd
ccc update MBnd and dMBTnd in case of
ccc exhausts released before last grid time
ccc processed, according to exhaust
ccc characteristics at each grid location,
ccc ix, at last grid time
do 140 ix=1,m
ccc calc oldest appropriate gridded age
if (ix.le.mx) then
n = idint(agend(ix)/Tresnd)
else
n = iTgs - 1
endif
if (n.ge.iTgs) then
print *,'algorithm fails: prog ends'
stop
endif
do 130 i=iolagt,n
j = iTgs - i
if (j.lt.1 .or. j.ge.iTgs) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
if (ix.le.mx) then
dum = znd(ix)
else
dum = Zhnd
endif
ccc obtain mean height exhausts have
ccc risen buoyantly to, as
ccc nondimensionalised on Monin-Obukhov
ccc length
dum = dum*zsndMO
ccc when centroid in lower ambient
ccc thermal layer
if (dabs(dum).lt.one) then
del = karho
ccc when centroid in upper ambient
ccc thermal layer (in stable
ccc stratification)
else if (dum.ge.one) then
del = karhou
else
del = zero
endif
dum1 = del*Tresnd
ccc when stable ambient
ccc stratification is absent
if (invMO.le.zero) then
if (dum1.lt.sig) then
dum = Tresnd
else
dum = dsinh(dum1)/del
endif
dum1 = dcosh(dum1)
k = 1
else
if (dum1.lt.sig) then
dum = Tresnd
else
dum = dsin(dum1)/del
endif
dum1 = dcos(dum1)
k = -1
endif
del1 = dum*dMBTnd(j) + dum1*MBnd(j)
dMBTnd(j) = dum1*dMBTnd(j) +
+ k*(del**2)*dum*MBnd(j)
MBnd(j) = del1
130 continue
ccc update iolagt as necessary
if (ix.le.mx) iolagt = n + 1
140 continue
dum = karho*Tresnd
ccc update MBnd and dMBTnd in case of
ccc exhausts released at last grid time
ccc processed
if (dum.lt.sig) then
MBnd(iTgs) = Tresnd
else if (invMO.le.zero) then
MBnd(iTgs) = dsinh(dum)/karho
else
MBnd(iTgs) = dsin(dum)/karho
endif
if (invMO.le.zero) then
dMBTnd(iTgs) = dcosh(dum)
else
dMBTnd(iTgs) = dcos(dum)
endif
ccc initialise value of mx to be inherited at
ccc next grid time
mxnew = mx + 1
ccc initialise grid location to be processed
ix = 0
xnd = zero
ccc step through grid locations, from
ccc nearest to furthest from source in
ccc direction of release
200 ix = ix + 1
xnd = xnd + xresnd
ccc when first processing a grid location
ccc at which exhausts were not present
ccc at last grid time processed
if (ix.eq.mx+1) then
i = 0
ccc set iteration invariants
c1 = mx*xresnd - (Tnd**2)/2
c2 = Tnd - agend(mx)
c3 = cX/Tresnd - Tnd
ccc set initial estimate of delXnd for
ccc improvement by Newton's method
delXnd = Tnd*Tresnd - cX
ccc step through Newton's method
ccc iterations
250 i = i + 1
dum = c1 + delXnd
dum1 = c2/delXnd
del1 = delXnd/Tresnd + c3
ccc calc function whose root with
ccc respect to dimensionless
ccc location of exhaust-plume head
ccc is to be obtained
del = corat*(dum**2)*del1 - dum1
ccc calc derivative of this function
ccc with respect to dimensionless
ccc location of exhaust-plume head
del1 =
+ corat*dum*(dum/Tresnd +
+ 2*del1) + dum1/delXnd
if (del1.lt.sig) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
ccc when function is
ccc significantly valued
if (dabs(del).ge.sig) then
ccc calc improvement on root
delXnd = delXnd - del/del1
if (delXnd.lt.
+ Tnd*Tresnd-cX) then
print *,
+ 'alg fails: prog ends'
stop
endif
endif
if (i.lt.nitmax .and.
+ dabs(del).ge.sig) goto 250
if (i.gt.nitmax .or.
+ (i.eq.nitmax .and.
+ dabs(del).ge.sig)) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
ccc calc value of mx to be inherited at
ccc next grid time
mxnew = mx + idint(delXnd/xresnd)
if (mxnew.gt.nx) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
dum = (delXnd + cX -
+ Tnd*Tresnd)/cokF
ccc update mean height of exhaust-plume
ccc head resulting from buoyant rise
Zhnd = Zhnd + dum*MBnd(1)
concf = dum/Tresnd
if (iTgs.eq.mTgs .and.
+ iTlas.ge.iTrot) dagxnd = dum1
endif
ccc when exhausts are present at least
ccc one grid location further from source
ccc than at last grid time processed,
ccc solve at associated grid points by
ccc linear extrapolation
if (ix.gt.mx .and. mxnew.gt.mx) then
dum = (ix - mx)*xresnd/delXnd
agend(ix) = agend(mx) + dum*c2
znd(ix) = znd(mx) +
+ dum*(Zhnd - znd(mx))
ccc when processing grid locations at
ccc which exhausts were present at last
ccc grid time processed
else if (ix.le.mx) then
i = 0
ccc set initial estimate of
ccc dimensionless exhaust age for
ccc improvement by Newton's method
agennd = agend(ix-1)
ccc step through Newton's method
ccc iterations
300 i = i + 1
dum = (agennd -
+ agend(ix-1))/xresnd
dum1 = one -
+ (agennd - agend(ix))/
+ Tresnd - dum*Tnd
del = xnd -
+ (Tnd - agennd/2)*agennd
ccc calc derivative of function
ccc whose root with respect to
ccc dimensionless exhaust age is
ccc sought
del1 = -corat*del*
+ (del*(Tnd/xresnd +
+ 1/Tresnd) +
+ 2*dum1*
+ (Tnd - agennd)) -
+ 2*dum/xresnd
ccc calc function
dum = corat*dum1*del**2 - dum**2
if (dum1.lt.sig .or.
+ del.lt.sig .or.
+ del1.gt.-sig) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
ccc when function is significantly
ccc valued, calc improvement on root
if (dabs(dum).ge.sig)
+ agennd = agennd - dum/del1
if (i.lt.nitmax .and.
+ dabs(dum).ge.sig) goto 300
if (i.gt.nitmax .or.
+ (i.eq.nitmax .and.
+ dabs(dum).ge.sig) .or.
+ agennd.lt.
+ agend(ix-1)-sig1 .or.
+ agennd.gt.agend(ix)+sig1)
+ then
print *,
+ 'algorithm fails: prog ends'
print *,'try increasing resn'
stop
endif
dagTnd = (agennd - agend(ix))/Tresnd
agend(ix) = agennd
dagxnd = (agend(ix) -
+ agend(ix-1))/xresnd
ccc establish appropriate value of MBnd
ccc by linear extrapolation, for
ccc calculation of mean height of
ccc exhausts resulting from buoyant rise
del = agend(ix)/Tresnd
i = idint(del)
j = iTgs - i + 1
if (j.lt.1 .or. j.gt.iTgs+1) then
print *,'alg fails: prog ends'
stop
endif
if (j.gt.iTgs) then
dum = zero
else
dum = MBnd(j)
endif
if (j.eq.1) then
dum1 = MBnd(1)
else
dum1 = MBnd(j-1)
endif
del1 = dum +
+ (del - dble(i))*(dum1 - dum)
dum = one - dagTnd
dum1 = dum/xresnd
dum = (dum - dagxnd*Tnd)/cokF
del = dagxnd/Tresnd
c1 = znd(ix)*del + znd(ix-1)*dum1
del = del + dum1
if (del.lt.sig) then
print *,
+ 'algorithm fails: prog ends'
stop
endif
dum1 = (dum*del1 + c1)/del
if ((invMO.le.zero .and.
+ (dum1.lt.znd(ix-1)-sig1 .or.
+ dum1.gt.znd(ix)+sig1)) .or.
+ (invMO.gt.zero .and.
+ dum1.lt.-sig1)) then
print *,
+ 'algorithm fails: prog ends'
print *,'try increasing resn'
stop
endif
znd(ix) = dum1
concf = dum
endif
ccc when grid time furthest into ground
ccc run currently for processing is
ccc reached, calc dimensional exhaust
ccc characteristics as function exhaust
ccc age, with ages as separated by 1 s
ccc intervals: this is for purposes of
ccc output, or for initialisation of
ccc dispersion regime after upward
ccc rotation of aircraft is begun
if (iTgs.eq.mTgs) then
ccc initialise youngest appropriate age,
ccc in integral s
if (ix.eq.1) magnex = 1
ccc when considering a grid location
if (ix.le.mxnew) then
ccc calc oldest appropriate age, in
ccc integral s
iage = idint(agend(ix)*ts)
if (iage.lt.0 .or.
+ iage.gt.iTgs/ngsps .or.
+ iage.gt.iTrot) then
print *,
+ 'alg fails: prog ends'
stop
endif
do 350 i=magnex,iage
j = i - (nage - iTfir)
if (j.gt.0) xLol(j) = xL(j)
dum = (xnd -
+ (Tnd - agend(ix)/2)*
+ agend(ix))*xs
xL(i) = dum
area(i) = (coke*dum)**2
z(i) = znd(ix)*zs
dum = concf*concs
if (ix.le.mx)
+ dum = dum/dagxnd
conc(i) = dum
ccc when upward rotation of
ccc aircraft has been begun,
ccc calc QB, dQBT and QBF so
ccc dispersion may be
ccc evaluated later
if (iTlas.ge.iTrot) then
j = mTgs - (i - 1)*ngsps
if (j.lt.1 .or.
+ j.gt.mTgs) then
print *,'alg fails:'
print *,'prog ends'
stop
endif
QB(i)=
+ MBnd(j)*dagxnd*QBs
dQBT(i)=
+ dMBTnd(j)*dagxnd*dQBTs
QLF(i) =
+ cokF*F0*area(i)*conc(i)
endif
350 continue
magnex = iage + 1
endif
ccc when there are no more grid
ccc locations to be processed,
ccc draw on characteristics of
ccc exhaust-plume head
if (ix.ge.mxnew) then
ccc calc oldest appropriate age, in
ccc integral s
iage = mTgs/ngsps
do 360 i=magnex,iage
j = i - (nage - iTfir)
if (j.gt.0) xLol(j) = xL(j)
dum = (mx*xresnd + delXnd -
+ (Tnd**2)/2)*xs
xL(i) = dum
area(i) = (coke*dum)**2
z(i) = Zhnd*zs
conc(i) = concf*concs
ccc when upward rotation of
ccc aircraft begins,
ccc calc QB, dQBT and QBF so
ccc dispersion may be
ccc evaluated later
if (iTlas.ge.iTrot) then
QB(i) =
+ MBnd(ngsps)*dagxnd*QBs
dQBT(i) =
+ dMBTnd(ngsps)*dagxnd*
+ dQBTs
QLF(i) =
+ cokF*F0*area(i)*conc(i)
endif
360 continue
endif
endif
ccc when grid locations remain to be
ccc processed, progress to next
if (ix.lt.mxnew) goto 200
500 continue
ccc when upward rotation of aircraft begins,
ccc prepare to evaluate dispersion of exhausts
ccc released earlier
if (iTlas.gt.iTrot) then
iTold = iTrot
goto 120
endif
ccc when upward rotation of aircraft has been begun
else
ccc calc number of grid times to be processed
mTgs = (iTlas - iTold)*ng2sps
ccc when rotation had been begun by last grid
ccc time at which outputs were generated, store
ccc val of xL at grid time of last output
if (iTfir.ge.iTrot) then
do 600 iage=1,iTrot
xLol(iage) = xL(iage)
600 continue
endif
ccc step through grid times
do 1000 iTgs=1,mTgs
ccc step through gridded ages, as separated
ccc by 1 s intervals
do 900 iage=1,iTrot
if (area(iage).lt.sig .or.
+ QLF(iage).lt.sig) then
print *,'algorithm fails: prog ends'
stop
endif
dum = z(iage)*invMO
ccc when mean height exhausts have risen
ccc to buoyantly lies in lower ambient
ccc thermal layer
if (dabs(dum).lt.one) then
del = karho2
ccc when mean height exhausts have risen
ccc to buoyantly lies in upper ambient
ccc thermal layer (in stable
ccc stratification)
else if (dum.ge.one) then
del = karhu2
else
del = zero
endif
ccc when stable ambient stratification
ccc is absent
if (invMO.le.zero) then
if (del.lt.sig) then
dum = one
else
dum = dsinh(del)/del
endif
dum1 = dcosh(del)
k = 1
else
if (del.lt.sig) then
dum = one
else
dum = dsin(del)/del
endif
dum1 = dcos(del)
k =-1
endif
ccc calc dimensional exhaust characteristics
del1 = dum*dQBT(iage)/ng2sps +
+ dum1*QB(iage)
dQBT(iage) =
+ dum1*dQBT(iage) +
+ k*ng2sps*(del**2)*dum*QB(iage)
c1 = QB(iage)
QB(iage) = del1
del = 2*ng2sps*rho
del1 = QLF(iage)/del
dum = 3*coke*del1/area(iage)**cthrhf
dum1 =
+ dsqrt(one + (c1/QLF(iage))**2) +
+ dsqrt(one + (QB(iage)/QLF(iage))**2)
dum = (one + dum*dum1)**twothd
c2 = area(iage)
area(iage) = dum*area(iage)
conc(iage) = conc(iage)/dum
del = 2*del
xL(iage) = xL(iage) +
+ del1*(1/c2 + 1/area(iage))
z(iage) =
+ z(iage) + (c1/c2 +
+ QB(iage)/area(iage))/del
900 continue
1000 continue
endif
write (50,*) 'old time since start of ground run (s)'
write (50,*) iTfir
write (50,*) 'new time since start of ground run (s)'
write (50,*) iTlas
write (50,*) 'number ages tracked'
ccc initialise number of gridded exhaust ages at which
ccc dispersive characteristics are available
ntr = 0
ntrold = 0
ccc set up output-loop invariants
c1 = coke/rho
ccc set nage to number of gridded ages
if (iTlas.gt.iTrot) then
nage = iTrot
else
nage = iTlas
c1 = c1*cokF*F0
endif
ccc when old time preceded upward rotation of aircraft,
ccc update itrac and include exhausts emitted after old
ccc time
if (iTfir.lt.iTrot) then
j = nage - iTfir
do 1050 i=nage,1,-1
if (i.gt.j) then
itrac(i) = itrac(i-j)
else
itrac(i) = 1
endif
1050 continue
endif
ccc set c2 to velocity scale of free convection (in m s-1)
dum = -iblht*invMO
if (dum.gt.zero) c2 = therc8*(dum**third)*fricv
ccc step through set of exhaust ages three times, first
ccc to obtain nunber of gridded ages at which dispersive
ccc characteristics are available, second to output
ccc these characteristics, third to output reason where
ccc characteristics are not available
do 1500 igo=1,3
if (igo.eq.2) then
write (50,*) ntr
if (ntr.gt.0) write (50,'(a57)')
+ ' t, Xs,delta_xL, xL, zb, zt,delta_y, conc'
else if (igo.eq.3 .and. ntr.lt.ntrold) then
write (50,*)
+ 'ages ceased tracking since old time'
write (50,'(a36)')
+ ' t, Xs, reason'
endif
ccc step through gridded ages
do 1100 iage=1,nage
ccc when dispersion model was applicable at old
ccc time
if (itrac(iage).eq.1) then
if (igo.eq.1) ntrold = ntrold + 1
if (area(iage).lt.sig) then
print *,'algorithm fails: prog ends'
stop
endif
del = dsqrt(area(iage))/2
ccc when plume is in contact with ground
if (z(iage).lt.del) then
del1 = 2*z(iage)
else
del1 = z(iage) + del
endif
if (del1.lt.sig) then
print *,'algorithm fails: prog ends'
stop
endif
ccc when exhausts have risen buoyantly to a
ccc mean height greater than magnitude of
ccc Monin-Obukhov length in ambient
ccc convection, set velocity scale of ambient
ccc diffusion to that of free convection
if (invMO*z(iage).lt.-one) then
dum = c2
ccc set velocity scale of ambient diffusion
ccc to friction velocity
else
dum = fricv
endif
if (iTlas.le.iTrot) then
dum1 = c1*conc(iage)
else
dum1 = c1*dsqrt(QB(iage)**2 +
+ QLF(iage)**2)/area(iage)
endif
ccc when exhaust dispersion reflects trapping
ccc at inversion capping boundary layer or
ccc action of ambient turbulence
if (del1.gt.dble(iblht) .or.
+ dum1.lt.dum) then
if (igo.eq.3) then
itrac(iage) = 0
i = iage + iTlas - nage
j = idnint(acc*((iTlas - i)**2)/2)
if (del1.gt.dble(iblht) .and.
+ dum1.ge.dum) then
write (50,'(i4,i5,a27)')i,j,
+ 'base of inversion reached'
else if (del1.le.
+ dble(iblht) .and.
+ dum1.lt.dum) then
write (50,'(i4,i5,a27)') i,j,
+ 'ambient diffusion underway'
else
write (50,'(i4,i5,a55)') i,j,
+ 'ambient diffusion underway / base of inversion reached'
endif
endif
else
if (igo.eq.1) then
ntr = ntr + 1
else if (igo.eq.2) then
i = iage + iTlas - nage
j = idnint(acc*((iTlas - i)**2)/2)
if (iTfir.ge.iTrot) then
k = 0
else
k = nage - iTfir
endif
if (iage.gt.k) then
dum = xL(iage) - xLol(iage-k)
else
dum = xL(iage)
endif
ccc when plume is in contact with
ccc ground
if (z(iage).lt.del) then
dum1 = zero
del = area(iage)/del1
else
dum1 = z(iage) - del
del = 2*del
endif
write
+ (50,'(i4,i5,f9.1,3f7.1,f8.1,e10.3)')
+ i,j,dum,xL(iage),dum1,del1,del,
+ conc(iage)
endif
endif
endif
1100 continue
1500 continue
if (ntrold.lt.1 .or. ntrold.lt.ntr) then
print *,'algorithm fails: prog ends'
stop
endif
if (ntr.gt.0) goto 100
endif
close (50)
end