phys_prop.f90

Go to the documentation of this file.

module phys_prop
    use constants
    use globals
    use foaming_globals_m
    use fmodena
    use modenastuff
    implicit none
    private
    ! interpolation variables
    logical :: rb_initialized
    integer :: rb_kx=4,rb_iknot=0,rb_inbvx
    real(dp), dimension(:), allocatable :: rb_tx,rb_coef
    public set_initial_physical_properties,physical_properties,rb,rderiv,&
        rb_initialized
contains
!********************************BEGINNING*************************************
subroutine set_initial_physical_properties
    time=tstart
    if (.not. firstrun) r0=rb(time)
    radius=r0
    temp=temp0
    conv=0.0_dp
    ! initial bubble contains only air
    xgas=0
    xgas(1)=1
    if (sum(xgas) /= 1) then
        write(*,*) 'Sum of initial molar fractions of gases in the bubble is &
            not equal to one. Normalizing...'
        xgas=xgas/sum(xgas)
        write(*,*) 'New initial molar fractions of gases in the bubble'
        write(*,*) xgas
    endif
    call createmodenamodels
    select case(rhop_model) !density is kept constant, calculate it only once
    case(1)
    case(2)
        call modena_inputs_set(rhopinputs, rhoptpos, temp)
        call modena_inputs_set(rhopinputs, rhopxohpos, 0.1_dp)
        ret = modena_model_call(rhopmodena, rhopinputs, rhopoutputs)
        if(ret /= 0) then
            call exit(ret)
        endif
        rhop = modena_outputs_get(rhopoutputs, 0_c_size_t)
        call modena_inputs_set(rhopinputs, rhoptpos, temp+100)
        call modena_inputs_set(rhopinputs, rhopxohpos, 0.9_dp)
        ret = modena_model_call(rhopmodena, rhopinputs, rhopoutputs)
        if(ret /= 0) then
            call exit(ret)
        endif
        !average density during foaming
        rhop=(rhop + modena_outputs_get(rhopoutputs, 0_c_size_t))/2
    case(3)
        call modena_inputs_set(rhopinputs, rhoptpos, temp)
        ret = modena_model_call(rhopmodena, rhopinputs, rhopoutputs)
        if(ret /= 0) then
            call exit(ret)
        endif
        rhop = modena_outputs_get(rhopoutputs, 0_c_size_t)
        call modena_inputs_set(rhopinputs, rhoptpos, temp+100)
        ret = modena_model_call(rhopmodena, rhopinputs, rhopoutputs)
        if(ret /= 0) then
            call exit(ret)
        endif
        !average density during foaming
        rhop=(rhop + modena_outputs_get(rhopoutputs, 0_c_size_t))/2
    end select
    call physical_properties(temp,conv,radius)
    d0=d
    surface_tension=sigma
    if (geometry=="3D") then
        pair0=(pamb+2*sigma/radius)*xgas(1)
        sn=(1._dp/(nb0*4._dp/3*pi*r0**3)+1)**(1._dp/3)
        s0=sn*radius
        vsh=4*pi/3*(s0**3-r0**3)
    elseif (geometry=="2D") then
        pair0=(pamb+sigma/radius)*xgas(1)
        sn=(1._dp/(nb0*pi*r0**2)+1)**(1._dp/2)
        s0=sn*radius
        vsh=pi*(s0**2-r0**2)
    endif
    gelpoint=.false.
    timestep=(tend-tstart)/its
    if (firstrun) then
        if (allocated(etat)) deallocate(etat)
        if (allocated(port)) deallocate(port)
        if (allocated(init_bub_rad)) deallocate(init_bub_rad)
        allocate(etat(0:its,2),port(0:its,2),init_bub_rad(0:its,2))
    endif
