-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcpd_hr_swelling.f
1733 lines (1662 loc) · 58.3 KB
/
cpd_hr_swelling.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C THIS IS THE CPDCP-NLG MODEL
C
C THIS MODEL WAS DEVELOPED BY SANDIA NATIONAL LABORATORIES UNDER
C FWP 0709 FOR THE DEPARTMENT OF ENERGY'S PITTSBURGH ENERGY
C TECHNOLOGY CENTER AND THE DOE DIVISION OF ENGINEERING AND GEOSCIENCES
C THROUGH THE OFFICE OF BASIC ENERGY SCIENCES;
C AND BY THE UNIVERSITY OF UTAH THROUGH FUNDING FROM
C THE ADVANCED COMBUSTION ENGINEERING RESEARCH CENTER (ACERC), WHICH
C IS PRINCIPALLY SPONSORED BY THE NATIONAL SCIENCE FOUNDATION, THE
C STATE OF UTAH, AND BY A CONSORTIUM OF INDUSTRIAL COMPANIES.
C THE CODE WILL NOT BE FORMALLY LICENSED. NEITHER THE U.S. OR THE
C DOE, NOR ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED,
C OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
C COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT, OR
C PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD INFRINGE PRIVATELY
C OWNED RIGHTS.
C
C The CPD model is intended to solve pyrolysis rate equations
C based on a percolative bond-breaking scheme. This version includes the
C flash distillation program to distinguish between tar and metaplast.
C This program also includes a crosslinking scheme. (January, 1991)
C
C
C Most recent modifications to this model include (a) the nitrogen release
C model of Perry, (b) the model of Genetti to break the light gas into
C species based on a correlation, and (C) slight modification to mdel
C to account for the mass associated with c0.
C These modifications were made at BYU by Dominic Genetti in his
C M.S. Thesis work (1999) and Steve Perry in his Ph.D. work (1999).
C
C
program CPDCP
C
C This version is coupled with a solver for the particle energy and
C momentum equations.
C UNITS ARE G,K,CM,S,CAL
C PARAMETER INITIALIZATION
C
C A BLOWING CORRECTION TO THE HEAT TRANSFER MODEL IS EMPLOYED.
C
C Merrick heat capacity correlations are used
C This version uses the heat capacity at the initial temperature as recommended
C by Maloney, D. J., R. Sampath and J. W. Zondlo, Combustion and Flame, 116, 94-104 (1999).
C
C Water is included through Omegaw
C
IMPLICIT real*4 (A-H,O-Z)
save
C Input parameters are found in PERKIN
C Outut parameters are found in PERKOT
DIMENSION XT(50),TGC(50),zv(50),vpz(50)
CHARACTER*80 A,vname,tgasin,output,infile,outgas,outnit,hrfile
REAL KG,NU,KMJ1,KMJ2
character*80 perkin,perkot,perkot2,perkot3
character*1 ans
dimension y(4),yp(4),ypp(4),ypred(4),tim(50),tem(50),gasf(9,10)
dimension yygas(5),ffgas(5)
real*4 l0,l,ma,kp,ft(35),mt(35),mgas,mtot
real*4 mwchar,machar,mwcharold,mw1,mdel,mb,trateold,volc
real*4 sp,smin,sigmdel,cscale,svar,cHR,cvisc,cext,sratio,sc
integer*2 lib
C intar = .true. calculates tar molecular weight distribution
C in subroutine perkp.
logical intar,idiff,imolw,ipred,inside
C ipred = .true. when on the predictor step
common/init/l0,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar
common/rate/ab,eb0,ebsig,ac,ec0,ecsig,ag,eg0,egsig,rg
common/tarmw/press
common/timer/tms,time
common/spheat/yelem(5)
common/nitmod/fnca0,fnca,an,en0,ensig
common/lgmod/yf,xoc,yhc
data xt/50*0./,tgc/50*0./,zv/50*0./,vpz/50*0./
data tim/50*0./,tem/50*0./,ft/35*0./,mt/35*0./
data y/4*0./,yp/4*0./,ypp/4*0./,ypred/4*0./
data nmax/10/,nt/1/,intar/.false./,idiff/.false./
data ftar0,fgas0,fchar0/0,0,1./,ip/30/,zero/0.0/
data small/1.e-7/ftarold/0.0/
C y1 = l labile bridges
C y2 = del ends
C y3 = C char links
C y4 = fnca mass fraction of nitrogen in site
data g0,g1,g2/0.,0.,0./
DATA PI/3.14159/,DELHV/0./,PR/0.7/,delhw/-540./,cpw/1.0/
data TMS,Xmm,V,FRAC,BLOW/4*0.0,1./
xwb=0.0
rtot=0.0
tratem=0.0
trateold=0.0
trcheck=0
volc=0.0
tpmax=0.0
C
C
IX = 1
IV = 1
X = 0.0
C
C input parameters
C
print*,' Enter main input file'
C read(*,1)infile
read(*,*)infile
print*,infile
open(unit=21,name=infile,type='unknown')
print*,' Enter input particle velocity file'
C read (21,1)vname
read (21,*)vname
open(unit=25,name=vname,type='unknown')
read(25,*)nv
C do 710 i=1,20
C read(25,1)a
C write(*,1)a
C if(a(1:1).eq.' ')go to 711
C710 continue
C711 CONTINUE
do 1001 i=1,nv
read(25,*)zv(i),vpz(i)
C write(*,*)zv(i),vpz(i)
zv(i) = zv(i)/10.
write(*,*)zv(i),vpz(i)
1001 continue
C print*,' Error: Too many velocity data points'
C1112 continue
C nv = i-1
C Enter input gas temperature file
read (21,*)tgasin
C read (21,1)tgasin
open(unit=22,name=tgasin,type='unknown')
C read in gas temperature profile
read(22,*)nx
C do 10 i=1,20
C read(22,1)a
C1 format(a)
C if(a(1:1).eq.' ')go to 11
C10 continue
C11 CONTINUE
do 1111 i=1,nx
read(22,*)xt(i),tgc(i)
xt(i) = xt(i)/10.
write(*,*)xt(i),tgc(i)
1111 continue
C print*,'too many data points for gas temperature'
C stop
C1002 continue
C nx = i-1
print*,' Enter name for output file'
read (21,*)output
write(*,*)output
if(output(1:1).eq.' ')output(1:11)='ubcp.out'
open(unit=1,name=output,type='unknown')
print*,' Enter name for nitrogen release output file'
read (21,*)outnit
write(*,*)outnit
open(unit=23,name=outnit,type='unknown')
read (21,*)outgas
write(*,*)outgas
open(unit=27,name=outgas,type='unknown')
C
read(21,*)TIMAX
C read(21,*)TG0
C read(21,*)VG0 !cm/s
TG0=tgc(1)
VG0=vpz(1)
write(*,*)Tg0,vg0
read(21,*)RHOP !G/CM**3
read(21,*)DP !CM
C read(21,*)Tau0 !adjustable to match particle temperatures (or mass release)
C Rhop=rhop*tau0 !this is one way to implement the correction
read(21,*)Swell !fraction final diameter increase
swell=swell-1.0
if((swell+1.0).gt.0.5) then
sratio=swell+1.0 !use experimental swelling ratio
else
swell=0.0
sratio=0.0 !triggers use of swelling correlation by Randy Shurtz
endif
sc=sratio
read(21,*)DELHV !CAL/G (- MEANS ENDOTHERMIC)
read(21,*)OMEGAW
read(21,*)OMEGAA
read(21,*)EMIS
read(21,*)TWALL
read(21,*)dt,dtmax,iprint
C coal input parameters
C read in measured p0 from NMR data
read(21,*)p0
print*,'p0',p0
C estimate c0
read(21,*)c0
C read in sigma+1 (coordination number) from NMR data
read(21,*)sigp1
C read in real MW/Cluster from NMR data (includes side chains)
read(21,*)mw1
C read in real mdel from NMR data
read(21,*)mdel
print*,'p0,c0,sigp1,mw1,mdel',p0,c0,sigp1,mw1,mdel
C kinetic parameters
read(21,*)ab
read(21,*)eb0
read(21,*)ebsig
read(21,*)ac
read(21,*)ec0
read(21,*)ag
read(21,*)eg0
read(21,*)egsig
read(21,*)acr
read(21,*)ecr
C kinetic parameters for light gas nitrogen release from char:
C A. release by reaction with free radicals (fast)
read(21,*)arad
read(21,*)erad
read(21,*)fstable
C B. release by high T decomposition (slow)
read(21,*)an
read(21,*)en0
read(21,*)ensig
C pressure in atmospheres
read(21,*)press
C read in DAF composition of the coal (CHNOS)
do 12 i=1,5
read(21,*)yelem(i)
12 continue
read(21,*)vastm
read(21,*)cpswitch
read(21,*)hrfile
C READ(21,*)THTR !These temperatures are for a Sandia drop-tube (CDL);
C READ(21,*)TTUB !setting them to the the wall temperature makes the view factors irrelavent
THTR=TWALL
TTUB=TWALL
open(unit=29,name=hrfile,type='unknown')
write(29,*) "Pressure (atm) ",press
write(29,*) "Particle Diameter (microns) ",(DP*1.e4)
write(29,*) "Moisture Fraction ", OMEGAW
write(29,*) "Dry Ash Fraction ",OMEGAA
write(29,*) "CHNOS daf Percentages ",
& yelem(1),yelem(2),yelem(3),yelem(4),yelem(5)
write(29,*)"ASTM daf Volatiles Percent",vastm
C Predict NMR parameters from Genetti's correlation if not provided, convert composition from percentages to fractions
if((p0.lt.0.0).or.(mdel.lt.0.0).or.
& (mw1.lt.0.0).or.(sigp1.lt.0.0).or.(c0.lt.0.0)) then
mdel=(421.957)+(-8.64692)*yelem(1)
& +(0.0463894)*yelem(1)*yelem(1)+(-8.47272)*yelem(2)
& +(1.18173)*yelem(2)*yelem(2)+(1.15366)*yelem(4)
& +(-0.0434024)*yelem(4)*yelem(4)
& +(0.556772)*vastm+(-0.00654575)*vastm*vastm
mw1=(1301.41)+(16.3879)*yelem(1)+(-0.187493)*yelem(1)*yelem(1)
& +(-454.773)*yelem(2)+(51.7109)*yelem(2)*yelem(2)
& +(-10.072)*yelem(4)+(0.0760827)*yelem(4)*yelem(4)
& +(1.36022)*vastm+(-0.0313561)*vastm*vastm
p0=(0.489809)+(-0.00981566)*yelem(1)
& +(0.000133046)*yelem(1)*yelem(1)+(0.155483)*yelem(2)
& +(-0.0243873)*yelem(2)*yelem(2)+(0.00705248)*yelem(4)
& +(0.000219163)*yelem(4)*yelem(4)+(-0.0110498)*vastm
& +(0.000100939)*vastm*vastm
sigp1=(-52.1054)+(1.63872)*yelem(1)
& +(-0.0107548)*yelem(1)*yelem(1)+(-1.23688)*yelem(2)
& +(0.0931937)*yelem(2)*yelem(2)+(-0.165673)*yelem(4)
& +(0.00409556)*yelem(4)*yelem(4)+(0.00926097)*vastm
& +(-8.26717E-05)*vastm*vastm
if (yelem(1).gt.85.9) then
c0=min((0.1183*yelem(1)-10.16),0.36)
elseif (yelem(4).gt.12.5) then
c0=min((0.014*yelem(4)-0.175),0.15)
else
c0=0.0
endif
endif
C calculate terms in swelling model as developed by Randy Shurtz
write(29,*)"c0 ",c0
write(29,*)"p0 ",p0
write(29,*)"(sigma+1) ",sigp1
write(29,*)"MW ",mw1
write(29,*)"Mdel ",mdel
smin=((1.-0.01*vastm)*(1.-omegaa)+omegaa)**(0.3333333333)
sigmdel=sigp1/mdel
if((sigmdel.le.0.01831).or.(sigmdel.ge.0.30084)) then
svar=0.0
elseif (sigmdel.le.0.2067) then
svar=1.6852*sigmdel-0.030856
else
svar=-3.3722*sigmdel+1.0145
endif
write(29,*)"(sigma+1)/mdel ",sigmdel
write(29,*)"smin ",smin
write(29,*)"svar ",svar
if ((sigmdel.le.0.1062).or.(sigmdel.ge.0.2541)) then
cHR=0.0
else
cHR=-191.293*sigmdel*sigmdel+68.921*sigmdel-5.161
endif
write(29,*)"cHR ",cHR
cvisc=1.0
cext=1.0
cscale=0.0
if(press.le.1.0) then
sp=1.0
else
cvisc=7.77
cext=3.47
cscale=38.89*sigp1-167.10
if(sigp1.lt.4.297) then
cscale = 0.0
endif
sp=1.0+cscale*((log(press))**cvisc)/(press**cext)
endif
write(29,*)"cvisc ",cvisc
write(29,*)"cext ",cext
write(29,*)"cscale ",cscale
write(29,*)"sp (swelling pressure factor) ",sp
yelem(1)=yelem(1)*.01
yelem(2)=yelem(2)*.01
yelem(3)=yelem(3)*.01
yelem(4)=yelem(4)*.01
yelem(5)=yelem(5)*.01
print*,'p0,c0,sigp1,mw1,mdel',p0,c0,sigp1,mw1,mdel
print*,'C,H,N,O,S fracs',yelem(1),yelem(2),yelem(3),yelem(4),
& yelem(5)
C save initial char NMR parameters as coal NMR parameters
C (char parameters calculated independent of those using
C empirical correlation for mdel)
mwchar = mw1
mwcharold = mwchar
machar = mwchar-sigp1*mdel
C adjust mdel to correct for c0 (Steve Perry, May 99)
mdel = mdel/(1.0-c0)
C empirical correlation to allow a small portion of alpha-carbon to stay with
C the aromatic cluster
mdel = mdel-7
C now calculate other chemical structure coefficients
l0 = p0 - c0
mb = 2.*mdel
ma = mw1-sigp1*mdel
sma = mw1-sigp1*(mdel+7)
sig = sigp1-1
finf = 1./(2.*ma/(mb*sigp1*(1-c0))+1)
rba = mb/ma
print*,'rba=',rba
fnit = yelem(3)
fnt = fnit
fnca0 = fnit*mw1/machar
Grav = 980.
dp0 = dp
C initialize variables
y(1) = l0
y(2) = 2.*(1.-c0-l0)
y(3) = c0
y(4) = fnca0
aind0 = l0 + (1.-c0-l0)
siginv = 1./sig
pstar = 0.5*siginv
ynhcn = 0.0
yntar = 0.0
yytar = 0.0
yf = 0.0
inside=.true.
C START calculation DO LOOP
Rg = 1.987 !CAL/GMOLE K
iii = 0
fcross = 0.0
C CALCULATE INITIAL PARTICLE VELOCITY
CALL PROPS(TG0,TG0,RHOG,UG,KG,CPG,diffw,press)
C VP = VG0 - DP**2*RHOP*Grav/(18*UG)
VP = Vg0
VG = VG0
C START DO LOOP
V = 0.
Tp = tg0
Tg = tg0
MOVE = 1
TIME = 0.0
C for now, assume that the apparent density is indicative of the as
C received coal.
ALPHAP = (4./3.)*PI*(DP/2.)**3*RHOP
alfa = alphap*omegaa
alfw = alphap*omegaw
ALFC = ALPHAP*(1-OMEGAA-Omegaw)
alfc0 = alfc
ALFV = 0.
ALPHA = Alphap
FRAC = ALFC/ALPHAP
W = ALFW/ALPHAP*100.
AP = PI*DP**2
RE = RHOG*ABS(VG-VP)*DP/UG
NU = 2. + 0.6*RE**.5*PR**.333
H = NU*KG/DP
SIGMA = 1.335E-12 !CAL/S CM**2 K**4
HEAT = 0.0
v = 0.0
C calculate O/C and H/C ratios for light gas model
xoc = (yelem(4)/16)/(yelem(1)/12)
yhc = yelem(2)/(yelem(1)/12)
C
R = 1.987 !CAL/GMOLE K
dpm = dp*1.e4
write(1,5000)dpm
5000 format('C Coal Calculations'/
x'C Initial Particle Diameter = ',f6.2,' microns'/'C')
write(6,5010)
5010 format(' Tim (ms) Dis (mm) Tp (K) Tg (K)',
x ' Vol. Char Blow')
write(1,202)
202 format('C time(ms) dis(mm) Tp(K) Tg(K) % Vol.',
x' fchar fcross ftar fmet Trate')
write(23,204)
204 format('C time(ms) temp(K) MW Nsite Nchar fnchar'
x,' fntar fnhcn',
x' fntot')
write(27,205)
205 format('C time(ms) fgas fh20 fco2 fch4',
x' fco fother yh20 yco2 ych4',
x' yco yother Xgas')
call heatcp(tp,cpc)
IF (CPSWITCH.GT.0.0) THEN
CPC=CPC*CPSWITCH
ENDIF
cpcold=cpc
call heatap(tp,cpa)
fntar = 0.0
C CALL HEATWP(tp,cpw)
cp = (alfc*cpc + alfa*cpa + alfw*cpw)/(alpha)
WRITE(6,1000)TMS,Xmm,Tp,TG0,v,fchar,BLOW
write(1,1000)TMS,Xmm,Tp,TG,v,fchar,fcross,
x ftar,fmet,trate
DO 100 I=1,100000 !goes to continue statement 100, exits to statement 3333
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PREDICTOR
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
fvolold = fvol
fcrossold = fcross
fmetold = fmet
XP = X + VP*DT
C CALCULATE GAS VELOCITY FROM Thermocouple Measurements
IF(XP.LE.XT(nx))THEN
TG = (XT(IX+1)-XP)/(XT(IX+1)-XT(IX))*(TGC(IX)-TGC(IX+1))
X + TGC(IX+1)
ELSE
PRINT*,'REACHED END OF GAS TEMPERATURE CORRELATION'
if(inside.eqv..false.)then
print*,'!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!'
print*,'O/C and H/C ratios are outside the '
print*,'bounds of the library coals. '
print*,'Estimation of Light gas distribution '
print*,'is based on library coal no.',lib
endif
go to 3333
ENDIF
if(tg.gt.4000.)then
print*,' >>>> WARNING <<<<<'
print*,' gas temperature too high----',tg
endif
C interpolate to get particle velocity
IF(XP.LE.zv(nv))THEN
vpp = (zv(iv+1)-xP)/(zv(iv+1)-zv(iv))*(vpz(iv)-vpz(iv+1))
X + vpz(iv+1)
ELSE
PRINT*,'REACHED END OF PARTICLE VELOCITY CORRELATION'
go to 3333
ENDIF
C CALCULATE GAS PROPERTIES
CALL PROPS(TG,Tp,RHOG,UG,KG,CPG,diffw,press)
C INTEGRATE PARTICLE VELOCITY
C Use this section if gas velocity .NE. particle velocity
C VLAG = (VP-VG)
C RE = RHOG*DP*ABS(VLAG)
C RE = RE/UG
C SUMP = -18.*UG*VLAG*(1.+.15*RE**.687)/(DP**2*RHOP)+Grav
C SUMP = -18.*UG*VLAG/(DP**2*RHOP)+Grav
C VPP = VP + SUMP*DT
RE = 0.
C ENERGY EQUATION
C
C-- CONVECTION
NU = 2. + .6*RE**.5*PR**.333
B = CPG*(rtot)/(2.*PI*DP*KG)
IF(B.GE.1.E-4)THEN
BLOW = B/(EXP(B)-1)
ELSE
BLOW = 1.0
ENDIF
H = BLOW*NU*KG/DP
QCONV = H*AP*(TG-Tp)
C-- mass transfer
if(alfw.gt.0.)then
Bw = (rtot)/(2.*PI*DP*diffw*rhog)
IF(Bw.GE.1.E-4)THEN
BLOwW = Bw/(EXP(Bw)-1)
ELSE
BLOwW = 1.0
ENDIF
else
bloww = 1.0
endif
C-- RADIATION
Z = XP
RAD = 2.54 ! 1" radius burner or injector probe?
if(z.gt.0.)then
R2 = RAD/Z
R3 = .125*2.54/Z !1/8" radius feeder tube?
else
r2 = 1.e12
r3 = 1.e12
endif
f12 = .5*(1.-1./sqrt(1.+r2**2)) !view factor to burner/drop tube injector
f13 = .5*(1.-1./sqrt(1.+r3**2)) !view factor to some other tube that is the same distance; feeder tube?
FAC = AP*SIGMA*EMIS
QTOP = (f12-F13)*(THTR**4-Tp**4)*FAC
QTUBE = F13*(TTUB**4-Tp**4)*FAC
QWALL = (1.-F12)*(TWALL**4-Tp**4)*FAC
QRAD = QWALL + QTOP + QTUBE
C-- Water Evaporation Rate
if(alfw.gt.0.)then
xw0 = exp(18.3036-3816.44/(Tp-46.13))/760.
xw0 = min(xw0/press,1.) !It should be the total pressure to get the mole fraction, not just 1 atm
xw0 = max(xw0,0.)
rwp = bloww*2.*rhog*diffw*pi*dp*(xw0-xwb)/(1.-xw0)
else
rwp = 0.
endif
C-- Coal pyrolysis rate
call perks(y,ypp,Tp)
C free radical light gas nitrogen release mechanism
if((mw1-mwchar)/mw1.GT.fstable)ypp(4)=ypp(4)-y(4)*arad*
x exp(-erad/rg/Tp)*(mwcharold-mwchar)/mwchar*
x machar/mwchar/dt
C COMPONENT MASS CONSERVATION
do 50 j=1,4
ypred(j) = y(j) + dt*ypp(j)
ypred(j) = max(ypred(j),zero)
50 continue
C crosslinking rate
fracr = 1.
if(fmetold.gt.small.and.acr.gt.0.)then
ratecr = acr*exp(-ecr/rg/tp)*fmetold*dt
fracr = 1.-ratecr/fmetold
fmet = fmetold-ratecr
fcross = fcrossold+ratecr
if(fmet.lt.0.)then
fcross = fcrossold-ratecr+fmet+ratecr
fmet = 0.0
fracr = 0.
endif
endif
C get product distribution from ypred array
if(ypred(1).gt.small)intar = .true.
call perkp(ypred,mgas,ftar,ftart,fgas,fchart,ft,mt,intar)
intar = .false.
gasmw = rba*ma/2.
C flash distillation
if(fgas.ge.1.e-5)then
imolw = .false.
if(mod(iii,iprint).eq.0)imolw=.true.
ipred = .true.
call flash(fgas,gasmw,ft,mt,fracr,ftar,fmet,
x tp,press,nmax,imolw,ipred)
elseif(fgas.lt.1.e-5)then
fmet = ftart
ftar = 0.
endif
intar = .false.
fvol = fgas+ftar
fchar = 1.-fvol
dvdt = (fvol-fvolold)/dt*alfc0
C done with pyrolysis rate!
rtotp = dvdt + rwp
delhw= 6.7845E-08*(tp**4) - 1.1123E-04*(tp**3)
x + 6.8595E-02*(tp**2) - 1.8093E+01*tp + 1.1229E+03
HEAT = dvdt*delhv + rwp*delhw
C calculate heat capacity
fmetmaxp=max(fmetmaxp,fmet)
if (CPSWITCH.ge.0.0) then
C call heatcp(tp,cpc)
C cpc=(fmet/fmetmaxp)*cpcold+(1.0-fmet/fmetmaxp)*cpc
cpc=cpcold
else
CALL HEATCP(TP,CPC)
endif
call heatap(tp,cpa)
cp = (alfc*cpc + alfa*cpa + alfw*cpw)/(alpha)
C
TRATEP = (QCONV+QRAD+HEAT)/(ALPHA*CP)
TPRED = Tp + TRATEP*DT
C print*,'<trate> qconv,qrad,heat,alpha,cp,tratep,tp,tg',
C x qconv,qrad,heat,alpha,cp,tratep,tp,tg
C COMPONENT MASS CONSERVATION
alfwp = alfw - rwp*dt
ALFVP = fvol*alfc0
ALFcP = fchar*alfc0
ALFCP = MAX(ALFCP,0.0)
alfwp = max(alfwp,0.)
ALPHA = (ALFCP + alfa + alfwp)
omegaa = alfa/alpha
C particle swelling
l = ypred(1)
dp = dp0*(1.+swell*(1.-l/l0))
AP = PI*DP**2
C PARTICLE DENSITY CHANGES DURING DEVOLATILIZATION
RHOP = ALPHA/((4./3.)*PI*(DP/2)**3)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C CORRECTOR
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
X = X + 0.5*(VPP+VP)*DT
C Interpolate to get gas temperature
IF(X.LE.XT(nx))THEN
TG = (XT(IX+1)-X)/(XT(IX+1)-XT(IX))*
X (TGC(IX)-TGC(IX+1)) + TGC(IX+1)
ELSE
PRINT*,' REACHED END OF GAS TEMPERATURE CORRELATION'
go to 3333
ENDIF
C interpolate to get particle velocity
IF(X.LE.zv(nv))THEN
vp = (zv(iv+1)-x)/(zv(iv+1)-zv(iv))*(vpz(iv)-vpz(iv+1))
X + vpz(iv+1)
ELSE
PRINT*,'REACHED END OF PARTICLE VELOCITY CORRELATION'
go to 3333
ENDIF
C CALCULATE GAS PROPERTIES
CALL PROPS(TG,TPRED,RHOG,UG,KG,CPG,diffw,press)
C INTEGRATE PARTICLE VELOCITY
C Use this section if gas velocity .NE. particle velocity
C VLAG = (VP-VG)
C RE = rhog*DP*ABS(VLAG)
C RE = RE/UG
C SUM = -18.*UG*VLAG*(1.+.15*RE**.687)/(DP**2*RHOP)+
C SUM = -18.*UG*VLAG/(DP**2*RHOP)+Grav
C VP = VP + 0.5*(SUMP+SUM)*DT
RE = 0
C ENERGY EQUATION
C
C-- CONVECTION
NU = 2. + 0.6*RE**.5*PR**.333
B = CPG*(rtotp)/(2.*PI*DP*KG)
IF(B.GE.1.E-4)THEN
BLOW = B/(EXP(B)-1)
ELSE
BLOW = 1.0
ENDIF
H = BLOW*NU*KG/DP
QCONV = H*AP*(TG-TPRED)
C-- mass transfer
if(alfw.gt.0.)then
Bw = (rtotp)/(2.*PI*DP*diffw*rhog)
IF(Bw.GE.1.E-4)THEN
BLOwW = Bw/(EXP(Bw)-1)
ELSE
BLOwW = 1.0
ENDIF
else
bloww = 1.
endif
C-- RADIATION
Z = X
RAD = 2.54
if(z.gt.0.)then
R2 = RAD/Z
R3 = .125*2.54/Z
else
r2 = 1.e12
r3 = 1.e12
endif
f12 = .5*(1.-1./sqrt(1.+r2**2))
f13 = .5*(1.-1./sqrt(1.+r3**2))
FAC = AP*SIGMA*EMIS
QTOP = (f12-F13)*(THTR**4-TPRED**4)*FAC
QTUBE = F13*(TTUB**4-TPRED**4)*FAC
QWALL = (1.-F12)*(TWALL**4-TPRED**4)*FAC
QRAD = QWALL + QTOP + QTUBE
C-- water evaporation rate
if(alfw.gt.0.)then
xw0 = exp(18.3036-3816.44/(Tpred-46.13))/760.
xw0 = min(xw0/press,1.)
xw0 = max(xw0,0.)
rw = bloww*2.*rhog*diffw*pi*dp*(xw0-xwb)/(1.-xw0)
else
rw = 0.
endif
C-- coal pyrolysis rate
Call perks(ypred,yp,Tpred)
C time step control
if(y(1).gt.5.e-3)then
dy1 = dt*0.5*(yp(1)+ypp(1))
else
dy1 = dt*0.5*(yp(3)+ypp(3))
endif
if(abs(dy1).lt.0.001)then
dt = dt*2.
if(dt.lt.dtmax)print*,'at time=',time,
x ' dt changed to ',dt
elseif(abs(dy1).gt.0.02)then
dt = 0.01/abs(dy1)*dt
print*,'at time=',time, 'dt changed to ',dt
endif
dt = min(dt,dtmax)
C free radical light gas nitrogen release mechanism
if((mw1-mwchar)/mw1.GT.fstable)yp(4)=yp(4)-y(4)*arad*
x exp(-erad/rg/Tpred)*(mwcharold-mwchar)/mwchar*
x machar/mwchar/dt
C COMPONENT MASS CONSERVATION
do 70 j=1,4
y(j) = y(j) + dt*0.5*(yp(j)+ypp(j))
y(j) = max(zero,y(j))
70 continue
C update current and old mwchar
mwcharold = mwchar
g = 2.*(1.-y(1)-c0)-y(2)
mwchar = mw1-g*mdel*sigp1/2.
C crosslinking rate
fracr = 1.
if(fmetold.gt.small.and.acr.gt.0.)then
ratecr = acr*exp(-ecr/rg/tpred)*fmetold*dt
fracr = 1.-ratecr/fmetold
fmet = fmetold-ratecr
fcross = fcrossold+ratecr
if(fmet.lt.0.)then
fcross = fcrossold-ratecr+fmet+ratecr
fmet = 0.0
fracr = 0.
endif
endif
C get product distribution from y array
if(y(1).gt.small)intar = .true.
call perkp(y,mgas,ftar,ftart,fgas,fchart,ft,mt,intar)
intar = .false.
gasmw = rba*ma/2.
C flash distillation
if(fgas.ge.small)then
imolw = .false.
if(mod(iii,iprint).eq.0)imolw=.true.
ipred = .false.
call flash(fgas,gasmw,ft,mt,fracr,ftar,fmet,
x tpred,press,nmax,imolw,ipred)
elseif(fgas.lt.1.e-5)then
fmet = ftart
ftar = 0.
endif
intar = .false.
fvol = fgas+ftar
fchar = 1.-fvol
dvdt = (fvol-fvolold)/dt*alfc0
C done with pyrolysis rate!
C rtot = dvdt + rwp !this looks wrong
rtot = dvdt + rw
delhw= 6.7845E-08*(tpred**4) - 1.1123E-04*(tpred**3)
x + 6.8595E-02*(tpred**2) - 1.8093E+01*tpred + 1.1229E+03
C HEAT = dvdt*delhv + rwp*delhw
HEAT = dvdt*delhv + rw*delhw
C
C calculate heat capacity
fmetmaxc=max(fmetmaxc,fmet)
if (CPSWITCH.ge.0.0) then
C call heatcp(tpred,cpc)
C cpc=(fmet/fmetmaxc)*cpcold+(1.0-fmet/fmetmaxc)*cpc
cpc=cpcold
else
CALL HEATCP(TPRED,CPC)
endif
call heatap(tpred,cpa)
cp = (alfcp*cpc + alfa*cpa + alfw*cpw)/(alpha)
C
if (trate.gt.tratem) trateold=tratem
TRATE = (QCONV+QRAD+HEAT)/(ALPHA*CP)
TRATEM=max(trate,tratem)
Tp = Tp + 0.5*(TRATE+TRATEP)*DT
tpmax=max(tp,tpmax)
C print*,'<trate> qconv,qrad,heat,alpha,cp,trate,tp,tg',
C x qconv,qrad,heat,alpha,cp,trate,tp,tg
Alfw = alfw - rw*dt
ALFV = fvol*alfc0
ALFC = fchar*alfc0
ALFC = MAX(ALFC,0.0)
ALFw = MAX(ALFw,0.0)
ALPHA = (ALFC + alfa + alfw)
omegaa = alfa/alpha
C particle diameter changes due to swelling during devolatilization
volc=fvol-fvolold
if((alfw.eq.0.0).and.(fvol.ge.1.e-7)) then !check for true maximum in heating rate
if ((sratio.lt.0.5).and.(tratem.gt.trate)) then
sratio=smin+svar*sp*((5.79e4/tratem)**cHR)
write(29,*) "Max HR Time ",(time+dt)*1000.
write(29,*)"HR for Swelling Correlation ",tratem
write(29,*)"Correlated Swelling Ratio ",sratio
swell=sratio-1.0 !use swelling correlation from Randy Shurtz
trateold=tratem
elseif ((sc.lt.0.5).and.(trateold.ne.tratem)) then
if (((volc).gt.1.e-6).and.(tratem.gt.trate)) then
sratio=smin+svar*sp*((5.79e4/tratem)**cHR)
write(29,*) "Max HR Time ",(time+dt)*1000.
write(29,*)"HR for Swelling Correlation ",tratem
write(29,*)"Correlated Swelling Ratio ",sratio
swell=sratio-1.0 !update swelling ratio with new maximum heating rate
trateold=tratem
endif
endif
endif
l = y(1)
fnca = y(4)
del2 = y(2)/2.
aind = del2 + l
dp = dp0*(1.+swell*(1.-l/l0))
AP = PI*DP**2
C
C NITROGEN RELEASE CALCULATIONS
C
C calculate tar nitrogen yield
C (assumes tar is released before light gas and HCN)
yntar = yntar + (ftar-ftarold)*fnt
ftarold = ftar
C nitrogen content of char and metaplast
fnt = y(4)*machar/mwchar
C nitrogen remaining in char
ynchar = fchar*fnt
C fraction of original nitrogen remaining in char
fnchar = ynchar/fnit
C fraction of original nitrogen released as tar
fntar = yntar/fnit
C fraction of original nitrogen released as light gas (diff.)
fnhcn = 1.0 - fnchar - fntar
C total fractional release of nitrogen
fntot = (fnit - fnt*fchar)/fnit
C DISTRIBUTE LIGHT GAS INTO H20, CO2, CO, CH4, & other HC's
C yf is a CPD indicator of the fraction of total light gas
C that has been released. The look up table on light gas
C composition is based on yf.
yf = 1-aind/aind0
call lightgas(yygas,inside,lib)
C calculate fraction of total mass release that is h2o, co2, ch4,
C co, and other light gases
do 218 ik=1,5
ffgas(ik)=fgas*yygas(ik)
218 continue
C
C PARTICLE DENSITY CHANGES DURING DEVOLATILIZATION
RHOP = ALPHA/((4./3.)*PI*(DP/2)**3)
C
FRAC = ALFC/ALPHAP
omegawf = alfaw/alpha
w = (alfw/alphap)*100.
TIME = TIME + DT
tms = time*1000.
IF(MOD(I,iprint).EQ.0)THEN
xmm = x*10.
v = fvol*100.
write(1,1000)TMS,Xmm,Tp,TG,v,fchar,fcross,ftar,fmet,trate
write(23,1021)tms,tp,mwchar,y(4),fnt,fnchar,fntar,fnhcn,fntot
write(27,1022)tms,fgas,ffgas(1),ffgas(2),ffgas(3),ffgas(4),
x ffgas(5),yygas(1),yygas(2),yygas(3),yygas(4),
x yygas(5),yf
if(mod(i,iprint).eq.0)then
WRITE(6,1000)TMS,Xmm,Tp,TG,v,fchar,BLOW
endif
ENDIF
1000 FORMAT(' ',F9.3,2X,F9.4,2X,F6.0,5X,F6.0,5X,1PE10.3,5X,1PE10.3
X ,3X,0PF6.3,3X,0PF6.2,3X,0PF6.3,3x,1pe10.3)
1021 format(' ',2(1pe10.3),0pf8.2,6(0pf10.6))
1022 format(' ',1pe10.3,0pf9.4,11(0pf9.4))
IF(TIME.GE.timax)go to 3333
C CHECK TO SEE IF INTERPOLATION FACTORS NEED UPDATE
IF(X.GT.XT(IX+1))THEN
IX = IX + 1
IF(IX.GE.50)go to 3333
ENDIF
IF(X.GT.zv(Iv+1))THEN
Iv = Iv + 1
IF(Iv.GE.50)go to 3333
ENDIF
100 CONTINUE
3333 continue
v = fvol*100.
xmm = x*10.
WRITE(6,1000)TMS,Xmm,Tp,TG,v,fchar,BLOW
write(1,1000)TMS,Xmm,Tp,TG,v,fchar,fcross,
x ftar,fmet,trate
write(29,*)"Final heat capacity (cal/gm-K) ",cpc
write(29,*)"Maximum particle temperature (K) ",tpmax
write(29,*)"Final Volatiles yield (%daf) ",v
write(29,*)"Maximum Heating Rate ", TRATEM
write(29,*)"Swelling Ratio Used ",(swell+1.0)
STOP
END
C
SUBROUTINE PROPS(TG,TP,RHOG,UG,KG,CPG,diffw,press)
C CALCULATION OF GAS PROPERTIES OF N2 AT ATMOSPHERIC PRESSURE; upgraded for pressure by Randy Shurtz 2010
C FROM HOLMAN, P. 505
PARAMETER (NDAT=10)
REAL TG0(NDAT),UG0(NDAT),KG,KG0(NDAT)
DATA TG0/300,400,500,600,700,800,900,1000,1100,1200/
DATA RHO/.3412E-3/,TRHO/1000./
DATA UG0/1.784E-4,2.198E-4,2.57E-4,2.911E-4,3.213E-4,3.484E-4,
X3.749E-4,4.E-4,4.228E-4,4.450E-4/
C UNITS ON KG0 ARE W/M/C, AND MUST BE DIVIDED BY 418.4 TO GET CAL/CM/S/C
DATA KG0/.0262,.03335,.03984,.0458,.0512,.0561,.0607,.0648,.0685,
X.07184/
C HEAT CAPACITY DATA
DATA CPA,CPB,CPC,CPD/6.529,.1488e-2,-0.02771e-5,0./
C CALCULATE GAS DENSITY FROM IDEAL GAS LAW
RHOG = RHO*TRHO*press/TG
C T = FILM TEMPERATURE (TP+TG)/2
T = (TP+TG)/2.
C CALCULATE CPG FROM HEAT CAPACITY EXPRESSION IN Himmelblau
CPG = (CPA+CPB*T+CPC*T**2)/28.
C CALCULATE UG AND KG FROM INTERPOLATION OF TABLE VALUES FROM HOLMAN
C FIND INTERVAL WHERE TEMPERATURE LIES.
IF(T.GT.1200.)THEN
UG = UG0(NDAT)*(T/TG0(NDAT))**.67
KG = KG0(NDAT)*(T/TG0(NDAT))**.58/418.4
RETURN
ENDIF
J=1
DO 10 I=1,9
IF(T.GT.TG0(I))J=J+1
10 CONTINUE
FAC = (TG0(J)-T)/(TG0(J)-TG0(J+1))
UG = -FAC*(UG0(J)-UG0(J+1)) + UG0(J)
KG = (-FAC*(KG0(J)-KG0(J+1)) + KG0(J))/418.4
C p = 1.
mw = 28.
C diffw = 1.905e-5*t**1.67/p
diffw = 4.3991e-7*t**2.334/press !Randy Shurtz 2010 from Bird, Stewart and Lightfoot 2002, H2O and N2
RETURN
END
C this program calculates the heat capacity of a particle from Merrick's
C correlations.
Subroutine heatcp(TP,cp)
C calculates daf coal heat capacity
C tp Particle Temperature (K)
C cp Particle heat capacity (cal/g/K)
C u(i) Atomic Weights of elements
C y(i) Mass fractions of elements (C,H,N,O,S)
C RGAS Gas constant (1.987 cal/gmole/K)
C RGASJ Gas constant (8314.3 J/kg/K)
dimension y(5),u(5)
common/spheat/yelem(5)
data u/12.,1.,14.,16.,32./
data RGAS/1.987/,RGASJ/8314.3/
C calculate mean atomic weight, a
a = 0.0
do 10 i=1,5
a = a+(yelem(i)/u(i))
10 continue
a = 1./a
f1 = 380./tp
f2 = 1800./tp
cp = rgas/a *(g1(f1)+2.*g1(f2))
return
end
C
Function g1(z)
g1=exp(z)/((exp(z)-1)/z)**2
return
end
Subroutine heatap(TP,cpa)
C calculates ash heat capacity
C tp Particle Temperature (K)
C cpa Ash heat capacity (cal/g/K)
C RGAS Gas constant (1.987 cal/g/K)
cpa = (754+.586*(tp-273))/(4.184e3)
return
end
C
C------------------------------------------------------------------
C
subroutine perks(y,yp,T)
IMPLICIT real*4 (A-H,O-Z)
save
C
C this subroutine is the meat of the devolatilization model
C
C y1 = l labile bridges
C y2 = del ends
C y3 = C char links
C y4 = fnca mass nitrogen per mass (aromatic) site
C yp(i) = derivative of y(i) in time
C
C nmax = number of terms in expansion for mol. wt. distribution
real*4 l0,l,ma,kb,kc,kg,kn,kp,ft(35),mt(35),y(4),yp(4)
common/init/l0,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar
common/rate/ab,eb0,ebsig,ac,ec0,ecsig,ag,eg0,egsig,rg
common/nitmod/fnca0,fnca,an,en0,ensig
data fx/0./,ft/35*0./,mt/35*0./
C
l = y(1)
del = y(2)