!================================================================================================ !================================================================================================= ! Paris-0.1 ! ! Boiling extensions ! written by Leon Malan ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of the ! License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA ! 02111-1307, USA. !================================================================================================= !================================================================================================= Module module_boil use module_grid implicit none !variables real(8), allocatable, dimension(:,:,:) :: Te, dTe, Te_old, s_v, kc, kc_old, Cp, Cp_old real(8), allocatable, dimension(:,:,:) :: energy real(8), allocatable, dimension(:,:,:,:) :: dPg, dug real(8) :: Cp1, Cp2, h_fg, kc1, kc2, T_sat, T0_1, T0_2 real(8) :: mdot, qdot real(8) :: ec1, ec2 ! Thermal expansion coeffs for Boussinesq real(8) :: BC_T(1:6) integer :: BDRY_T(1:6) contains !Read parameters & initialize subroutine ReadHeatParameters use module_VOF use module_flow implicit none integer :: i,j,k,ierr,BC namelist /heatparameters/ Cp1, Cp2, h_fg, kc1, kc2, T_sat, T0_1, T0_2,& BDRY_T, BC_T, ec1, ec2 Cp1 = 1.0d0; Cp2 = 10.0; h_fg=1.0d30; kc1=1.0d0; kc2=1.0d3; T_sat=1.0d0 T0_1=1.0d1; T0_2=1.0d1 BDRY_T = 0 BC_T = 1.0d0 open(unit=10, file='inputheat', status='old', action='read', iostat=ierr) if (ierr .ne. 0) call err_no_out_dir("ReadHeatParameters: error opening 'inputheat' file --- perhaps it does not exist ?") read(10,NML=heatparameters) close(10) do BC=1,6 if (.not.(BDRY_T(BC)==0 .or. BDRY_T(BC)==1)) then write(*,'("Error, Temp BC must be either 0 or 1")') call pariserror("Invalid temperature BC") endif enddo ! Allocate variables allocate (Te(imin:imax,jmin:jmax,kmin:kmax),dTe(imin:imax,jmin:jmax,kmin:kmax),& s_v(imin:imax,jmin:jmax,kmin:kmax),kc(imin:imax,jmin:jmax,kmin:kmax),& Cp(imin:imax,jmin:jmax,kmin:kmax), energy(imin:imax,jmin:jmax,kmin:kmax),& dPg(imin:imax,jmin:jmax,kmin:kmax,1:3),dug(imin:imax,jmin:jmax,kmin:kmax,1:3) ) if (itime_scheme==2) then allocate( Te_old(imin:imax,jmin:jmax,kmin:kmax), kc_old(imin:imax,jmin:jmax,kmin:kmax),& Cp_old(imin:imax,jmin:jmax,kmin:kmax)) endif !Initcondition if (TwoPhase) then if (DoVof) then do i=is,ie; do j=js,je; do k=ks,ke Te(i,j,k) = cvof(i,j,k)*T0_2 + (1.0d0-cvof(i,j,k))*T0_1 enddo; enddo; enddo endif else Te = T0_1 endif end subroutine ReadHeatParameters !get volume source at interface !get heat source at interface !solve energy, conjugate heat transfer !Heat diffusion subroutine TempDiffusion use module_flow implicit none integer :: i,j,k do i=is,ie; do j=js,je; do k=ks,ke dTe(i,j,k) = dTe(i,j,k) + 1.d0/(rho(i,j,k)*cp(i,j,k))*& ( (kc(i+1,j,k)+kc(i,j,k))*(Te(i+1,j,k)-Te(i,j,k)) -& (kc(i,j,k)+kc(i-1,j,k))*(Te(i,j,k)-Te(i-1,j,k)) )/( dxh(i)*dx(i)*2.0d0 ) + & ( (kc(i,j+1,k)+kc(i,j,k))*(Te(i,j+1,k)-Te(i,j,k)) -& (kc(i,j,k)+kc(i,j-1,k))*(Te(i,j,k)-Te(i,j-1,k)) )/( dyh(j)*dy(j)*2.0d0 ) + & ( (kc(i,j,k+1)+kc(i,j,k))*(Te(i,j,k+1)-Te(i,j,k)) -& (kc(i,j,k)+kc(i,j,k-1))*(Te(i,j,k)-Te(i,j,k-1)) )/( dzh(k)*dz(k)*2.0d0 ) enddo; enddo; enddo end subroutine TempDiffusion subroutine get_heat_source implicit none !Calculate temp gradients at interface in normal direction qdot=1.0d2 mdot=qdot/h_fg end subroutine get_heat_source subroutine vofandenergysweeps(tswap) use module_VOF use module_flow use module_tmpvar implicit none integer, dimension(imin:imax,jmin:jmax,kmin:kmax) :: dummy_flag integer, intent(in) :: tswap integer :: i,j,k !get cell energies, energies at s-1 and e+1 needed for BC's do i=is-1,ie+1; do j=js-1,je+1; do k=ks-1,ke+1 energy(i,j,k) = rho(i,j,k)*cp(i,j,k)*Te(i,j,k) enddo; enddo; enddo tmp = cvof dummy_flag = vof_flag if (MOD(tswap,3).eq.0) then ! do z then x then y call swpz_energy(w,tmp,3) call swp(w,tmp,dummy_flag,3,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(u,tmp,1) call swp(u,tmp,dummy_flag,1,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(v,tmp,2) call swp(v,tmp,dummy_flag,2,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) elseif (MOD(tswap,3).eq.1) then ! do y z x call swpz_energy(v,tmp,2) call swp(v,tmp,dummy_flag,2,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(w,tmp,3) call swp(w,tmp,dummy_flag,3,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(u,tmp,1) call swp(u,tmp,dummy_flag,1,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) else ! do x y z call swpz_energy(u,tmp,1) call swpz(u,tmp,dummy_flag,1,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(v,tmp,2) call swpz(v,tmp,dummy_flag,2,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) call swpz_energy(w,tmp,3) call swpz(w,tmp,dummy_flag,3,work(:,:,:,1),work(:,:,:,2),work(:,:,:,3)) endif call get_dT_from_energy(tmp) end subroutine vofandenergysweeps subroutine swpz_energy(us,c,d) use module_VOF use module_flow use module_BC implicit none logical error integer i,j,k integer i0,j0,k0 integer i1,j1,k1 integer, intent(in) :: d real(8), DIMENSION(imin:imax,jmin:jmax,kmin:kmax), intent(in) :: us real(8), DIMENSION(imin:imax,jmin:jmax,kmin:kmax), intent(inout) :: c real(8), DIMENSION(imin:imax,jmin:jmax,kmin:kmax) :: en1,en2,en3 real(8) mm1,mm2,vof real(8) a1,a2,alpha,T_adv real(8) rhocp_ll, rhocp_l, rhocp_r, rhocp_rr, rhocp_c real(8) T_l, T_r,d_int REAL(8) deltax(3),x0(3),fl3d real(8) nr(3), stencil3x3(-1:1,-1:1,-1:1) intrinsic dmax1,dmin1 !*** call init_i0j0k0 (d,i0,j0,k0) if(ng.lt.2) call pariserror("wrong ng") do k=ks-1,ke+1 do j=js-1,je+1 do i=is-1,ie+1 if (d.eq.1) then a2 = us(i,j,k)*dt/dxh(i) a1 = us(i-1,j,k)*dt/dxh(i-1) d_int=dx(i)/2.d0 elseif (d.eq.2) then a2 = us(i,j,k)*dt/dyh(j) a1 = us(i,j-1,k)*dt/dyh(j-1) d_int=dy(j)/2.d0 elseif (d.eq.3) then a2 = us(i,j,k)*dt/dzh(k) a1 = us(i,j,k-1)*dt/dzh(k-1) d_int=dz(k)/2.d0 endif ! Get new density and heat capacity after vofsweep rhocp_ll=(rho2*cp2*c(i-i0*2,j-j0*2,k-k0*2)+rho1*cp1*(1.d0-c(i-i0*2,j-j0*2,k-k0*2))) rhocp_l=(rho2*cp2*c(i-i0,j-j0,k-k0)+rho1*cp1*(1.d0-c(i-i0,j-j0,k-k0))) rhocp_c=(rho2*cp2*c(i ,j ,k )+rho1*cp1*(1.d0-c(i ,j ,k ))) rhocp_r=(rho2*cp2*c(i+i0,j+j0,k+k0)+rho1*cp1*(1.d0-c(i+i0,j+j0,k+k0))) rhocp_rr=(rho2*cp2*c(i+i0*2,j+j0*2,k+k0*2)+rho1*cp1*(1.d0-c(i+i0*2,j+j0*2,k+k0*2))) ! Use new cell temperature after previous energy and vof sweep T_adv = energy(i,j,k)/(rhocp_c) ! Interpolate for T-1/2 (T_l) and T+1/2 (T_r) using chosen AdvectionScheme ! For UpWind, T_l=T_r=T_adv if (.not.(i==is-1 .or. i==ie+1 .or. j==js-1 .or. j==je+1 .or. k==ks-1 .or. k==ke+1)) then if (us(i-i0,j-j0,k-k0)>0.d0) then T_l=interpole3(energy(i-i0*2,j-j0*2,k-k0*2)/(rhocp_ll),energy(i-i0,j-j0,k-k0)/(rhocp_l),& T_adv,AdvectionScheme,d_int) else T_l=interpole3(energy(i+i0,j+j0,k+k0)/(rhocp_r),T_adv,energy(i-i0,j-j0,k-k0)/(rhocp_l),& AdvectionScheme,d_int) endif if (us(i,j,k)>0.d0) then T_r=interpole3(energy(i-i0,j-j0,k-k0)/(rhocp_l),T_adv,energy(i+i0,j+j0,k+k0)/(rhocp_r),& AdvectionScheme,d_int) else T_r=interpole3(energy(i+i0*2,j+j0*2,k+k0*2)/(rhocp_rr),energy(i+i0,j+j0,k+k0)/(rhocp_r),& T_adv,AdvectionScheme,d_int) endif else T_l=T_adv T_r=T_adv endif ! Energy for full cells mm1 = dmax1(a1,0.0d0) mm2 = 1.d0 - mm1 + dmin1(0.d0,a2) if (T_adv > 0.d0) then en1(i,j,k) = dmax1(-a1,0.d0)*energy(i,j,k)/T_adv*T_l en3(i,j,k) = dmax1( a2,0.d0)*energy(i,j,k)/T_adv*T_r en2(i,j,k) = energy(i,j,k) - en1(i,j,k) - en3(i,j,k) else en1(i,j,k) = dmax1(-a1,0.d0)*energy(i,j,k) en3(i,j,k) = dmax1( a2,0.d0)*energy(i,j,k) en2(i,j,k) = energy(i,j,k) - en1(i,j,k) - en3(i,j,k) endif if ((c(i,j,k) .gt. 0.d0).and.(c(i,j,k) .lt. 1.d0)) then do i1=-1,1; do j1=-1,1; do k1=-1,1 stencil3x3(i1,j1,k1) = c(i+i1,j+j1,k+k1) enddo;enddo;enddo call fit_plane_new(c(i,j,k),d,a1,a2,stencil3x3,nr,alpha,error) if(error) then !write(*,*)"WARNING: new plane error!" cycle endif x0=0.d0 deltax=1.d0 if(a1<0.d0) then x0(d)=a1 deltax(d)=-a1 vof = fl3d(nr,alpha,x0,deltax) en1(i,j,k) = ( rho2*cp2*vof+rho1*cp1*(-a1-vof) )*T_l !if (k==3 .and. i==20) write(*,'("Cut cell en1: ",e14.5)')en1(i,j,k) endif if(a2>0.d0) then x0(d)=1.d0 deltax(d)=a2 vof = fl3d(nr,alpha,x0,deltax) en3(i,j,k) = ( rho2*cp2*vof+rho1*cp1*(a2-vof) )*T_r !if (k==3 .and. i==20) then ! write(*,'("Cut cell en3, vof, a2: ",3e14.5)')en3(i,j,k),vof,a2 ! write(*,'("en3 alpha, normals",4e14.5)')alpha,nr(1:3) !endif endif x0(d)=mm1 deltax(d)=mm2 vof = fl3d(nr,alpha,x0,deltax) en2(i,j,k) = ( rho2*cp2*vof+rho1*cp1*(mm2-vof) )*T_adv !if (k==3 .and. i==20) then ! write(*,'("Cut cell en2, vof, a1,a2: ",4e14.5)')en2(i,j,k),vof,a1,a2 ! write(*,'("en2 alpha, normals",4e14.5)')alpha,nr(1:3) !endif endif enddo enddo enddo do k=ks,ke do j=js,je do i=is,ie energy(i,j,k) = en1(i+i0,j+j0,k+k0)+en2(i,j,k)+en3(i-i0,j-j0,k-k0) !if (k==3 .and. i==20 .and. c(i,j,k)<1.d0 .and. c(i,j,k)>0.d0 ) then ! write(*,'("Cut cell energy, en1 right, en2, en3 left: ",4e14.5)')& ! energy(i,j,k),en1(i+i0,j+j0,k+k0),en2(i,j,k),en3(i-i0,j-j0,k-k0) !endif enddo enddo enddo call do_all_ghost(energy) end subroutine swpz_energy subroutine get_dT_from_energy(c) use module_flow implicit none real(8), DIMENSION(imin:imax,jmin:jmax,kmin:kmax), intent(in) :: c integer :: i,j,k do k=ks,ke do j=js,je do i=is,ie dTe(i,j,k) = dTe(i,j,k) + 1.d0/dt*(energy(i,j,k)/& ((rho2*c(i,j,k)+rho1*(1.d0-c(i,j,k)))*(cp2*c(i,j,k)+cp1*(1.d0-c(i,j,k)))) - Te(i,j,k)) enddo enddo enddo end subroutine get_dT_from_energy !Boundary values for temperature subroutine SetTempBc implicit none integer :: i,j,k !BC on x- if (coords(1)==0) then do j=js,je; do k=ks,ke if (BDRY_T(1)==0) then !Dirichlet Te(is-1,j,k)=2.0d0*BC_T(1)-Te(is,j,k) else !Neumann Te(is-1,j,k)=Te(is,j,k)-BC_T(1)*dxh(is-1) endif enddo; enddo endif !BC on x+ if (coords(1)==nPx-1) then do j=js,je; do k=ks,ke if (BDRY_T(4)==0) then Te(ie+1,j,k)=2.0d0*BC_T(4)-Te(ie,j,k) else Te(ie+1,j,k)=Te(ie,j,k)+BC_T(4)*dxh(ie) endif enddo; enddo endif !BC on y- if (coords(2)==0) then do i=is,ie; do k=ks,ke if (BDRY_T(2)==0) then !Dirichlet Te(i,js-1,k)=2.0d0*BC_T(2)-Te(i,js,k) else !Neumann Te(i,js-1,k)=Te(i,js,k)-BC_T(2)*dyh(js-1) endif enddo; enddo endif !BC on y+ if (coords(2)==nPy-1) then do i=is,ie; do k=ks,ke if (BDRY_T(5)==0) then Te(i,je+1,k)=2.0d0*BC_T(5)-Te(i,je,k) else Te(i,je+1,k)=Te(i,je,k)+BC_T(5)*dyh(je) endif enddo; enddo endif !BC on z- if (coords(3)==0) then do i=is,ie; do j=js,je if (BDRY_T(3)==0) then !Dirichlet Te(i,j,ks-1)=2.0d0*BC_T(3)-Te(i,j,ks) else !Neumann Te(i,j,ks-1)=Te(i,j,ks)-BC_T(3)*dzh(ks-1) endif enddo; enddo endif !BC on z+ if (coords(3)==nPz-1) then do i=is,ie; do j=js,je if (BDRY_T(6)==0) then Te(i,j,ke+1)=2.0d0*BC_T(6)-Te(i,j,ke) else Te(i,j,ke+1)=Te(i,j,ke)+BC_T(6)*dzh(ke) endif enddo; enddo endif end subroutine SetTempBC !microlayer end Module module_boil