15 magmit=1000,& !<maximum number of multigrid iterations
21 integer,
dimension(:),
allocatable :: &
24 real(dp),
dimension(:),
allocatable :: &
25 ta,& !<coefficients of temperature sparse matrix
26 trhs,& !<righthandside for temperature matrix
28 integer,
dimension(:,:,:),
allocatable :: &
29 sfiel,& !<structure field (0=solid, 1=pore)
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
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
49 kgas=cond1*por/(por+(1-por)*f)
50 ksol=cond2*f*(1-por)/(por+(1-por)*f)
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,
'%' 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
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' 77 end subroutine effcond
86 subroutine loadstructure_vtk(filename)
88 character(len=80),
intent(in) :: filename
89 character(len=80) :: dummy_char,datatype
91 integer,
dimension(:),
allocatable :: bmat
93 open(newunit(fi),file=filename)
98 read(fi,*) dummy_char,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))
110 if (datatype==
'SCALARS')
then 114 read(fi,*) (sfiel(i,j,k),k=1,dimx)
118 elseif (datatype==
'COLOR_SCALARS')
then 119 allocate(bmat(tnode))
129 bmat(2*i+1)=int(b+0.5)
133 bmat(tnode)=int(a+0.5)
145 stop
'unknown VTK format' 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)
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)
172 end subroutine loadstructure_vtk
178 subroutine savestructure_vtk(filename)
180 character(len=80),
intent(in) :: filename
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' 195 write(fi,
'(1000I2)') (sfiel(i,j,k),k=1,dimx)
200 end subroutine savestructure_vtk
208 subroutine save3dfield_vtk(filename,field)
210 character(len=80),
intent(in) :: filename
211 real(dp),
dimension(:,:,:),
intent(in) :: field
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' 226 write(fi,
'(1000E12.5)') (field(i,j,k),k=1,dimx)
231 end subroutine save3dfield_vtk
237 subroutine initialization_oc_ss
240 allocate(cfiel(dimz,0:dimy+1,0:dimx+1))
244 if (sfiel(i,j,k)==0)
then 254 cfiel(i,j,0)=cfiel(i,j,dimx)
255 cfiel(i,j,dimx+1)=cfiel(i,j,1)
260 cfiel(i,0,j)=cfiel(i,dimy,j)
261 cfiel(i,dimy+1,j)=cfiel(i,1,j)
265 allocate(tfiel(dimz,dimy,dimx))
268 allocate(tvec(dimx*dimy*dimz))
279 allocate(ind(dimz,0:dimy+1,0:dimx+1))
291 ind(i,j,0)=ind(i,j,dimx)
292 ind(i,j,dimx+1)=ind(i,j,1)
297 ind(i,0,j)=ind(i,dimy,j)
298 ind(i,dimy+1,j)=ind(i,1,j)
301 end subroutine initialization_oc_ss
307 subroutine calculation_oc_ss
309 allocate(ta(7*tnode),tja(7*tnode),tia(tnode+1),trhs(tnode))
310 call make_tmtrx_oc_ss
314 call dagmg(tnode,ta(1:tp),tja(1:tp),tia,trhs,tvec,10,1,10,iter,mtol)
316 write(*,*)
'temperature was not calculated' 319 deallocate(ta,tja,tia,trhs,ind)
320 end subroutine calculation_oc_ss
326 subroutine make_tmtrx_oc_ss
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)
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)
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)
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)
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)
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)
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)
409 end subroutine make_tmtrx_oc_ss
415 subroutine make_trhs_oc_ss
438 end subroutine make_trhs_oc_ss
444 subroutine heatflux_oc_ss
447 character(len=80) :: filename
448 allocate(chfz(dimz,dimy,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)
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)
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)
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)
489 xxx=xxx+effc_num/area/(temp1-temp2)*dfoam
492 filename=
'../results/Temp.vtk' 494 deallocate(chfz,tfiel,tvec,cfiel,sfiel,dx,dy,dz)
495 end subroutine heatflux_oc_ss
497 end module conduction