MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
conduction.f90
Go to the documentation of this file.
1 
8 module conduction
9  use constants
10  use ioutils
11  implicit none
12  private
13  public effcond
14  integer :: &
15  magmit=1000,& !<maximum number of multigrid iterations
16  tnode,&
17  tp,&
18  dimz,dimy,dimx
19  real(dp) :: &
20  mtol=1e-30_dp
21  integer, dimension(:), allocatable :: &
22  tja,& !<column index
23  tia
24  real(dp), dimension(:), allocatable :: &
25  ta,& !<coefficients of temperature sparse matrix
26  trhs,& !<righthandside for temperature matrix
27  tvec
28  integer, dimension(:,:,:), allocatable :: &
29  sfiel,& !<structure field (0=solid, 1=pore)
30  ind
31  real(dp), dimension(:,:,:), allocatable :: &
32  dx,dy,dz,& !<voxel sizes in x,y,z
33  cfiel,& !<conductivity field
34  chfz,& !<conductive heat flux in z
35  tfiel
36 contains
37 !********************************BEGINNING*************************************
41 subroutine effcond
42  real(dp) :: xw,xs,f
43  character(len=30) :: fmt='(2x,A,1x,es9.3,1x,A)'
44  write(*,*) 'Conduction - algebraic:'
45  write(mfi,*) 'Conduction - algebraic:'
46  xw=2*(1+cond1/(2*cond2))/3
47  xs=(1+4*cond1/(cond1+cond2))/3
48  f=(1-fs)*xw+fs*xs
49  kgas=cond1*por/(por+(1-por)*f)
50  ksol=cond2*f*(1-por)/(por+(1-por)*f)
51  effc=kgas+ksol
52  eqc_ross=krad+effc
53  gcontr=kgas/eqc_ross
54  scontr=ksol/eqc_ross
55  rcontr=krad/eqc_ross
56  write(*,fmt) 'effective conductivity:',effc*1e3_dp,'mW/m/K'
57  write(*,fmt) 'equivalent conductivity:',eqc_ross*1e3_dp,'mW/m/K'
58  write(*,fmt) 'contribution of gas:',gcontr*1e2_dp,'%'
59  write(*,fmt) 'contribution of solid:',scontr*1e2_dp,'%'
60  write(*,fmt) 'contribution of radiation:',rcontr*1e2_dp,'%'
61  write(mfi,fmt) 'effective conductivity:',effc*1e3_dp,'mW/m/K'
62  write(mfi,fmt) 'equivalent conductivity:',eqc_ross*1e3_dp,'mW/m/K'
63  write(mfi,fmt) 'contribution of gas:',gcontr*1e2_dp,'%'
64  write(mfi,fmt) 'contribution of solid:',scontr*1e2_dp,'%'
65  write(mfi,fmt) 'contribution of radiation:',rcontr*1e2_dp,'%'
66  if (numcond) then
67  write(*,*) 'Conduction - numerical:'
68  write(mfi,*) 'Conduction - numerical:'
69  structurename='../'//structurename
70  call loadstructure_vtk(structurename)
71  call initialization_oc_ss
72  call calculation_oc_ss
73  call heatflux_oc_ss
74  write(*,fmt) 'effective conductivity:',effc_num*1e3_dp,'mW/m/K'
75  write(mfi,fmt) 'effective conductivity:',effc_num*1e3_dp,'mW/m/K'
76  endif
77 end subroutine effcond
78 !***********************************END****************************************
79 
80 
81 !********************************BEGINNING*************************************
86 subroutine loadstructure_vtk(filename)
87  implicit none
88  character(len=80), intent(in) :: filename
89  character(len=80) :: dummy_char,datatype
90  real(dp) :: a,b,c
91  integer, dimension(:), allocatable :: bmat
92  integer :: i,j,k,l,fi
93  open(newunit(fi),file=filename)
94  read(fi,*)
95  read(fi,*)
96  read(fi,*)
97  read(fi,*)
98  read(fi,*) dummy_char,dimz,dimy,dimx
99  tnode=dimz*dimy*dimx
100  allocate(sfiel(dimz+1,0:dimy+1,0:dimx+1),dx(dimz,0:dimy+1,0:dimx+1),&
101  dy(dimz,0:dimy+1,0:dimx+1),dz(dimz,0:dimy+1,0:dimx+1))
102  dz=dfoam/dimz
103  dy=dz
104  dx=dz
105  sfiel=1
106  read(fi,*)
107  read(fi,*)
108  read(fi,*)
109  read(fi,*) datatype
110  if (datatype=='SCALARS') then
111  read(fi,*)
112  do i=1,dimz
113  do j=1,dimy
114  read(fi,*) (sfiel(i,j,k),k=1,dimx)
115  enddo
116  read(fi,*)
117  enddo
118  elseif (datatype=='COLOR_SCALARS') then
119  allocate(bmat(tnode))
120  ! 3 values on first data line
121  read(fi,*) a,b,c
122  bmat(1)=int(a+0.5)
123  bmat(2)=int(b+0.5)
124  bmat(3)=int(c+0.5)
125  do i=1,tnode/2-2
126  ! 2 values on next lines\
127  read(fi,*) a,b
128  bmat(2*i)=int(a+0.5)
129  bmat(2*i+1)=int(b+0.5)
130  enddo
131  ! 1 value on last line
132  read(fi,*) a
133  bmat(tnode)=int(a+0.5)
134  l=1
135  do i=1,dimz
136  do j=1,dimy
137  do k=1,dimx
138  sfiel(i,j,k)=bmat(l)
139  l=l+1
140  enddo
141  enddo
142  enddo
143  deallocate(bmat)
144  else
145  stop 'unknown VTK format'
146  endif
147  close(fi)
148  do i=1,dimz
149  do j=1,dimy
150  sfiel(i,j,0)=sfiel(i,j,dimx)
151  sfiel(i,j,dimx+1)=sfiel(i,j,1)
152  dx(i,j,0)=dx(i,j,dimx)
153  dx(i,j,dimx+1)=dx(i,j,1)
154  dy(i,j,0)=dy(i,j,dimx)
155  dy(i,j,dimx+1)=dy(i,j,1)
156  dz(i,j,0)=dz(i,j,dimx)
157  dz(i,j,dimx+1)=dz(i,j,1)
158  enddo
159  enddo
160  do i=1,dimz
161  do j=1,dimx
162  sfiel(i,0,j)=sfiel(i,dimy,j)
163  sfiel(i,dimy+1,j)=sfiel(i,1,j)
164  dx(i,0,j)=dx(i,dimy,j)
165  dx(i,dimy+1,j)=dx(i,1,j)
166  dy(i,0,j)=dy(i,dimy,j)
167  dy(i,dimy+1,j)=dy(i,1,j)
168  dz(i,0,j)=dz(i,dimy,j)
169  dz(i,dimy+1,j)=dz(i,1,j)
170  enddo
171  enddo
172 end subroutine loadstructure_vtk
173 !***********************************END****************************************
174 
175 
176 !********************************BEGINNING*************************************
178 subroutine savestructure_vtk(filename)
179  implicit none
180  character(len=80), intent(in) :: filename
181  integer :: i,j,k,fi
182  open(newunit(fi),file=filename)
183  write(fi,'(A26)') '# vtk DataFile Version 3.0'
184  write(fi,'(A7)') 'vtkfile'
185  write(fi,'(A5)') 'ASCII'
186  write(fi,'(A25)') 'DATASET STRUCTURED_POINTS'
187  write(fi,'(A10,2x,I8,2x,I4,2x,I4)') 'DIMENSIONS',dimz,dimy,dimx
188  write(fi,'(A12)') 'ORIGIN 0 0 0'
189  write(fi,'(A13)') 'SPACING 1 1 1'
190  write(fi,'(A10,2x,I8)') 'POINT_DATA',tnode
191  write(fi,'(A18)') 'SCALARS values int'
192  write(fi,'(A20)') 'LOOKUP_TABLE default'
193  do i=1,dimz
194  do j=1,dimy
195  write(fi,'(1000I2)') (sfiel(i,j,k),k=1,dimx)
196  enddo
197  write(fi,*)
198  enddo
199  close(fi)
200 end subroutine savestructure_vtk
201 !***********************************END****************************************
202 
203 
204 !********************************BEGINNING*************************************
208 subroutine save3dfield_vtk(filename,field)
209  implicit none
210  character(len=80), intent(in) :: filename
211  real(dp), dimension(:,:,:), intent(in) :: field
212  integer :: i,j,k,fi
213  open(newunit(fi),file=filename)
214  write(fi,'(A26)') '# vtk DataFile Version 3.0'
215  write(fi,'(A7)') 'vtkfile'
216  write(fi,'(A5)') 'ASCII'
217  write(fi,'(A25)') 'DATASET STRUCTURED_POINTS'
218  write(fi,'(A10,2x,I8,2x,I4,2x,I4)') 'DIMENSIONS',dimz,dimy,dimx
219  write(fi,'(A12)') 'ORIGIN 0 0 0'
220  write(fi,'(A13)') 'SPACING 1 1 1'
221  write(fi,'(A10,2x,I8)') 'POINT_DATA',tnode
222  write(fi,'(A20)') 'SCALARS values float'
223  write(fi,'(A20)') 'LOOKUP_TABLE default'
224  do i=1,dimz
225  do j=1,dimy
226  write(fi,'(1000E12.5)') (field(i,j,k),k=1,dimx)
227  enddo
228  write(fi,*)
229  enddo
230  close(fi)
231 end subroutine save3dfield_vtk
232 !***********************************END****************************************
233 
234 
235 !********************************BEGINNING*************************************
237 subroutine initialization_oc_ss
238  integer :: i,j,k,l
239  !assign conductivity to voxels
240  allocate(cfiel(dimz,0:dimy+1,0:dimx+1))
241  do i=1,dimz
242  do j=1,dimy
243  do k=1,dimx
244  if (sfiel(i,j,k)==0) then
245  cfiel(i,j,k)=cond2
246  else
247  cfiel(i,j,k)=cond1
248  endif
249  enddo
250  enddo
251  enddo
252  do i=1,dimz
253  do j=1,dimy
254  cfiel(i,j,0)=cfiel(i,j,dimx)
255  cfiel(i,j,dimx+1)=cfiel(i,j,1)
256  enddo
257  enddo
258  do i=1,dimz
259  do j=1,dimx
260  cfiel(i,0,j)=cfiel(i,dimy,j)
261  cfiel(i,dimy+1,j)=cfiel(i,1,j)
262  enddo
263  enddo
264  !initial temperature field not needed
265  allocate(tfiel(dimz,dimy,dimx))
266  tfiel=0
267  !vector from field
268  allocate(tvec(dimx*dimy*dimz))
269  l=0
270  do i=1,dimz
271  do j=1,dimy
272  do k=1,dimx
273  l=l+1
274  tvec(l)=tfiel(i,j,k)
275  enddo
276  enddo
277  enddo
278  !index
279  allocate(ind(dimz,0:dimy+1,0:dimx+1))
280  l=0
281  do i=1,dimz
282  do j=1,dimy
283  do k=1,dimx
284  l=l+1
285  ind(i,j,k)=l
286  enddo
287  enddo
288  enddo
289  do i=1,dimz
290  do j=1,dimy
291  ind(i,j,0)=ind(i,j,dimx)
292  ind(i,j,dimx+1)=ind(i,j,1)
293  enddo
294  enddo
295  do i=1,dimz
296  do j=1,dimx
297  ind(i,0,j)=ind(i,dimy,j)
298  ind(i,dimy+1,j)=ind(i,1,j)
299  enddo
300  enddo
301 end subroutine initialization_oc_ss
302 !***********************************END****************************************
303 
304 
305 !********************************BEGINNING*************************************
307 subroutine calculation_oc_ss
308  integer :: iter !maximum number of multigrid iterations
309  allocate(ta(7*tnode),tja(7*tnode),tia(tnode+1),trhs(tnode))
310  call make_tmtrx_oc_ss !fill temperature matrix
311  call make_trhs_oc_ss !make right hand side
312  !set maximum number of iterations iter for dagmg (it changes in dagmg)
313  iter=magmit
314  call dagmg(tnode,ta(1:tp),tja(1:tp),tia,trhs,tvec,10,1,10,iter,mtol)
315  if (iter<0) then
316  write(*,*) 'temperature was not calculated'
317  stop
318  endif
319  deallocate(ta,tja,tia,trhs,ind)
320 end subroutine calculation_oc_ss
321 !***********************************END****************************************
322 
323 
324 !********************************BEGINNING*************************************
326 subroutine make_tmtrx_oc_ss
327  integer :: i,j,k,l,m
328  l=1
329  m=1
330  do i=1,dimy
331  do j=1,dimx
332  ta(l)=1
333  tja(l)=ind(1,i,j)
334  tia(m)=l
335  l=l+1
336  m=m+1
337  enddo
338  enddo
339  do i=2,dimz-1
340  do j=1,dimy
341  do k=1,dimx
342  ta(l)=-cfiel(i-1,j,k)*cfiel(i,j,k)*(dz(i-1,j,k)+dz(i,j,k))/&
343  (cfiel(i-1,j,k)*dz(i,j,k)+cfiel(i,j,k)*dz(i-1,j,k))/&
344  (dz(i-1,j,k)/2+dz(i,j,k)/2)*dx(i,j,k)*dy(i,j,k)
345  tja(l)=ind(i-1,j,k)
346  tia(m)=l
347  l=l+1
348  m=m+1
349  ta(l)=-cfiel(i,j-1,k)*cfiel(i,j,k)*(dy(i,j-1,k)+dy(i,j,k))/&
350  (cfiel(i,j-1,k)*dy(i,j,k)+cfiel(i,j,k)*dy(i,j-1,k))/&
351  (dy(i,j-1,k)/2+dy(i,j,k)/2)*dx(i,j,k)*dz(i,j,k)
352  tja(l)=ind(i,j-1,k)
353  l=l+1
354  ta(l)=-cfiel(i,j,k-1)*cfiel(i,j,k)*(dx(i,j,k-1)+dx(i,j,k))/&
355  (cfiel(i,j,k-1)*dx(i,j,k)+cfiel(i,j,k)*dx(i,j,k-1))/&
356  (dx(i,j,k-1)/2+dx(i,j,k)/2)*dy(i,j,k)*dz(i,j,k)
357  tja(l)=ind(i,j,k-1)
358  l=l+1
359  ta(l)=0 &
360  +cfiel(i-1,j,k)*cfiel(i,j,k)*(dz(i-1,j,k)+dz(i,j,k))/&
361  (cfiel(i-1,j,k)*dz(i,j,k)+cfiel(i,j,k)*dz(i-1,j,k))/&
362  (dz(i-1,j,k)/2+dz(i,j,k)/2)*dx(i,j,k)*dy(i,j,k) &
363  +cfiel(i,j-1,k)*cfiel(i,j,k)*(dy(i,j-1,k)+dy(i,j,k))/&
364  (cfiel(i,j-1,k)*dy(i,j,k)+cfiel(i,j,k)*dy(i,j-1,k))/&
365  (dy(i,j-1,k)/2+dy(i,j,k)/2)*dx(i,j,k)*dz(i,j,k) &
366  +cfiel(i,j,k-1)*cfiel(i,j,k)*(dx(i,j,k-1)+dx(i,j,k))/&
367  (cfiel(i,j,k-1)*dx(i,j,k)+cfiel(i,j,k)*dx(i,j,k-1))/&
368  (dx(i,j,k-1)/2+dx(i,j,k)/2)*dy(i,j,k)*dz(i,j,k) &
369  +cfiel(i,j,k)*cfiel(i,j,k+1)*(dx(i,j,k)+dx(i,j,k+1))/&
370  (cfiel(i,j,k)*dx(i,j,k+1)+cfiel(i,j,k+1)*dx(i,j,k))/&
371  (dx(i,j,k)/2+dx(i,j,k+1)/2)*dy(i,j,k)*dz(i,j,k) &
372  +cfiel(i,j,k)*cfiel(i,j+1,k)*(dy(i,j,k)+dy(i,j+1,k))/&
373  (cfiel(i,j,k)*dy(i,j+1,k)+cfiel(i,j+1,k)*dy(i,j,k))/&
374  (dy(i,j,k)/2+dy(i,j+1,k)/2)*dx(i,j,k)*dz(i,j,k) &
375  +cfiel(i,j,k)*cfiel(i+1,j,k)*(dz(i,j,k)+dz(i+1,j,k))/&
376  (cfiel(i,j,k)*dz(i+1,j,k)+cfiel(i+1,j,k)*dz(i,j,k))/&
377  (dz(i,j,k)/2+dz(i+1,j,k)/2)*dx(i,j,k)*dy(i,j,k)
378  tja(l)=ind(i,j,k)
379  l=l+1
380  ta(l)=-cfiel(i,j,k)*cfiel(i,j,k+1)*(dx(i,j,k)+dx(i,j,k+1))/&
381  (cfiel(i,j,k)*dx(i,j,k+1)+cfiel(i,j,k+1)*dx(i,j,k))/&
382  (dx(i,j,k)/2+dx(i,j,k+1)/2)*dy(i,j,k)*dz(i,j,k)
383  tja(l)=ind(i,j,k+1)
384  l=l+1
385  ta(l)=-cfiel(i,j,k)*cfiel(i,j+1,k)*(dy(i,j,k)+dy(i,j+1,k))/&
386  (cfiel(i,j,k)*dy(i,j+1,k)+cfiel(i,j+1,k)*dy(i,j,k))/&
387  (dy(i,j,k)/2+dy(i,j+1,k)/2)*dx(i,j,k)*dz(i,j,k)
388  tja(l)=ind(i,j+1,k)
389  l=l+1
390  ta(l)=-cfiel(i,j,k)*cfiel(i+1,j,k)*(dz(i,j,k)+dz(i+1,j,k))/&
391  (cfiel(i,j,k)*dz(i+1,j,k)+cfiel(i+1,j,k)*dz(i,j,k))/&
392  (dz(i,j,k)/2+dz(i+1,j,k)/2)*dx(i,j,k)*dy(i,j,k)
393  tja(l)=ind(i+1,j,k)
394  l=l+1
395  enddo
396  enddo
397  enddo
398  do i=1,dimy
399  do j=1,dimx
400  ta(l)=1
401  tja(l)=ind(dimz,i,j)
402  tia(m)=l
403  l=l+1
404  m=m+1
405  enddo
406  enddo
407  tia(m)=l
408  tp=l-1
409 end subroutine make_tmtrx_oc_ss
410 !***********************************END****************************************
411 
412 
413 !********************************BEGINNING*************************************
415 subroutine make_trhs_oc_ss
416  integer :: i,j,k,l
417  l=1
418  do i=1,dimy
419  do j=1,dimx
420  trhs(l)=temp2
421  l=l+1
422  enddo
423  enddo
424  do i=2,dimz-1
425  do j=1,dimy
426  do k=1,dimx
427  trhs(l)=0
428  l=l+1
429  enddo
430  enddo
431  enddo
432  do i=1,dimy
433  do j=1,dimx
434  trhs(l)=temp1
435  l=l+1
436  enddo
437  enddo
438 end subroutine make_trhs_oc_ss
439 !***********************************END****************************************
440 
441 
442 !********************************BEGINNING*************************************
444 subroutine heatflux_oc_ss
445  integer :: i,j,k,l
446  real(dp) :: xxx,area
447  character(len=80) :: filename
448  allocate(chfz(dimz,dimy,dimx))
449  l=1
450  do i=1,dimz
451  do j=1,dimy
452  do k=1,dimx
453  tfiel(i,j,k)=tvec(l)
454  l=l+1
455  enddo
456  enddo
457  enddo
458  do i=2,dimz-1
459  do j=1,dimy
460  do k=1,dimx
461  if (sfiel(i,j,k)==sfiel(i+1,j,k)) then
462  chfz(i,j,k)=-(tfiel(i,j,k)-tfiel(i+1,j,k))/&
463  (dz(i,j,k)/2+dz(i+1,j,k)/2)*cfiel(i,j,k)
464  else
465  chfz(i,j,k)=(tfiel(i,j,k)-tfiel(i-1,j,k))/&
466  (dz(i,j,k)/2+dz(i-1,j,k)/2)*cfiel(i,j,k)
467  endif
468  enddo
469  enddo
470  enddo
471  do j=1,dimy
472  do k=1,dimx
473  chfz(1,j,k)=-(tfiel(1,j,k)-tfiel(2,j,k))/&
474  (dz(1,j,k)/2+dz(2,j,k)/2)*cfiel(1,j,k)
475  chfz(dimz,j,k)=(tfiel(dimz,j,k)-tfiel(dimz-1,j,k))/&
476  (dz(dimz,j,k)/2+dz(dimz-1,j,k)/2)*cfiel(dimz,j,k)
477  enddo
478  enddo
479  xxx=0
480  do i=1,dimz
481  effc_num=0
482  area=0
483  do j=1,dimy
484  do k=1,dimx
485  effc_num=effc_num+chfz(i,j,k)*dx(i,j,k)*dy(i,j,k)
486  area=area+dx(i,j,k)*dy(i,j,k)
487  enddo
488  enddo
489  xxx=xxx+effc_num/area/(temp1-temp2)*dfoam
490  enddo
491  effc_num=xxx/dimz
492  filename='../results/Temp.vtk'
493  ! call save3DField_vtk(filename,tfiel)
494  deallocate(chfz,tfiel,tvec,cfiel,sfiel,dx,dy,dz)
495 end subroutine heatflux_oc_ss
496 !***********************************END****************************************
497 end module conduction