end subroutine set_initial_physical_properties
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine physical_properties(temp,conv,radius)
    real(dp), intent(in) :: temp
    real(dp), intent(in) :: conv
    real(dp), intent(in) :: radius
    integer :: i
    real(dp) :: Aeta,Eeta,Cg,AA,B !viscosity model constants
    if (temp>500) then
        print*, 'temperature > 500', temp
        return
    endif
    if (.not. gelpoint) then
        select case(visc_model)
        case(1)
        case(2)
            aeta=4.1e-8_dp
            eeta=38.3e3_dp
            aa=4.0_dp
            b=-2.0_dp
            eta=aeta*exp(eeta/(rg*temp))*&
                (gelpointconv/(gelpointconv-conv))**(aa+b*conv)
            if (eta>1e10) eta=1.0e10_dp
        case(3)
            !set input vector
            call modena_inputs_set(viscinputs, visctpos, temp);
            call modena_inputs_set(viscinputs, viscxpos, conv);
            !call model
            ret = modena_model_call(viscmodena, viscinputs, viscoutputs)
            if(ret /= 0) then
                call exit(ret)
            endif
            !fetch results
            eta = modena_outputs_get(viscoutputs, 0_c_size_t);
        end select
        if (conv>gelpointconv*0.9_dp) then
            gelpoint=.true.
        endif
    endif
    select case(itens_model)
    case(1)
    case(2)
        call modena_inputs_set(itensinputs, itenstpos, temp)
        ret = modena_model_call(itensmodena, itensinputs, itensoutputs)
        if(ret /= 0) then
            call exit(ret)
        endif
        sigma = modena_outputs_get(itensoutputs, 0_c_size_t)*1e-3_dp
    end select
    do i=1,ngas
        select case(diff_model(i))
        case(2)
            call modena_inputs_set(diffinputs(i), difftpos(i), temp)
            ret = modena_model_call(diffmodena(i),diffinputs(i),diffoutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            d(i) = modena_outputs_get(diffoutputs(i), 0_c_size_t)
        end select
        select case(sol_model(i))
        case(1) !constant Baser 10.1002/pen.760340805
            call modena_inputs_set(solinputs(i), soltpos(i), temp)
            ret = modena_model_call(solmodena(i), solinputs(i), soloutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            kh(i) = modena_outputs_get(soloutputs(i), 0_c_size_t)
        case(2) !pcsaft
            ! TODO: implement properly
            call modena_inputs_set(solinputs(i), soltpos(i), temp)
            call modena_inputs_set(solinputs(i), solxgaspos(i), 1.0e-4_dp)
            call modena_inputs_set(solinputs(i), solxmdipos(i), 0.5_dp)
            call modena_inputs_set(solinputs(i), solxpolyolpos(i), 0.5_dp)
            ret = modena_model_call(solmodena(i), solinputs(i), soloutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            kh(i) = modena_outputs_get(soloutputs(i), 0_c_size_t)
            ! KH(i)=rhop/Mbl(i)/KH(i)
        case(3) !n-pentane, Gupta, 10.1002/pen.11405
            call modena_inputs_set(solinputs(i), soltpos(i), temp)
            ret = modena_model_call(solmodena(i), solinputs(i), soloutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            kh(i) = modena_outputs_get(soloutputs(i), 0_c_size_t)
        case(4) !pentane, Winkler Ph.D.
            call modena_inputs_set(solinputs(i), soltpos(i), temp)
            ret = modena_model_call(solmodena(i), solinputs(i), soloutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            kh(i) = modena_outputs_get(soloutputs(i), 0_c_size_t)
        case(6) !R11, Baser, 10.1002/pen.760340804
            call modena_inputs_set(solinputs(i), soltpos(i), temp)
            ret = modena_model_call(solmodena(i), solinputs(i), soloutputs(i))
            if(ret /= 0) then
                call exit(ret)
            endif
            kh(i) = modena_outputs_get(soloutputs(i), 0_c_size_t)
        case(7) !pentane, Winkler Ph.D.
            kh(i)=rhop/mbl(i)/pamb*(0.0064_dp+0.0551_dp*exp(-(temp-298)**2/&
                (2*17.8_dp**2)))
        case(8) !constant Baser 10.1002/pen.760340805
        end select
    enddo
    if (solcorr) then
        if (geometry=="3D") then
            kh=kh*exp(2*sigma*mbl/(rhop*rg*temp*radius))
        elseif (geometry=="2D") then
            kh=kh*exp(sigma*mbl/(rhop*rg*temp*radius))
        endif
    endif
    cp=cppol+sum(cbl*mbl*cpbll)/rhop
end subroutine physical_properties
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function rderiv(t)
    use bspline_module
    real(dp), intent(in) :: t
    real(dp) :: dt
    integer :: nx,idx,iflag
    ! dt=timestep/100
    ! Rderiv=(Rb(t+dt)-Rb(t))/dt
    nx=size(bub_rad(:,1))
    if (.not. rb_initialized) then
        call rb_spline_ini
    endif
    idx=1
    call db1val(t,idx,rb_tx,nx,rb_kx,rb_coef,rderiv,iflag,rb_inbvx)
    if (iflag /= 0) then
        print*, 'evaluation of bubble radius derivative from spline failed',&
            iflag
        stop
    endif
endfunction rderiv
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function rb(t)
    use bspline_module
    use interpolation
    real(dp), intent(in) :: t
    integer :: nx,idx,iflag
    integer :: ni=1   !number of points, where we want to interpolate
    real(dp) :: xi(1)   !x-values of points, where we want to interpolate
    real(dp) :: yi(1)   !interpolated y-values
    nx=size(bub_rad(:,1))
    ! xi(1)=t
    ! call pwl_interp_1d ( nx, bub_rad(:,1), bub_rad(:,bub_inx+1), ni, xi, yi )
    ! Rb=yi(1)
    if (.not. rb_initialized) then
        call rb_spline_ini
    endif
    idx=0
    call db1val(t,idx,rb_tx,nx,rb_kx,rb_coef,rb,iflag,rb_inbvx)
    if (iflag /= 0) then
        print*, 'evaluation of bubble radius from spline failed',iflag
        stop
    endif
endfunction rb
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine rb_spline_ini
    use bspline_module
    integer :: nx,iflag
    nx=size(bub_rad(:,1))
    if (allocated(rb_tx)) deallocate(rb_tx)
    if (allocated(rb_coef)) deallocate(rb_coef)
    allocate(rb_tx(nx+rb_kx),rb_coef(nx))
    rb_tx=0
    call db1ink(bub_rad(:,1),nx,bub_rad(:,bub_inx+1),&
        rb_kx,rb_iknot,rb_tx,rb_coef,iflag)
    rb_inbvx=1
    rb_initialized=.true.
    if (iflag /= 0) then
        print*, 'initialization of spline failed'
        stop
    endif
end subroutine rb_spline_ini
!***********************************END****************************************
end module phys_prop