MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
foamprop.f90
Go to the documentation of this file.
1 
7 module foamprop
8  use constants
9  use ioutils
10  use gasprop
11  use filmprop
12  use cylprop
13  implicit none
14  private
15  public effrad,fbep
16  real(dp) :: unin=1.57_dp
17  real(dp), dimension(:), allocatable :: &
18  lambdaf,&
19  !foam radiative properties - wall contribution
20  kappafwall,sigmafwall,betafwall,omegafwall,betatrfwall,&
21  !foam radiative properties - strut contribution
22  kappafstrut,sigmafstrut,betafstrut,omegafstrut,betatrfstrut,&
23  !foam radiative properties - gas contribution
24  kappafgas,&
25  !foam radiative properties - overall
26  kappaf,sigmaf,betaf,omegaf,betatrf
27 contains
28 !********************************BEGINNING*************************************
30 subroutine effrad(spectra)
31  use quadpack
32  character(len=*), intent(in) :: spectra
33  integer :: i,j,fi,nwawel=100
34 ! real(dp) :: theta !incident angle
35 ! real(dp) :: Rwin !reflectance
36 ! real(dp) :: Twin !transmittance
37 ! real(dp) :: Awin !absorptance
38 ! real(dp) :: rn !random number
39 ! real(dp) :: absp !absorption parameter
40 ! real(dp) :: scap !scattering parameter
41 ! real(dp) :: npart !number of particlesper unit volume
42  real(dp) :: dwmin,dwmax,lambdamin,lambdamax
43  !qags variables
44  real(dp) :: a !start point of integration
45  real(dp) :: abserr
46  real(dp) :: b !end point of integration
47  real(dp), parameter :: epsabs = 0.0e0_dp
48  real(dp) :: epsrel = 1e-3_dp
49  integer :: ier
50  integer :: neval
51  real(dp) :: res
52  !end of qags variables
53  write(*,*) 'Radiative properties of the foam:'
54  write(mfi,*) 'Radiative properties of the foam:'
55 ! call cylconst(2*pi,pi/4,1.0e-1_dp,Qs,Qt,mu)
56 ! stop
57  !monte carlo method, does not work properly
58 ! swin=dcell**2
59 ! absp=0
60 ! scap=0
61 ! npart=3/dcell**3
62 ! lambda=10e-6
63 ! call srand(int(abs(sin(dble(time()))*1e6)))
64 ! do i=1,nrays
65 ! theta=acos(rand())
66 !! theta=rand()*pi/2
67 ! call filmconst(lambda,theta,dwall,Rwin,Twin,Awin)
68 ! absp=absp+Swin*Awin
69 ! scap=scap+Swin*Rwin
70 !! absp=absp+cos(theta)*sin(theta)*Swin*Awin*pi/2
71 !! scap=scap+cos(theta)*sin(theta)*Swin*Rwin*pi/2
72 !! rn=rand()
73 !! if (rn>1-Awin) then
74 !! absp=absp+Swin
75 !! elseif(rn<Rwin) then
76 !! scap=scap+Swin
77 !! endif
78 ! enddo
79 ! kappaf=npart*absp/nrays
80 ! sigmaf=npart*scap/nrays
81 ! betaf=kappaf+sigmaf
82 ! omegaf=sigmaf/betaf
83 ! write(*,*) kappaf
84 ! write(*,*) sigmaf
85 ! write(*,*) betaf
86 ! write(*,*) omegaf
87 ! a=0
88 ! b=pi/2
89 ! call qags ( reflectance, a, b, epsabs, epsrel, res, abserr, neval, ier )
90 ! write(*,*)
91 ! write(*,*) sigmaf
92 ! sigmaf=3*res/dcell
93 ! write(*,*) sigmaf
94 
95  allocate(lambdaf(nwawel),kappafwall(nwawel),sigmafwall(nwawel),&
96  betafwall(nwawel),omegafwall(nwawel),betatrfwall(nwawel),&
97  kappafstrut(nwawel),sigmafstrut(nwawel),betafstrut(nwawel),&
98  omegafstrut(nwawel),betatrfstrut(nwawel),kappafgas(nwawel),&
99  kappaf(nwawel),sigmaf(nwawel),betaf(nwawel),omegaf(nwawel),&
100  betatrf(nwawel))
101  lambdamin=lambdan(1)
102  lambdamax=lambdan(size(lambdan))
103  do i=1,nwawel
104  lambda=lambdamin+(i-1)*(lambdamax-lambdamin)/nwawel
105 ! write(*,*) lambda
106  lambdaf(i)=lambda
107  kappafgas(i)=abscoeffgas(lambda)
108  if (fs<1-struttol) then
109  epsrel=1e-3_dp
110  if (wdist) then
111  dwmax=dwall*10 !if you choose larger interval, it can fail for
112  ! small standard deviation;
113  ! this should work for wsdev between 0.01 and 1
114  dwmin=dwall/10
115  call qag ( scattwall, dwmin, dwmax, epsabs, epsrel, 1, res, &
116  abserr, neval, ier )
117  if (ier /= 0) then
118  write(*,*) 'qag returned',ier
119  write(*,*) 'wavelength',lambda
120  write(*,*) 'scattering coefficient of walls not calculated'
121  write(mfi,*) 'qag returned',ier
122  write(mfi,*) 'wavelength',lambda
123  write(mfi,*) 'scattering coefficient of walls not calculated'
124  stop
125  endif
126  sigmafwall(i)=res
127  call qag ( trextwall, dwmin, dwmax, epsabs, epsrel, 1, res, &
128  abserr, neval, ier )
129  if (ier /= 0) then
130  write(*,*) 'qag returned',ier
131  write(*,*) 'wavelength',lambda
132  write(*,*) 'transport extinction coefficient of walls &
133  not calculated'
134  write(mfi,*) 'qag returned',ier
135  write(mfi,*) 'wavelength',lambda
136  write(mfi,*) 'transport extinction coefficient of walls &
137  not calculated'
138  stop
139  endif
140  betatrfwall(i)=res
141  call qag ( absorwall, dwmin, dwmax, epsabs, epsrel, 1, res, &
142  abserr, neval, ier )
143  if (ier /= 0) then
144  write(*,*) 'qag returned',ier
145  write(*,*) 'wavelength',lambda
146  write(*,*) 'absorption coefficient of walls not calculated'
147  write(mfi,*) 'qag returned',ier
148  write(mfi,*) 'wavelength',lambda
149  write(mfi,*) 'absorption coefficient of walls not calculated'
150  stop
151  endif
152  kappafwall(i)=res
153  else
154  a=0
155  b=pi/2
156  call qags ( scattcoeffintwall, a, b, epsabs, epsrel, res, &
157  abserr, neval, ier )
158  if (ier /= 0) then
159  write(*,*) 'qags returned',ier
160  write(*,*) 'wavelength',lambda
161  write(*,*) 'scattering coefficient of walls not calculated'
162  write(mfi,*) 'qags returned',ier
163  write(mfi,*) 'wavelength',lambda
164  write(mfi,*) 'scattering coefficient of walls not calculated'
165  stop
166  endif
167  sigmafwall(i)=res*(1-por)/dwall
168  call qags ( trextcoeffintwall, a, b, epsabs, epsrel, res, &
169  abserr, neval, ier )
170  if (ier /= 0) then
171  write(*,*) 'qags returned',ier
172  write(*,*) 'wavelength',lambda
173  write(*,*) 'transport extinction coefficient of walls &
174  not calculated'
175  write(mfi,*) 'qags returned',ier
176  write(mfi,*) 'wavelength',lambda
177  write(mfi,*) 'transport extinction coefficient of walls &
178  not calculated'
179  stop
180  endif
181  betatrfwall(i)=res*(1-por)/dwall
182  call qags ( abscoeffintwall, a, b, epsabs, epsrel, res, &
183  abserr, neval, ier )
184  if (ier /= 0) then
185  write(*,*) 'qags returned',ier
186  write(*,*) 'wavelength',lambda
187  write(*,*) 'absorption coefficient of walls not calculated'
188  write(mfi,*) 'qags returned',ier
189  write(mfi,*) 'wavelength',lambda
190  write(mfi,*) 'absorption coefficient of walls not calculated'
191  stop
192  endif
193  kappafwall(i)=res*(1-por)/dwall
194  endif
195  betafwall(i)=kappafwall(i)+sigmafwall(i)
196  omegafwall(i)=sigmafwall(i)/betafwall(i)
197  else
198  kappafwall=0
199  sigmafwall=0
200  betafwall=0
201  omegafwall=0
202  betatrfwall=0
203  endif
204  if (fs>struttol) then
205  a=0
206  b=pi/2
207  if (pi*dstrut/lambda>10) then
208  epsrel=1e-1_dp
209  elseif(pi*dstrut/lambda>1) then
210  epsrel=1e-2_dp
211  else
212  epsrel=1e-3_dp
213  endif
214  call qags ( trextcoeffintstrut, a, b, epsabs, epsrel, res, abserr, &
215  neval, ier )
216  if (ier /= 0) then
217  write(*,*) 'qags returned',ier,neval
218  write(*,*) 'qags returned',res,abserr
219  write(*,*) 'wavelength',lambda
220  write(*,*) 'transport extinction coefficient of struts &
221  not calculated'
222  stop
223  endif
224  betatrfstrut(i)=res*4/pi/dstrut*(1-por)
225  call qags ( extcoeffintstrut, a, b, epsabs, epsrel, res, abserr, &
226  neval, ier )
227  if (ier /= 0) then
228  write(*,*) 'qags returned',ier
229  write(*,*) 'wavelength',lambda
230  write(*,*) 'extinction coefficient of struts not calculated'
231  stop
232  endif
233  betafstrut(i)=res*4/pi/dstrut*(1-por)
234  call qags ( scattcoeffintstrut, a, b, epsabs, epsrel, res, abserr, &
235  neval, ier )
236  if (ier /= 0) then
237  write(*,*) 'qags returned',ier
238  write(*,*) 'wavelength',lambda
239  write(*,*) 'scattering coefficient of struts not calculated'
240  stop
241  endif
242  sigmafstrut(i)=res*4/pi/dstrut*(1-por)
243  kappafstrut(i)=betafstrut(i)-sigmafstrut(i)
244  omegafstrut(i)=sigmafstrut(i)/betafstrut(i)
245  else
246  kappafstrut=0
247  sigmafstrut=0
248  betafstrut=0
249  omegafstrut=0
250  betatrfstrut=0
251  endif
252  enddo
253  kappafwall=(1-fs)*kappafwall
254  sigmafwall=(1-fs)*sigmafwall
255  betafwall=(1-fs)*betafwall
256  betatrfwall=(1-fs)*betatrfwall
257  kappafstrut=fs*kappafstrut
258  sigmafstrut=fs*sigmafstrut
259  betafstrut=fs*betafstrut
260  betatrfstrut=fs*betatrfstrut
261  kappaf=kappafwall+kappafstrut
262  sigmaf=sigmafwall+sigmafstrut
263  betaf=betafwall+betafstrut
264  betatrf=betatrfwall+betatrfstrut
265 ! kappaf=(1-fs)*kappafwall+fs*kappafstrut
266 ! sigmaf=(1-fs)*sigmafwall+fs*sigmafstrut
267 ! betaf=(1-fs)*betafwall+fs*betafstrut
268 ! betatrf=(1-fs)*betatrfwall+fs*betatrfstrut
269  omegaf=sigmaf/betaf
270  open(newunit(fi),file=trim(adjustl(spectra)))
271  write(fi,'(1000A23)') '#wavelength','abs.coeff.wall','scatt.coeff.wall',&
272  'ext.coeff.wall','albedo wall','tr.ext.coeff.wall','abs.coeff.strut',&
273  'scatt.coeff.strut','ext.coeff.strut','albedo strut',&
274  'tr.ext.coeff.strut','abs.coeff','scatt.coeff','ext.coeff','albedo ',&
275  'tr.ext.coeff'
276  do i=1,nwawel
277  write(fi,'(1000es23.15)') lambdaf(i),kappafwall(i),sigmafwall(i),&
278  betafwall(i),omegafwall(i),betatrfwall(i),kappafstrut(i),&
279  sigmafstrut(i),betafstrut(i),omegafstrut(i),betatrfstrut(i),&
280  kappaf(i),sigmaf(i),betaf(i),omegaf(i),betatrf(i)
281  enddo
282  close(fi)
283  !calculate gray properties
284  a=lambdaf(1)
285  b=lambdaf(size(lambdaf))
286  b=1e-4_dp !100 um is relatively enough (if you use lower value, you should
287  ! incorporate fraction of blackbody radiation function)
288  effn=por*n1+(1-por)*unin !this is not perfect, use of some average value
289  ! of n2 instead of unin would be better
290  call qags ( rosextc, a, b, epsabs, epsrel, res, abserr, neval, ier )
291  if (ier /= 0) then
292  write(*,*) 'qags returned',ier
293  write(*,*) 'rosseland extinction coefficient not calculated'
294  stop
295  endif
296  rossextcoeff=(4*effn**2*sigmab*tmean**3)/res
297  call qags ( planckextc, a, b, epsabs, epsrel, res, abserr, neval, ier )
298  if (ier /= 0) then
299  write(*,*) 'qags returned',ier
300  write(*,*) 'planck extinction coefficient not calculated'
301  stop
302  endif
303  planckextcoeff=res/(effn**2*sigmab*tmean**4)
304  call qags ( planckalbedo, a, b, epsabs, epsrel, res, abserr, neval, ier )
305  if (ier /= 0) then
306  write(*,*) 'qags returned',ier
307  write(*,*) 'scattering albedo not calculated'
308  stop
309  endif
310  albedo=res/(effn**2*sigmab*tmean**4)
311  krad=16*sigmab*tmean**3/(3*rossextcoeff)
312  !calculate actual (usually non-gray) properties
313  if (nbox<1) then
314  stop 'choose number of gray boxes 1 for gray approximation, greater &
315  than 1 for non-gray simulation'
316  else
317  allocate(lambdabox(nbox+1),trextcoeffbox(nbox),albedobox(nbox),&
318  abscoeffbox(nbox),scattcoeffbox(nbox))
319  if (.not. allocated(fbepbox)) allocate(fbepbox(nbox))
320  lambdamin=2e-6_dp
321  lambdamax=25e-6_dp
322  lambdabox(nbox+1)=100e-6_dp
323  lambdabox(nbox)=lambdamax
324  lambdabox(1)=lambdamin !if nbox equals 1, then we have only one
325  ! box: 2-100 um
326  do i=2,nbox-1
327  lambdabox(i)=lambdamin+(i-1)*(lambdamax-lambdamin)/(nbox-1._dp)
328  enddo
329  do i=1,nbox
330  fbepbox(i)=fbep(effn,lambdabox(i+1),tmean)-fbep(effn,lambdabox(i)&
331  ,tmean)
332  call qags(planckextc,lambdabox(i),lambdabox(i+1), epsabs, epsrel, &
333  res, abserr, neval, ier)
334  if (ier /= 0) then
335  write(*,*) 'qags returned',ier
336  write(*,*) 'gray box',i
337  write(*,*) 'planck extinction coefficient not calculated'
338  stop
339  endif
340  trextcoeffbox(i)=res/(effn**2*sigmab*tmean**4*fbepbox(i))
341  call qags(planckalbedo,lambdabox(i),lambdabox(i+1), epsabs, &
342  epsrel, res, abserr, neval, ier)
343  if (ier /= 0) then
344  write(*,*) 'qags returned',ier
345  write(*,*) 'gray box',i
346  write(*,*) 'scattering albedo not calculated'
347  stop
348  endif
349  albedobox(i)=res/(effn**2*sigmab*tmean**4*fbepbox(i))
350  enddo
351  scattcoeffbox=albedobox*trextcoeffbox
352  abscoeffbox=trextcoeffbox-scattcoeffbox
353  endif
354  write(*,'(2x,A,1x,es9.3)') 'effective index of refraction:',effn
355  write(*,'(2x,A,1x,es9.3,1x,A)') 'Planck mean extinction coefficient:',&
356  planckextcoeff,'m^-1'
357  write(*,'(2x,A,1x,es9.3,1x,A)') 'Rosseland extinction coefficient:',&
358  rossextcoeff,'m^-1'
359  write(*,'(2x,A,1x,es9.3,1x,A)') 'Rosseland extinction coefficient:',&
360  rossextcoeff/(1-por)/rhos,'kg/m^2'
361  write(*,'(2x,A,1x,es9.3,1x,A)') 'radiative conductivity:',&
362  krad*1e3_dp,'mW/m/K'
363  write(mfi,'(2x,A,1x,es9.3)') 'effective index of refraction:',effn
364  write(mfi,'(2x,A,1x,es9.3,1x,A)') 'Planck mean extinction coefficient:',&
365  planckextcoeff,'m^-1'
366  write(mfi,'(2x,A,1x,es9.3,1x,A)') 'Rosseland extinction coefficient:',&
367  rossextcoeff,'m^-1'
368  write(mfi,'(2x,A,1x,es9.3,1x,A)') 'Rosseland extinction coefficient:',&
369  rossextcoeff/(1-por)/rhos,'kg/m^2'
370  write(mfi,'(2x,A,1x,es9.3,1x,A)') 'radiative conductivity:',&
371  krad*1e3_dp,'mW/m/K'
372  deallocate(lambdaf)
373  deallocate(kappafwall,sigmafwall,betafwall,omegafwall,betatrfwall,&
374  kappafstrut,sigmafstrut,betafstrut,omegafstrut,betatrfstrut,&
375  kappaf,sigmaf,betaf,omegaf,betatrf,kappafgas,lambdabox,trextcoeffbox,&
376  albedobox)
377 end subroutine effrad
378 !***********************************END****************************************
379 
380 
381 
382 !********************************BEGINNING*************************************
384 real(dp) function rosextc(lambda)
385  use interpolation
386  real(dp), intent(in) :: lambda
387  real(dp) :: beta
388  real(dp) :: n !refractive index
389  integer :: ni=1 !number of points, where we want to interpolate
390  real(dp) :: xi(1) !x-values of points, where we want to interpolate
391  real(dp) :: yi(1) !interpolated y-values
392 ! call optconst(lambda,n,k) !could yield recursive call to qags
393  n=effn !used Planck law is valid only constant index of refraction
394  if (lambda<lambdaf(1)) then
395  write(*,*) 'No data for such low wavelength.'
396  write(*,*) 'Minimum wavelength is',lambdaf(1)
397  stop
398  elseif (lambda>lambdaf(size(lambdaf))) then
399 ! write(*,*) 'No data for such high wavelength.'
400 ! write(*,*) 'Maximum wavelength is',lambdaf(size(lambdaf))
401 ! stop
402  beta=betatrf(size(betatrf))
403  else
404  xi(1)=lambda
405  call pwl_interp_1d ( size(lambdaf), lambdaf, betatrf, ni, xi, yi )
406  beta=yi(1)
407  endif
408  rosextc=2*c0**3*exp(hpc*c0/(kb*n*tmean*lambda))*hpc**2*pi/&
409  ((exp(hpc*c0/(kb*n*tmean*lambda))-1)**2*kb*n**3*tmean**2*lambda**6)/beta
410 end function rosextc
411 !***********************************END****************************************
412 
413 
414 !********************************BEGINNING*************************************
416 real(dp) function planckextc(lambda)
417  use interpolation
418  real(dp), intent(in) :: lambda
419  real(dp) :: beta
420  real(dp) :: n !refractive index
421  integer :: ni=1 !number of points, where we want to interpolate
422  real(dp) :: xi(1) !x-values of points, where we want to interpolate
423  real(dp) :: yi(1) !interpolated y-values
424 ! call optconst(lambda,n,k) !could yield recursive call to qags
425  n=effn !used Planck law is valid only constant index of refraction
426  if (lambda<lambdaf(1)) then
427  write(*,*) 'No data for such low wavelength.'
428  write(*,*) 'Minimum wavelength is',lambdaf(1)
429  stop
430  elseif (lambda>lambdaf(size(lambdaf))) then
431 ! write(*,*) 'No data for such high wavelength.'
432 ! write(*,*) 'Maximum wavelength is',lambdaf(size(lambdaf))
433 ! stop
434  beta=betatrf(size(betatrf))
435  else
436  xi(1)=lambda
437  call pwl_interp_1d ( size(lambdaf), lambdaf, betatrf, ni, xi, yi )
438  beta=yi(1)+abscoeffgas(lambda)
439  endif
440  planckextc=2*pi*hpc*c0**2/&
441  (n**2*lambda**5*(exp(hpc*c0/(n*lambda*kb*tmean))-1))*beta
442 end function planckextc
443 !***********************************END****************************************
444 
445 
446 !********************************BEGINNING*************************************
448 real(dp) function planckalbedo(lambda)
449  use interpolation
450  real(dp), intent(in) :: lambda
451  real(dp) :: omega
452  real(dp) :: n !refractive index
453  integer :: ni=1 !number of points, where we want to interpolate
454  real(dp) :: xi(1) !x-values of points, where we want to interpolate
455  real(dp) :: yi(1) !interpolated y-values
456 ! call optconst(lambda,n,k) !could yield recursive call to qags
457  n=effn !used Planck law is valid only constant index of refraction
458  if (lambda<lambdaf(1)) then
459  write(*,*) 'No data for such low wavelength.'
460  write(*,*) 'Minimum wavelength is',lambdaf(1)
461  stop
462  elseif (lambda>lambdaf(size(lambdaf))) then
463 ! write(*,*) 'No data for such high wavelength.'
464 ! write(*,*) 'Maximum wavelength is',lambdaf(size(lambdaf))
465 ! stop
466  omega=omegaf(size(omegaf))
467  else
468  xi(1)=lambda
469  call pwl_interp_1d ( size(lambdaf), lambdaf, omegaf, ni, xi, yi )
470  omega=yi(1)
471  endif
472  planckalbedo=2*pi*hpc*c0**2/&
473  (n**2*lambda**5*(exp(hpc*c0/(n*lambda*kb*tmean))-1))*omega
474 end function planckalbedo
475 !***********************************END****************************************
476 
477 
478 !********************************BEGINNING*************************************
483 real(dp) function fbep(n,lambda,T)
484  real(dp), intent(in) :: n
485  real(dp), intent(in) :: lambda
486  real(dp), intent(in) :: T
487  integer :: i,maxit=100
488  real(dp) :: dzeta,old,res,tol=1e-12_dp
489  dzeta=c2/(n*lambda*t)
490  old=0
491  do i=1,maxit
492  fbep=old+15/pi**4*(exp(-i*dzeta)/i*&
493  (dzeta**3+3*dzeta**2/i+6*dzeta/i**2+6e0_dp/i**3))
494  res=abs(fbep-old)
495  if (res<tol) exit
496  old=fbep
497  if (i==maxit) stop 'fraction of blackbody radiation was &
498  not calculated (did not converge)'
499  enddo
500 end function fbep
501 !***********************************END****************************************
502 end module foamprop
lambda
double lambda
Latent heat of blowing agent, J/kg.
Definition: experimentalInputs.h:118